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

[Olivier Mattelaer]

Hi guys,

With Rikkert and Stefano, we agreed to release 2.5.6 this friday.
They are two branch that need to be merged before that:
1) bias_function ( which i need to finish the review --already started--)
2) reweight_mass ( Which both Valentin and Rikkert have started to review)

Two branch have actually already been approved for merging but (as far as I know where not merged already):
HwU_allocated -> Rikkert can you merge it (if you want)
mg5hpc -> That's for me.

I would like to see if it is possible to merge:
maddm_dev (would be review by Valentin)
but if that does not make it, it is not critical.

Anything else that you want to include in this version?
I would actually propose to call it 2.6.0 since this version includes a lot of new features.

Cheers,

Olivier

[Stefano Frixione (stefano-frixione)]

Hi Olivier,

> I would actually propose to call it 2.6.0 since this version includes a
> lot of new features.
I think we have advertised a few times that 2.6.0 would contain
EW corrections, so I personally prefer 2.5.6 (not a vry strong
feeling; I'll be happy to go with the majority).

Concerning the review: after the MS bug, I've kept on testing the NLO
bias stuff with cases more involved than single V, and I didn't
find any issues.

Cheers, Stefano.

>
> Cheers,
>
> Olivier
> --
> The attached diff has been truncated due to its size.
> Your team MadTeam is requested to review the proposed merge of lp:~mg5core1/mg5amcnlo/2.5.6 into lp:mg5amcnlo.
>

[Rikkert Frederix (frederix)]

Hi Olivier, Stefano,

I think EW should be 3.0.0.

I don't mind calling it 2.6.0, but 2.5.6 is also fine for me. No strong feelings either way.

HwU_allocated has already been merged a couple of versions ago.

I'm currently testing the reweight_mass with the complete H+2j FxFx. I've got ~500k events. It seems to be a bit faster as the old hand-written reweighting. I'll make some plots this afternoon and check with our previous results.

Cheers,
Rik

[lp:~mg5core1/mg5amcnlo/2.5.6]

317. By Rikkert Frederix on 2017-08-07

removed the s-hat as dynamical-scale-choice option 4. This is not IR safe.

318. By Rikkert Frederix on 2017-08-08

merge with the bias_function branch

319. By Rikkert Frederix on 2017-08-08

Fixed a problem with the bias in split-event generation and the cross
section printed to the screen at the end of the run.

320. By olivier-mattelaer on 2017-08-09

merge reweight_mass

321. By olivier-mattelaer on 2017-08-10

allow syntax set nevents 100k and 1M (in general for all integer variable)

322. By marco zaro on 2017-08-11

fix in combine_root.C; variables were declared inside loops but
used also outside

323. By olivier-mattelaer on 2017-08-14

update iotest

324. By olivier-mattelaer on 2017-08-14

fix to put support for lhapdf6.2

325. By olivier-mattelaer on 2017-08-14

merge maddm_dev (v.298)

326. By olivier-mattelaer on 2017-08-14

merge 2.5.6_cflow_fix

327. By olivier-mattelaer on 2017-08-14

update Update note and documentation

328. By olivier-mattelaer on 2017-08-15

fixing tests

329. By olivier-mattelaer on 2017-08-15

do the reweigthing always in the center of mass frame to avoid loop stability issue

330. By olivier-mattelaer on 2017-08-15

fixing test issue + one bug in the re-weighting in output mode 2.0

331. By Valentin Hirschi on 2017-08-16

1. Fixed issue in the PY8 histogramming related to
         https://bugs.launchpad.net/mg5amcnlo/+bug/1706139
2. Updated ninja tarball.

332. By olivier-mattelaer on 2017-08-16

fixing one parralel test

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2017-05-22 08:05:00 +0000
+++ MadSpin/decay.py	2017-08-16 21:26:30 +0000
@@ -207,6 +207,7 @@
         line+="\n"
         return line
 
+
     def reshuffle_resonances(self,mother):
         """ reset the momentum of each resonance in the production event
                 to the sum of momenta of the daughters
@@ -222,9 +223,9 @@
                 if self.resonance[index]["mothup1"]==mother:
                     daughters.append(index)
 
-        if len(daughters)!=2:
-            logger.info("Got more than 2 (%s) daughters for one particles" % len(daughters))
-            logger.info("in one production event (before decay)")
+#        if len(daughters)!=2:
+#            logger.info("Got more than 2 (%s) daughters for one particles" % len(daughters))
+#            logger.info("in one production event (before decay)")
 
 
         if daughters[0]>0:
@@ -1633,11 +1634,15 @@
                 self.ask_edit_cards(['param_card'],[], plot=False)
         
         
+        commandline = 'import model %s' % model.get('modelpath+restriction') 
+        if not model.mg5_name:
+            commandline += ' --modelname'
+        cmd.exec_cmd(commandline)
+
         line = 'compute_widths %s %s' % \
                 (' '.join([str(i) for i in opts['particles']]),
                  ' '.join('--%s=%s' % (key,value) for (key,value) in opts.items()
                         if key not in ['model', 'force', 'particles'] and value))
-        cmd.exec_cmd('import model %s' % model.get('modelpath+restriction'))
 
         #pattern for checking complex mass scheme.
         has_cms = re.compile(r'''set\s+complex_mass_scheme\s*(True|T|1|true|$|;)''', re.M)
 
=== modified file 'MadSpin/src/driver.f'
--- MadSpin/src/driver.f	2017-03-13 14:10:42 +0000
+++ MadSpin/src/driver.f	2017-08-16 21:26:30 +0000
@@ -384,6 +384,7 @@
           do k=1,nb_mc_masses
              m(indices_mc_masses(k))=values_mc_masses(k)
           enddo
+          ivar=0
           call generate_momenta_conf(jac,x,itree,qmass,qwidth,ptrial,pprod,map_external2res) 
           
           if (jac.lt.0d0) then
 
=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2016-11-21 14:12:27 +0000
+++ Template/LO/Cards/run_card.dat	2017-08-16 21:26:30 +0000
@@ -114,7 +114,6 @@
  %(pta)s  = pta       ! minimum pt for the photons 
  %(ptl)s  = ptl       ! minimum pt for the charged leptons 
  %(misset)s  = misset    ! minimum missing Et (sum of neutrino's momenta)
- %(ptheavy)s  = ptheavy   ! minimum pt for one heavy final state
  %(ptjmax)s  = ptjmax    ! maximum pt for the jets
  %(ptbmax)s  = ptbmax    ! maximum pt for the b
  %(ptamax)s  = ptamax    ! maximum pt for the photons
@@ -191,6 +190,7 @@
 #*********************************************************************
 # Inclusive cuts                                                     *
 #*********************************************************************
+ %(ptheavy)s  = ptheavy   ! minimum pt for at least one heavy final state
  %(xptj)s  = xptj ! minimum pt for at least one jet  
  %(xptb)s  = xptb ! minimum pt for at least one b 
  %(xpta)s  = xpta ! minimum pt for at least one photon 
 
=== modified file 'Template/LO/Source/.make_opts'
--- Template/LO/Source/.make_opts	2016-12-02 21:56:34 +0000
+++ Template/LO/Source/.make_opts	2017-08-16 21:26:30 +0000
@@ -37,9 +37,6 @@
 CFLAGS= -O $(STDLIB_FLAG) $(MACFLAG)
 endif
 
-# Increase the number of allowed charcters in a Fortran line
-FFLAGS+= -ffixed-line-length-132
-
 # Set FC unless it's defined by an environment variable
 ifeq ($(origin FC),default)
 FC=$(DEFAULT_F_COMPILER)
@@ -48,6 +45,15 @@
 F2PY=$(DEFAULT_F2PY_COMPILER)
 endif
 
+# Increase the number of allowed charcters in a Fortran line
+ifeq ($(FC), ftn)
+FFLAGS+= -extend-source # for ifort type of compiler
+else
+FFLAGS+= -ffixed-line-length-132
+endif
+
+
+
 UNAME := $(shell uname -s)
 ifeq ($(origin LDFLAGS), undefined)
 LDFLAGS=$(STDLIB) $(MACFLAG)
@@ -89,12 +95,20 @@
 # Option lhapdf
 
 ifneq ($(lhapdf),)
-CXXFLAGS += $(shell $(lhapdf) --cppflags)
-alfas_functions=alfas_functions_lhapdf
-llhapdf+= -lLHAPDF
+  CXXFLAGS += $(shell $(lhapdf) --cppflags)
+  alfas_functions=alfas_functions_lhapdf
+  llhapdf+= $(shell $(lhapdf) --cflags --libs) -lLHAPDF
+# check if we need to activate c++11 (for lhapdf6.2)
+  ifeq ($(origin CXX),default)
+    ifeq ($lhapdfversion$lhapdfsubversion,62)
+      CXX=$(DEFAULT_CPP_COMPILER) -std=c++11
+    else
+      CXX=$(DEFAULT_CPP_COMPILER)
+    endif
+  endif	     
 else
-alfas_functions=alfas_functions
-llhapdf=
+  alfas_functions=alfas_functions
+  llhapdf=
 endif
 
 # Helper function to check MG5 version
 
=== modified file 'Template/LO/Source/PDF/pdf.f'
--- Template/LO/Source/PDF/pdf.f	2014-10-20 07:53:19 +0000
+++ Template/LO/Source/PDF/pdf.f	2017-08-16 21:26:30 +0000
@@ -14,7 +14,7 @@
       include 'pdf.inc'
 C      
       call fdist(ih,x, q, pdf)
-      
+
       return	
       end
 
 
=== modified file 'Template/LO/Source/PDF/pdg2pdf.f'
--- Template/LO/Source/PDF/pdg2pdf.f	2014-01-07 10:50:42 +0000
+++ Template/LO/Source/PDF/pdg2pdf.f	2017-08-16 21:26:30 +0000
@@ -26,6 +26,11 @@
       data pdlabellast/2*'abcdefg'/
       data ihlast/2*-99/
 
+      if (ih.eq.9) then
+         pdg2pdf = 1d0
+         return
+      endif
+
 c     Make sure we have a reasonable Bjorken x. Note that even though
 c     x=0 is not reasonable, we prefer to simply return pdg2pdf=0
 c     instead of stopping the code, as this might accidentally happen.
 
=== modified file 'Template/LO/Source/PDF/pdg2pdf_lhapdf.f'
--- Template/LO/Source/PDF/pdg2pdf_lhapdf.f	2013-12-09 11:12:33 +0000
+++ Template/LO/Source/PDF/pdg2pdf_lhapdf.f	2017-08-16 21:26:30 +0000
@@ -22,6 +22,11 @@
       data pdflast/30*-99d9/
       data imemlast/2*-99/
 
+      if (ih.eq.9) then
+         pdg2pdf =1d0
+         return
+      endif
+
 c     Make sure we have a reasonable Bjorken x. Note that even though
 c     x=0 is not reasonable, we prefer to simply return pdg2pdf=0
 c     instead of stopping the code, as this might accidentally happen.
 
=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2015-10-01 10:08:45 +0000
+++ Template/LO/Source/dsample.f	2017-08-16 21:26:30 +0000
@@ -997,6 +997,43 @@
 
       end subroutine write_discrete_grids
 
+      subroutine write_grid(name)
+c************************************************************************
+c     Write out the grid
+c************************************************************************
+      implicit none
+
+      character*(*) name
+
+      include 'genps.inc'
+
+      double precision tmean, trmean, tsigma
+      integer             dim, events, itm, kn, cur_it, invar, configs
+      common /sample_common/
+     .     tmean, trmean, tsigma, dim, events, itm, kn, cur_it, invar, configs
+
+      double precision twgt, maxwgt,swgt(maxevents)
+      integer                             lun, nw, itmin
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
+
+      double precision    grid(2, ng, 0:maxinvar)
+      common /data_grid/ grid
+
+      double precision force_max_wgt
+      common/unwgt_secondary_max/force_max_wgt
+
+      integer i,j
+
+      open(26, file='ftn26',status='unknown')
+      write(26,fmt='(4f21.17)') ((grid(2,i,j),i=1,ng),j=1,invar)
+      write(26,*) twgt, force_max_wgt
+c                  write(26,fmt='(4f21.16)') (alpha(i),i=1,maxconfigs)
+      call write_discrete_grids(26,'ref')
+      close(26)
+      return
+      end
+
+
       subroutine read_discrete_grids(stream_id)
 c************************************************************************
 c     Write out the grid using the DiscreteSampler module
@@ -1714,7 +1751,7 @@
          if (kn .ge. max_events .and. non_zero .le. 5) then
             call none_pass(max_events)
          endif
-         if (non_zero .eq. events .or. (kn .gt. 200*events .and.
+         if (non_zero .ge. events .or. (kn .gt. 200*events .and.
      $        non_zero .gt. 5)) then
 
 c          # special mode where we store information to combine them
@@ -2179,13 +2216,7 @@
      .              /23X,11HCross sec =,e12.4/
      .              13X,21HChi**2 per DoF.     =,f12.4/1X,79(1H-))
                if (use_cut .ne. 0) then
-               open(26, file='ftn26',status='unknown')
-               write(26,fmt='(4f21.17)')
-     $              ((grid(2,i,j),i=1,ng),j=1,invar)
-               write(26,*) twgt, force_max_wgt
-               call write_discrete_grids(26,'ref')
-c               write(26,fmt='(4f21.17)') (alpha(i),i=1,maxconfigs)
-               close(26)
+                   call write_grid('ftn26')
                endif
                call sample_writehtm()
 c               open(unit=22,file=result_file,status='old',
 
=== added file 'Template/LO/SubProcesses/dummy_fct.f'
--- Template/LO/SubProcesses/dummy_fct.f	1970-01-01 00:00:00 +0000
+++ Template/LO/SubProcesses/dummy_fct.f	2017-08-16 21:26:30 +0000
@@ -0,0 +1,47 @@
+      subroutine get_dummy_x1(sjac, X1, R, pbeam1, pbeam2, stot, shat)
+      implicit none
+      include 'maxparticles.inc'
+      include 'run.inc'
+c      include 'genps.inc'
+      double precision sjac ! jacobian. should be updated not reinit
+      double precision X1   ! bjorken X. output
+      double precision R    ! random value after grid transfrormation. between 0 and 1
+      double precision pbeam1(0:3) ! momentum of the first beam (input and/or output)
+      double precision pbeam2(0:3) ! momentum of the second beam (input and/or output)
+      double precision stot        ! total energy  (input and /or output)
+      double precision shat        ! output
+
+c     global variable to set (or not)
+      double precision cm_rap
+      logical set_cm_rap
+      common/to_cm_rap/set_cm_rap,cm_rap
+      
+      set_cm_rap=.false. ! then cm_rap will be set as .5d0*dlog(xbk(1)*ebeam(1)/(xbk(2)*ebeam(2)))
+                         ! ebeam(1) and ebeam(2) are defined here thanks to 'run.inc'
+      shat = x1*ebeam(1)*ebeam(2)
+      return 
+      end
+
+      subroutine get_dummy_x1_x2(sjac, X, R, pbeam1, pbeam2, stot,shat)
+      implicit none
+      include 'maxparticles.inc'
+      include 'run.inc'
+c      include 'genps.inc'
+      double precision sjac ! jacobian. should be updated not reinit
+      double precision X(2)   ! bjorken X. output
+      double precision R(2)    ! random value after grid transfrormation. between 0 and 1
+      double precision pbeam1(0:3) ! momentum of the first beam
+      double precision pbeam2(0:3) ! momentum of the second beam
+      double precision stot        ! total energy
+      double precision shat        ! output
+
+c     global variable to set (or not)
+      double precision cm_rap
+      logical set_cm_rap
+      common/to_cm_rap/set_cm_rap,cm_rap
+      
+      set_cm_rap=.false. ! then cm_rap will be set as .5d0*dlog(xbk(1)*ebeam(1)/(xbk(2)*ebeam(2)))
+                         ! ebeam(1) and ebeam(2) are defined here thanks to 'run.inc'
+      shat = x(1)*x(2)*ebeam(1)*ebeam(2)
+      return 
+      end
 
=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2016-09-06 15:01:18 +0000
+++ Template/LO/SubProcesses/genps.f	2017-08-16 21:26:30 +0000
@@ -124,7 +124,6 @@
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       common/to_sprop/sprop,tprid
-      integer          lwgt(0:maxconfigs,maxinvar)
       logical firsttime
 
       double precision xprop(3,nexternal),tprop(3,nexternal)
@@ -132,7 +131,7 @@
       integer imatch
       save maxwgt
 
-      integer ninvar, nconfigs
+      integer ninvar
       
 c
 c     External
@@ -142,6 +141,9 @@
 c
 c     Global
 c
+      integer          lwgt(0:maxconfigs,maxinvar)
+      common/to_lwgt/lwgt
+
       double precision pmass(nexternal)
       common/to_mass/  pmass
 
@@ -172,7 +174,6 @@
       double precision stot,m1,m2
       common/to_stot/stot,m1,m2
 
-      save lwgt
       save ndim,nfinal,nbranch,nparticles
 
       integer jfig,k
@@ -200,40 +201,7 @@
 c      write(*,*) 'using iconfig',iconfig
       if (firsttime) then
          firsttime=.false.
-         do i=1,nexternal
-            m(i)=pmass(i)
-         enddo
-c        Set stot
-         if (nincoming.eq.1) then
-            stot=m(1)**2
-         else
-            m1=m(1)
-            m2=m(2)
-            if (abs(lpp(1)) .eq. 1 .or. abs(lpp(1)) .eq. 2) m1 = 0.938d0
-            if (abs(lpp(2)) .eq. 1 .or. abs(lpp(2)) .eq. 2) m2 = 0.938d0
-            if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
-            if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
-            if(ebeam(1).lt.m1) ebeam(1)=m1
-            if(ebeam(2).lt.m2) ebeam(2)=m2
-            pi1(0)=ebeam(1)
-            pi1(3)=sqrt(max(ebeam(1)**2-m1**2, 0d0))
-            pi2(0)=ebeam(2)
-            pi2(3)=-sqrt(max(ebeam(2)**2-m2**2, 0d0))
-            stot=m1**2+m2**2+2*(pi1(0)*pi2(0)-pi1(3)*pi2(3))
-         endif
-         write(*,'(x,a,f13.2)') 'Set CM energy to ',sqrt(stot)
-c        Start graph mapping
-         do i=1,mapconfig(0)
-            if (mapconfig(i) .eq. iconfig) this_config=i
-         enddo
-         write(*,*) 'Mapping Graph',iconfig,' to config',this_config
-         iconfig = this_config
-         nconfigs = 1
-         mincfig=iconfig
-         maxcfig=iconfig
-         call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
-         maxwgt=0d0
-c         write(*,'(a,12i4)') 'Summing configs',(isym(i),i=1,isym(0))
+         call configure_integral(this_config,mincfig, maxcfig, invar,maxwgt)
          nparticles   = nexternal
          nfinal       = nparticles-nincoming
          nbranch      = nparticles-2
@@ -241,51 +209,10 @@
          if (ndim .lt. 0) ndim = 0   !For 2->1 processes  tjs 5/24/2010
          if (abs(lpp(1)) .ge. 1) ndim=ndim+1
          if (abs(lpp(2)) .ge. 1) ndim=ndim+1
-         call set_peaks
-         if (.false. ) then
-            call find_matches(iconfig,isym(0))
-            write(*,'(a,12i4)') 'Summing configs',(isym(i),i=1,isym(0))
-         endif
-         if (.false.) then
-            i=1
-            do while (mapconfig(i) .ne. iconfig
-     $          .and. i .lt. mapconfig(0))
-               i=i+1
-            enddo
-         endif
-
+         do i=1,nexternal
+            m(i)=pmass(i)
+         enddo
          write(*,'(a,12e10.3)') ' Masses:',(m(i),i=1,nparticles)
-         do j=1,invar
-            lwgt(0,j)=0
-         enddo
-c
-c     Here we set up which diagrams contribute to each variable
-c     in principle more than 1 diagram can contribute to a variable
-c     if we believe they will have identical structure.
-c
-c         do i=1,mapconfig(0)
-         do i=mincfig,maxcfig
-c         do k=1,isym(0)
-c            i = isym(k)
-            write(*,'(15i4)') i,(minvar(j,i),j=1,ndim)
-            do j=1,ndim
-               ipole = minvar(j,i)
-               if (ipole .ne. 0) then
-                  n = lwgt(0,ipole)+1
-                  lwgt(n,ipole)=mapconfig(i)  
-                  lwgt(0,ipole)=n
-               endif
-            enddo
-         enddo
-
-c     Initialize dsig (needed for subprocess group running mode)
-         dum=dsig(0,0,1)
-
-      else
-         do i=1,11
-c            swidth(i)=-5d0         !tells us to use the same point over again
-         enddo
-c         swidth(10)=0d0
       endif                          !First_time
 
       if (.false.) then
@@ -301,24 +228,48 @@
       xbk(2)   = 1d0
       sjac = 1d0
       if (abs(lpp(1)) .ge. 1 .and. abs(lpp(2)) .ge. 1) then
-         call sample_get_x(sjac,x(ndim-1),ndim-1,mincfig,0d0,1d0)
+         if (abs(lpp(1)).eq.9.or.abs(lpp(2)).eq.9)then
+            call sample_get_x(sjac,x(ndim),ndim,mincfig,0d0,1d0)
+            call sample_get_x(sjac,x(ndim-1),ndim-1,mincfig,0d0,1d0)
+            call get_dummy_x1_x2(sjac, Xbk(1), x(ndim-1),pi1, pi2, stot, s(-nbranch))
+            if (.not.set_cm_rap)then
+               cm_rap=.5d0*dlog(xbk(1)*ebeam(1)/(xbk(2)*ebeam(2)))
+               set_cm_rap=.true.
+            endif
+         else
+            call sample_get_x(sjac,x(ndim-1),ndim-1,mincfig,0d0,1d0)
 c-----
 c tjs 5/24/2010 for 2->1 process
 c-------
-         xtau = x(ndim-1)
-         if(nexternal .eq. 3) then
-            x(ndim-1) = pmass(3)*pmass(3)/stot
-            sjac=1 / stot     !for delta function in d_tau
+            xtau = x(ndim-1)
+            if(nexternal .eq. 3) then
+               x(ndim-1) = pmass(3)*pmass(3)/stot
+               sjac=1 / stot    !for delta function in d_tau
+            endif
+
+            call sample_get_x(sjac,x(ndim),ndim,mincfig,0d0,1d0)
+            CALL GENCMS(STOT,Xbk(1),Xbk(2),X(ndim-1), SMIN,SJAC)
+            x(ndim-1) = xtau    !Fix for 2->1 process
+c           Set CM rapidity for use in the rap() function
+            cm_rap=.5d0*dlog(xbk(1)*ebeam(1)/(xbk(2)*ebeam(2)))
+            set_cm_rap=.true.
+c           Set shat
+            s(-nbranch) = xbk(1)*xbk(2)*stot
          endif
 
+      elseif (lpp(1).eq.9.or.lpp(2).eq.9) then 
          call sample_get_x(sjac,x(ndim),ndim,mincfig,0d0,1d0)
-         CALL GENCMS(STOT,Xbk(1),Xbk(2),X(ndim-1), SMIN,SJAC)
-         x(ndim-1) = xtau                   !Fix for 2->1 process
-c        Set CM rapidity for use in the rap() function
-         cm_rap=.5d0*dlog(xbk(1)*ebeam(1)/(xbk(2)*ebeam(2)))
-         set_cm_rap=.true.
-c        Set shat
-         s(-nbranch) = xbk(1)*xbk(2)*stot
+         if (lpp(1).eq.9)then
+            call get_dummy_x1(sjac, xbk(1), x(ndim), pi1, pi2, stot, s(-nbranch))
+            xbk(2) = 1d0
+         else
+            call get_dummy_x1(sjac, xbk(2), x(ndim), pi1, pi2, stot, s(-nbranch))
+            xbk(1) = 1d0
+         endif
+         if (.not.set_cm_rap)then
+            cm_rap=.5d0*dlog(xbk(1)*ebeam(1)/(xbk(2)*ebeam(2)))
+            set_cm_rap=.true.
+         endif
       elseif (abs(lpp(1)) .ge. 1) then
          call sample_get_x(sjac,x(ndim),ndim,mincfig,0d0,1d0)
          xbk(1) = x(ndim)
@@ -360,7 +311,10 @@
 c
 c     First Generate Momentum for initial state particles
 c
-      if(nincoming.eq.2) then
+      if (lpp(1).eq.9.or.lpp(2).eq.9)then
+         p(:,1) = pi1(:)
+         p(:,2) = pi2(:)
+      else if(nincoming.eq.2) then
         call mom2cx(m(-nbranch),m(1),m(2),1d0,0d0,p(0,1),p(0,2))
       else
         do i=0,3
@@ -548,6 +502,138 @@
       end
 
 
+      subroutine configure_integral(iconfig,mincfig,maxcfig,invar,maxwgt)
+c**************************************************************************
+c     inputs  iconfig   == Current configuration working on
+c     output  m
+c**************************************************************************
+
+      implicit none
+
+      include 'genps.inc'
+      include 'maxconfigs.inc'
+      include 'nexternal.inc'
+      include 'maxamps.inc'
+      include 'run.inc'
+
+c     local
+      double precision pi1(0:3),pi2(0:3),p0,p3
+      double precision dum
+      integer i,j,ipole,n
+      integer nbranch,ndim,nconfigs
+      integer ninvar
+      integer nparticles,nfinal
+
+
+c
+c     Arguments
+c
+      integer iconfig,mincfig,maxcfig,invar
+      double precision maxwgt
+c
+c     External
+c
+      double precision lambda,dot,dsig
+      logical passcuts
+
+
+      logical firsttime
+      data firsttime/.true./
+      save firsttime
+c
+c     global
+c
+      double precision M(-max_branch:max_particles)
+
+      double precision pmass(nexternal)
+      common/to_mass/  pmass
+
+      double precision stot,m1,m2
+      common/to_stot/stot,m1,m2
+
+      integer            mapconfig(0:lmaxconfigs), this_config
+      common/to_mconfigs/mapconfig, this_config
+
+      integer           Minvar(maxdim,lmaxconfigs)
+      common /to_invar/ Minvar
+
+      integer          lwgt(0:maxconfigs,maxinvar)
+      common/to_lwgt/lwgt
+
+      if (firsttime)then
+         firsttime=.false.
+         do i=1,nexternal
+            m(i)=pmass(i)
+         enddo
+c        Set stot
+         if (nincoming.eq.1) then
+            stot=m(1)**2
+         else
+            m1=m(1)
+            m2=m(2)
+            if (abs(lpp(1)) .eq. 1 .or. abs(lpp(1)) .eq. 2) m1 = 0.938d0
+            if (abs(lpp(2)) .eq. 1 .or. abs(lpp(2)) .eq. 2) m2 = 0.938d0
+            if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
+            if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
+            if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1
+            if(ebeam(2).lt.m2.and.lpp(2).ne.9) ebeam(2)=m2
+            pi1(0)=ebeam(1)
+            pi1(3)=sqrt(max(ebeam(1)**2-m1**2, 0d0))
+            pi2(0)=ebeam(2)
+            pi2(3)=-sqrt(max(ebeam(2)**2-m2**2, 0d0))
+            stot=m1**2+m2**2+2*(pi1(0)*pi2(0)-pi1(3)*pi2(3))
+         endif
+         write(*,'(x,a,f13.2)') 'Set CM energy to ',sqrt(stot)
+         endif
+c        Start graph mapping
+         do i=1,mapconfig(0)
+            if (mapconfig(i) .eq. iconfig) this_config=i
+         enddo
+         write(*,*) 'Mapping Graph',iconfig,' to config',this_config
+         iconfig = this_config
+         nconfigs = 1
+         mincfig=iconfig
+         maxcfig=iconfig
+         call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
+         maxwgt=0d0
+c         write(*,'(a,12i4)') 'Summing configs',(isym(i),i=1,isym(0))
+         nparticles   = nexternal
+         nfinal       = nparticles-nincoming
+         nbranch      = nparticles-2
+         ndim         = 3*nfinal-4
+         if (ndim .lt. 0) ndim = 0   !For 2->1 processes  tjs 5/24/2010
+         if (abs(lpp(1)) .ge. 1) ndim=ndim+1
+         if (abs(lpp(2)) .ge. 1) ndim=ndim+1
+         call set_peaks
+         do j=1,invar
+            lwgt(0,j)=0
+         enddo
+c
+c     Here we set up which diagrams contribute to each variable
+c     in principle more than 1 diagram can contribute to a variable
+c     if we believe they will have identical structure.
+c
+c         do i=1,mapconfig(0)
+         do i=mincfig,maxcfig
+c         do k=1,isym(0)
+c            i = isym(k)
+            write(*,'(15i4)') i,(minvar(j,i),j=1,ndim)
+            do j=1,ndim
+               ipole = minvar(j,i)
+               if (ipole .ne. 0) then
+                  n = lwgt(0,ipole)+1
+                  lwgt(n,ipole)=mapconfig(i)
+                  lwgt(0,ipole)=n
+               endif
+            enddo
+         enddo
+
+c     Initialize dsig (needed for subprocess group running mode)
+         dum=dsig(0,0,1)
+
+      return
+      end
+
       subroutine one_tree(itree,iconfig,nbranch,P,M,S,X,jac,pswgt)
 c************************************************************************
 c     Calculates the momentum for everything below in the tree until
@@ -1180,3 +1266,7 @@
       X2 = SQRT(TAU)*EXP(-ETA)
 
       END
+      
+
+
+
 
=== modified file 'Template/LO/SubProcesses/makefile'
--- Template/LO/SubProcesses/makefile	2017-03-13 20:27:08 +0000
+++ Template/LO/SubProcesses/makefile	2017-08-16 21:26:30 +0000
@@ -32,7 +32,7 @@
 
 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 \
+	 idenparts.o dummy_fct.o \
          $(patsubst %.f,%.o,$(wildcard auto_dsig*.f)) \
          $(patsubst %.f,%.o,$(wildcard matrix*.f))
 
@@ -64,6 +64,7 @@
 driver.f: genps.inc
 symmetry.o: genps.inc nexternal.inc configs.inc run_config.inc
 genps.o: genps.inc nexternal.inc configs.inc
+dummy_fct.0: run.inc genps.inc
 cuts.o: genps.inc nexternal.inc pmass.inc
 setcuts.o: genps.inc run_config.inc
 invarients.o: genps.inc nexternal.inc
 
=== modified file 'Template/LO/SubProcesses/refine.sh'
--- Template/LO/SubProcesses/refine.sh	2016-06-12 00:01:23 +0000
+++ Template/LO/SubProcesses/refine.sh	2017-08-16 21:26:30 +0000
@@ -89,7 +89,7 @@
      if [ "$keeplog" = true ] ; then
 	 echo "" >> $k; echo "ls status:" >> $k; ls >> $k	 
      else
-	 rm ftn26 > /dev/null
+	 rm ftn26 &> /dev/null
      fi
 
 
 
=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2016-09-03 13:26:00 +0000
+++ Template/LO/SubProcesses/unwgt.f	2017-08-16 21:26:30 +0000
@@ -603,7 +603,7 @@
       if (nincoming.eq.2) then
          if (xbk(1) .gt. 0d0 .and. xbk(1) .le. 1d0 .and.
      $       xbk(2) .gt. 0d0 .and. xbk(2) .le. 1d0) then
-            if(xbk(1).eq.1d0.or.pmass(1).eq.0d0) then
+            if(lpp(2).ne.0.and.(xbk(1).eq.1d0.or.pmass(1).eq.0d0)) then
                ! construct the beam momenta in each frame and compute the related (z)boost
                ebi(0) = p(0,1)/xbk(1) ! this assumes that particle 1 is massless or mass equal to beam
                ebi(1) = 0
 
=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2016-09-08 23:15:34 +0000
+++ Template/NLO/Cards/run_card.dat	2017-08-16 21:26:30 +0000
@@ -36,7 +36,7 @@
 # Normalize the weights of LHE events such that they sum or average to *
 # the total cross section                                              *
 #***********************************************************************
- %(event_norm)s = event_norm      ! average or sum
+ %(event_norm)s = event_norm    ! valid settings: average, sum, bias
 #***********************************************************************
 # Number of points per itegration channel (ignored for aMC@NLO runs)   *
 #***********************************************************************
 
=== modified file 'Template/NLO/Cards/shower_card.dat'
--- Template/NLO/Cards/shower_card.dat	2016-02-19 13:21:36 +0000
+++ Template/NLO/Cards/shower_card.dat	2017-08-16 21:26:30 +0000
@@ -27,7 +27,7 @@
 #***********************************************************************
 # PDFs and non-perturbative modelling                                  *
 #***********************************************************************
-pdfcode      = 0      # 0 = internal, 1 = same as NLO, other = lhaglue
+pdfcode      = 1     # 0 = internal, 1 = same as NLO, other = lhaglue
 ue_enabled   = F      # underlying event
 hadronize    = T      # hadronisation on/off        !IGNORED BY HERWIG6!
 lambda_5     = -1     # Lambda_5 (< 0 = default)    !IGNORED BY PYTHIA8!
 
=== modified file 'Template/NLO/FixedOrderAnalysis/HwU.f'
--- Template/NLO/FixedOrderAnalysis/HwU.f	2017-02-03 13:37:36 +0000
+++ Template/NLO/FixedOrderAnalysis/HwU.f	2017-08-16 21:26:30 +0000
@@ -1,7 +1,7 @@
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C                                                                          C
 C                  HwU: Histograms with Uncertainties                      C
-C                   By Rikkert Frederix, Dec. 2014                         C
+C                 By Rikkert Frederix, 12-2014--05-2017                    C
 C                                                                          C
 C     Book, fill and write out histograms. Particularly suited for NLO     C
 C     computations with correlations between points (ie. event and         C
@@ -9,26 +9,41 @@
 C     and PDF uncertainties through reweighting).                          C
 C                                                                          C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
+c The module contains effectively the common block with allocatable
+c variables (something not possible in old fortran version)
+      module HwU_variables
+         implicit none
+         integer :: max_plots,max_points,max_bins,nwgts,np
+         integer :: error_estimation=3
+         logical, allocatable :: booked(:)
+         integer, allocatable :: nbin(:),histi(:,:),p_bin(:),p_label(:)
+         character(len=50), allocatable :: title(:)
+         character(len=50), allocatable :: wgts_info(:)
+         double precision, allocatable :: histy(:,:,:),histy_acc(:,:,:)
+     $        ,histy2(:,:),histy_err(:,:),histxl(:,:),histxm(:,:)
+     $        ,step(:),p_wgts(:,:)
+         save
+      end module HwU_variables
+
       
 c To be called once at the start of each run. Initialises the packages
 c and sets the number of weights that need to be included for each point.
       subroutine HwU_inithist(nweights,wgt_info)
+      use HwU_variables
       implicit none
-      include "HwU.inc"
       integer i,nweights
       character*(*) wgt_info(*)
-      do i=1,max_plots
-         booked(i)=.false.
-      enddo
+      call HwU_deallocate_all
+      max_plots=0
+      max_points=0
+      max_bins=0
+      np=0
 c     Number of weights associated to each point. Note that the first
 c     weight should always be the 'central value' and it should not be
 c     zero if any of the other weights are non-zero.
       nwgts=nweights
-      if (nwgts.gt.max_wgts) then
-         write (*,*) 'ERROR: increase max_wgts in HwU histogramming'
-     $        ,max_wgts,nwgts
-         stop 1
-      endif
+      allocate(wgts_info(nwgts))
       do i=1,nwgts
          wgts_info(i)=wgt_info(i)
       enddo
@@ -48,10 +63,9 @@
 c       square root of the sum of the squares. Perform a weighted average
 c       iteration-by-iteration
 c     input 3: Same as input 2, but weighted average is same as from MINT
+      use HwU_variables
       implicit none
       integer input
-      integer error_estimation
-      common /HwU_common2/ error_estimation
       if (input.ge.0 .and. input.le.3) then
          error_estimation=input
       else
@@ -60,39 +74,20 @@
       endif
       return
       end
-      
-      block data HwU
-c set the default for the error estimation method. To reduce the size of
-c the executable, put the error_estimation variable in a separate common
-c block. If we would have included the 'HwU.inc' file here, that
-c complete common block seems to be included in the size executable
-c (approx. 110 MB).
-      integer error_estimation
-      common /HwU_common2/ error_estimation
-      data error_estimation /3/
-      end
-      
+            
 c Book the histograms at the start of the run. Give a 'label' (an
 c integer) that identifies the plot when filling it and a title
 c ('title_l') for each plot. Also the number of bins ('nbin_l') and the
 c plot range (from 'xmin' to 'xmax') should be given.
       subroutine HwU_book(label,title_l,nbin_l,xmin,xmax)
+      use HwU_variables
       implicit none
-      include "HwU.inc"
       integer label,nbin_l,i,j
       character*(*) title_l
       double precision xmin,xmax
-c     Check that label and number of bins are reasonable
-      if (label.gt.max_plots) then
-         write (*,*) 'ERROR: increase max_plots in HwU histogramming'
-     $        ,max_plots, label
-         stop 1
-      endif
-      if (nbin_l.gt.max_bins) then
-         write (*,*) 'ERROR: increase max_bins in HwU histogramming'
-     $        ,max_bins,nbin_l
-         stop 1
-      endif
+c     Allocate space for new histograms if needed      
+      call HwU_allocate_histo(label,nbin_l)
+c     Setup the histogram
       booked(label)=.true.
       title(label)=title_l
       nbin(label)=nbin_l
@@ -111,7 +106,6 @@
          histy2(label,i)=0d0
          histy_err(label,i)=0d0
       enddo
-      np=0
       return
       end
       
@@ -121,8 +115,8 @@
 c the 'HwU_inithist' subroutine. That means that each point should have
 c the same number of weights.
       subroutine HwU_fill(label,x,wgts)
+      use HwU_variables
       implicit none
-      include "HwU.inc"
       integer label,i,j,bin
       double precision x, wgts(*)
 c     If central weight is zero do not add this point.
@@ -147,11 +141,7 @@
       enddo
 c     If a new bin, add it to the list of points
       np=np+1
-      if (np.gt.max_points) then
-         write (*,*) 'ERROR: increase max_points in HwU histogramming'
-     $        ,max_points
-         stop 1
-      endif
+      call HwU_allocate_p
       p_label(np)=label
       p_bin(np)=bin
       do j=1,nwgts
@@ -159,7 +149,7 @@
       enddo
       return
       end
-
+      
 c Call after all correlated contributions for a give phase-space
 c point. I.e., every time you get a new set of random numbers from
 c MINT/VEGAS. It adds the current list of points to the histograms. Add
@@ -168,8 +158,8 @@
 c this way, correlations between events and counter-events can be
 c correctly taken into account.
       subroutine HwU_add_points
+      use HwU_variables
       implicit none
-      include "HwU.inc"
       integer i,j
       do i=1,np
          do j=1,nwgts
@@ -192,12 +182,11 @@
 c the current iteration so that they can be filled with the next
 c iteration.
       subroutine HwU_accum_iter(inclde,nPSpoints,values)
+      use HwU_variables
       implicit none
-      include "HwU.inc"
       logical inclde
       integer nPSpoints,label,i,j
-      double precision nPSinv,etot,vtot(max_wgts),niter,y_squared
-     $     ,values(2)
+      double precision nPSinv,etot,niter,y_squared,values(2)
       data niter /0d0/
       nPSinv = 1d0/dble(nPSpoints)
       if (inclde) niter = niter+1d0
@@ -226,8 +215,8 @@
 c intermediate stages this function can be called (together with
 c HwU_output) to write intermediate plots to disk.
       subroutine finalize_histograms(nPSpoints)
+      use HwU_variables
       implicit none
-      include "HwU.inc"
       integer label,nPSpoints,i,j
       double precision nPSinv,niter,dummy(2)
       nPSinv=1d0/dble(nPSpoints)
@@ -257,13 +246,13 @@
 c histograms the central value should not be zero if any of the other
 c weights are non-zero.
       subroutine accumulate_results(label,nPSinv,niter,values)
+      use HwU_variables
       implicit none
-      include "HwU.inc"
       integer label,i,j
-      double precision nPSinv,etot,vtot(max_wgts),niter,y_squared
+      double precision nPSinv,etot,niter,y_squared
      $     ,values(2),a1,a2
-      integer error_estimation
-      common /HwU_common2/ error_estimation
+      double precision,allocatable :: vtot(:)
+      if (.not. allocated(vtot)) allocate(vtot(nwgts))
       if (error_estimation.eq.2) then
 c     Use the weighted average bin-by-bin. This is not really justified
 c     for fNLO computations, because for bins with low statistics, the
@@ -395,14 +384,15 @@
 c Write the histograms to disk at the end of the run, multiplying the
 c output by 'xnorm'
       subroutine HwU_output(unit,xnorm)
+      use HwU_variables
       implicit none
-      include "HwU.inc"
       integer unit,i,j,label
       integer max_length
-      parameter (max_length=(max_wgts+3)*17)
-      character*(max_length) buffer
+      character(len=:), allocatable :: buffer
       character*4 str_nbin
       double precision xnorm
+      if (.not. allocated(buffer))
+     &     allocate(character(len=(nwgts+3)*17) :: buffer)
 c     column info: x_min, x_max, y (central value), dy, {extra
 c     weights}.
       write (unit,'(a$)') '##& xmin'
@@ -440,9 +430,156 @@
          write (unit,'(a)') ''
          write (unit,'(a)') ''         
       enddo
-      return
-      end
-
+      deallocate(buffer)
+      return
+      end
+
+c Clean all the allocatable variables:
+      subroutine HwU_deallocate_all
+      use HwU_variables
+      implicit none
+      if (allocated(wgts_info)) deallocate(wgts_info)
+      if (allocated(booked)) deallocate(booked)
+      if (allocated(title)) deallocate(title)
+      if (allocated(nbin)) deallocate(nbin)
+      if (allocated(step)) deallocate(step)
+      if (allocated(histxl)) deallocate(histxl)
+      if (allocated(histxm)) deallocate(histxm)
+      if (allocated(histy)) deallocate(histy)
+      if (allocated(histy_acc)) deallocate(histy_acc)
+      if (allocated(histi)) deallocate(histi)
+      if (allocated(histy2)) deallocate(histy2)
+      if (allocated(histy_err)) deallocate(histy_err)
+      if (allocated(p_bin)) deallocate(p_bin)
+      if (allocated(p_label)) deallocate(p_label)
+      if (allocated(p_wgts)) deallocate(p_wgts)
+      return
+      end
+
+
+      subroutine HwU_allocate_p
+      use HwU_variables
+      implicit none
+      integer,allocatable :: itemp1(:)
+      double precision, allocatable :: temp2(:,:)
+      if (.not. allocated(p_bin)) then
+         allocate(p_bin(max_plots))
+         allocate(p_label(max_plots))
+         allocate(p_wgts(nwgts,max_plots))
+         max_points=max_plots
+      else
+         if (np.gt.max_points) then
+c p_bin
+            allocate(itemp1(np+max_plots))
+            itemp1(1:max_points)=p_bin
+            call move_alloc(itemp1,p_bin)
+      
+c p_label
+            allocate(itemp1(np+max_plots))
+            itemp1(1:max_points)=p_label
+            call move_alloc(itemp1,p_label)
+c p_wgts
+            allocate(temp2(nwgts,np+max_plots))
+            temp2(1:nwgts,1:max_points)=p_wgts
+            call move_alloc(temp2,p_wgts)
+            max_points=np+max_plots
+         endif
+      endif
+      return
+      end
+
+      subroutine HwU_allocate_histo(label,nbin_l)
+      use HwU_variables
+      implicit none
+      logical,allocatable :: ltemp(:)
+      integer,allocatable :: itemp1(:),itemp2(:,:)
+      character(len=50),allocatable :: ctemp(:)
+      double precision, allocatable :: temp1(:),temp2(:,:),temp3(:,:,:)
+      integer label,i,nbin_l,label_max,nbin_max
+      logical debug
+      parameter (debug=.true.)
+c Check if variables are already allocated. If not, simply allocate a
+c single histogram
+      if (.not. allocated(booked)) then
+         allocate(booked(1))
+         booked(1)=.false.
+         allocate(title(1))
+         allocate(nbin(1))
+         allocate(step(1))
+         allocate(histxl(1,nbin_l))
+         allocate(histxm(1,nbin_l))
+         allocate(histy(nwgts,1,nbin_l))
+         allocate(histy_acc(nwgts,1,nbin_l))
+         allocate(histi(1,nbin_l))
+         allocate(histy2(1,nbin_l))
+         allocate(histy_err(1,nbin_l))
+         max_plots=1
+         max_bins=nbin_l
+      endif
+c If current label is greater than the plots already allocated, increase
+c the size of the allocated arrays. This is kind of slow, but shouldn't
+c really matter since it's only done at the start of a run.
+      if (label.gt.max_plots .or. nbin_l.gt.max_bins) then
+         label_max=max(label,max_plots)
+         nbin_max=max(nbin_l,max_bins)
+c booked
+         allocate(ltemp(label_max))
+         ltemp(1:max_plots)=booked
+         call move_alloc(ltemp,booked)
+         do i=max_plots+1,label_max
+            booked(i)=.false. ! histos have not yet been setup
+         enddo
+c title
+         allocate(ctemp(label_max))
+         ctemp(1:max_plots)=title
+         call move_alloc(ctemp,title)
+c nbin
+         allocate(itemp1(label_max))
+         itemp1(1:max_plots)=nbin
+         call move_alloc(itemp1,nbin)
+c step
+         allocate(temp1(label_max))
+         temp1(1:max_plots)=step
+         call move_alloc(temp1,step)
+c histxl
+         allocate(temp2(label_max,nbin_max))
+         temp2(1:max_plots,1:max_bins)=histxl
+         call move_alloc(temp2,histxl)
+c histxm
+         allocate(temp2(label_max,nbin_max))
+         temp2(1:max_plots,1:max_bins)=histxm
+         call move_alloc(temp2,histxm)
+c histy
+         allocate(temp3(nwgts,label_max,nbin_max))
+         temp3(1:nwgts,1:max_plots,1:max_bins)=histy
+         call move_alloc(temp3,histy)
+c histy_acc
+         allocate(temp3(nwgts,label_max,nbin_max))
+         temp3(1:nwgts,1:max_plots,1:max_bins)=histy_acc
+         call move_alloc(temp3,histy_acc)
+c histi
+         allocate(itemp2(label_max,nbin_max))
+         itemp2(1:max_plots,1:max_bins)=histi
+         call move_alloc(itemp2,histi)
+c histy2
+         allocate(temp2(label_max,nbin_max))
+         temp2(1:max_plots,1:max_bins)=histy2
+         call move_alloc(temp2,histy2)
+c histy_err
+         allocate(temp2(label_max,nbin_max))
+         temp2(1:max_plots,1:max_bins)=histy_err
+         call move_alloc(temp2,histy_err)
+c Update maximums
+         max_plots=label_max
+         max_bins=nbin_max
+      elseif (booked(label)) then
+         write (*,*) 'ERROR in HwU.f: histogram already booked',label
+         stop
+      endif
+      return
+      end
+      
+      
 c dummy subroutine
       subroutine accum(idummy)
       integer idummy
@@ -451,3 +588,6 @@
       subroutine addfil(string)
       character*(*) string
       end
+
+
+
 
=== removed file 'Template/NLO/FixedOrderAnalysis/HwU.inc'
--- Template/NLO/FixedOrderAnalysis/HwU.inc	2016-02-18 14:05:45 +0000
+++ Template/NLO/FixedOrderAnalysis/HwU.inc	1970-01-01 00:00:00 +0000
@@ -1,22 +0,0 @@
-*     -*-fortran-*-
-
-      integer max_plots,max_bins,max_wgts,max_points
-      parameter (max_plots=200)
-      parameter (max_bins=100)
-      parameter (max_wgts=1024)
-      parameter (max_points=max_plots*40)
-
-      logical booked(max_plots)
-      integer nbin(max_plots),nwgts,np,p_bin(max_points)
-     &     ,p_label(max_points),histi(max_plots,max_bins)
-      character*50 title(max_plots)
-      character*50 wgts_info(max_wgts)
-      double precision histy(max_wgts,max_plots,max_bins)
-     $     ,histy_acc(max_wgts,max_plots,max_bins),histy2(max_plots
-     $     ,max_bins),histy_err(max_plots,max_bins),histxl(max_plots
-     $     ,max_bins),histxm(max_plots,max_bins),step(max_plots)
-     $     ,p_wgts(max_wgts,max_points)
-
-      common/HwU_common/histy,histy_acc,histy2,histy_err,histxl,histxm
-     &     ,p_wgts,step,histi,nbin,p_bin,p_label,np,nwgts
-     &     ,booked,title,wgts_info
 
=== modified file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_template.f'
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_template.f	2014-12-03 11:51:14 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_template.f	2017-08-16 21:26:30 +0000
@@ -17,12 +17,12 @@
 c     o) The first argument is an integer that labels the histogram. In
 c     the analysis_end and analysis_fill subroutines this label is used
 c     to keep track of the histogram. The label should be a number
-c     between 1 and max_plots=200 (can be increased in HwU.inc).
+c     starting at 1 and be increased for each plot.
 c     o) The second argument is a string that will apear above the
 c     histogram. Do not use brackets "(" or ")" inside this string.
 c     o) The third, forth and fifth arguments are the number of bis, the
 c     lower edge of the first bin and the upper edge of the last
-c     bin. There is a maximum of 100 bins per histogram.
+c     bin.
 c     o) When including scale and/or PDF uncertainties, declare a
 c     histogram for each weight, and compute the uncertainties from the
 c     final set of histograms
 
=== removed symlink 'Template/NLO/MCatNLO/include/HwU.inc'
=== target was u'../../FixedOrderAnalysis/HwU.inc'
=== modified symlink 'Template/NLO/MCatNLO/include/reweight0.inc'
=== target was u'../../SubProcesses/reweight0.inc'
--- Template/NLO/MCatNLO/include/reweight0.inc	1970-01-01 00:00:00 +0000
+++ Template/NLO/MCatNLO/include/reweight0.inc	2017-08-16 21:26:30 +0000
@@ -0,0 +1,3 @@
+      integer max_weight_shower,maxscales,maxPDFs
+      parameter (max_weight_shower=1024,maxscales=9,maxPDFs=200)
+      integer iwgtnumpartn,jwgtinfo,mexternal
 
=== modified file 'Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f'
--- Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2014-12-09 10:04:28 +0000
+++ Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2017-08-16 21:26:30 +0000
@@ -323,7 +323,7 @@
 C---EVENTS ARE NORMALIZED TO SUM OR AVERAGE TO THE TOTAL CROSS SECTION
       WRITE(*,*)'How are the events normalized ("average" or "sum")?'
       READ(*,*)NORM_EVENT
-      if (NORM_EVENT.eq.'average')MQQ=1
+      if (NORM_EVENT(1:3).ne.'sum')MQQ=1
 
       MSFLAG=0
       IF (LHSOFT) THEN
@@ -397,3 +397,21 @@
       CALL RCLOS()
  999  STOP
       END
+
+
+c     dummy routines for stdhep
+      SUBROUTINE PYERRM(MERR,CHMESS)
+      implicit none
+      integer MERR
+      CHARACTER CHMESS*(*)
+      write(*,*)'dummy PYERRM should never be called from HW6'
+      stop
+      end
+
+
+      FUNCTION PYCOMP(KF)
+      implicit none
+      integer KF,PYCOMP
+      write(*,*)'dummy PYCOMP should never be called from HW6'
+      stop
+      end
 
=== modified file 'Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f'
--- Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f	2017-01-17 07:54:00 +0000
+++ Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f	2017-08-16 21:26:30 +0000
@@ -80,7 +80,6 @@
              read(iunit,'(a)')string ! <rwgt>
           enddo
           if(index(string,'</event>').eq.0)then
-             wgtref=XWGTUP/MQQ
              do iww=2,nwgt      ! start at 2, because 'central value' is not part of the extra weights
                 call read_rwgt_line_wgt(iunit,ww(iww))
              enddo
 
=== modified file 'Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f'
--- Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f	2014-11-27 15:36:48 +0000
+++ Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f	2017-08-16 21:26:30 +0000
@@ -318,7 +318,7 @@
 C---EVENTS ARE NORMALIZED TO SUM OR AVERAGE TO THE TOTAL CROSS SECTION
       WRITE(*,*)'How are the events normalized ("average" or "sum")?'
       READ(*,*)NORM_EVENT
-      if (NORM_EVENT.eq.'average')MQQ=1
+      if (NORM_EVENT(1:3).ne.'sum')MQQ=1
 
 C---USER'S INITIAL CALCULATIONS
       CALL PYABEG
 
=== modified file 'Template/NLO/MCatNLO/srcPythia/madfks_pylhin.f'
--- Template/NLO/MCatNLO/srcPythia/madfks_pylhin.f	2016-03-21 08:45:42 +0000
+++ Template/NLO/MCatNLO/srcPythia/madfks_pylhin.f	2017-08-16 21:26:30 +0000
@@ -80,7 +80,6 @@
         if(jwgtinfo.eq.9)then
            if (nwgt.gt.1) then
               read(iunit,'(a)')string ! <rwgt>
-              wgtref=XWGTUP/MQQ
               do iww=2,nwgt     ! start at 2, because 'central value' is not part of the extra weights
                  call read_rwgt_line_wgt(iunit,ww(iww))
               enddo
 
=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia8.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia8.cc	2016-02-23 11:25:37 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia8.cc	2017-08-16 21:26:30 +0000
@@ -58,7 +58,7 @@
   int iEventshower=pythia.mode("Main:spareMode1");
   string evt_norm=pythia.word("Main:spareWord1");
   int iEventtot_norm=iEventtot;
-  if (evt_norm == "average"){
+  if (evt_norm != "sum"){
     iEventtot_norm=1;
   }
 
@@ -112,7 +112,7 @@
     double normhepmc;
     // Add the weight of the current event to the cross section.
     normhepmc = 1. / double(iEventshower);
-    if (evt_norm == "average") {
+    if (evt_norm != "sum") {
       sigmaTotal  += evtweight*normhepmc;
     } else {
       sigmaTotal  += evtweight*normhepmc*iEventtot;
 
=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia82.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia82.cc	2016-02-23 11:25:37 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia82.cc	2017-08-16 21:26:30 +0000
@@ -57,7 +57,7 @@
   int iEventshower=pythia.mode("Main:spareMode1");
   string evt_norm=pythia.word("Main:spareWord1");
   int iEventtot_norm=iEventtot;
-  if (evt_norm == "average"){
+  if (evt_norm != "sum"){
     iEventtot_norm=1;
   }
 
@@ -118,7 +118,7 @@
     double normhepmc;
     // Add the weight of the current event to the cross section.
     normhepmc = 1. / double(iEventshower);
-    if (evt_norm == "average") {
+    if (evt_norm != "sum") {
       sigmaTotal  += evtweight*normhepmc;
     } else {
       sigmaTotal  += evtweight*normhepmc*iEventtot;
 
=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia82_hep.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia82_hep.cc	2015-03-19 11:57:06 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia82_hep.cc	2017-08-16 21:26:30 +0000
@@ -70,7 +70,7 @@
     // the number of events read by Pythia so far
     int nSelected=pythia.info.nSelected();
 
-    if (nSelected >= iEventshower) break;
+    if (nSelected > iEventshower) break;
     if (pythia.info.isLHA() && iPrintLHA < nPrintLHA) {
       pythia.LHAeventList();
       pythia.info.list();
@@ -83,7 +83,7 @@
     double evtweight = pythia.info.weight();
     double normhepmc;
     // ALWAYS NORMALISE HEPMC WEIGHTS TO SUM TO THE CROSS SECTION
-    if (evt_norm == "average") {
+    if (evt_norm != "sum") {
       normhepmc = 1. / double(iEventshower);
     } else {
       normhepmc = double(iEventtot) / double(iEventshower);
 
=== modified file 'Template/NLO/MCatNLO/srcPythia8/Pythia8_hep.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia8_hep.cc	2015-03-19 12:02:08 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia8_hep.cc	2017-08-16 21:26:30 +0000
@@ -64,7 +64,7 @@
     // the number of events read by Pythia so far
     int nSelected=pythia.info.nSelected();
 
-    if (nSelected >= iEventshower) break;
+    if (nSelected > iEventshower) break;
     if (pythia.info.isLHA() && iPrintLHA < nPrintLHA) {
       pythia.LHAeventList();
       pythia.info.list();
@@ -77,7 +77,7 @@
     double evtweight = pythia.info.weight();
     double normhepmc;
     // ALWAYS NORMALISE HEPMC WEIGHTS TO SUM TO THE CROSS SECTION
-    if (evt_norm == "average") {
+    if (evt_norm != "sum") {
       normhepmc = 1. / double(iEventshower);
     } else {
       normhepmc = double(iEventtot) / double(iEventshower);
 
=== modified file 'Template/NLO/Source/alfas_functions_lhapdf.f'
--- Template/NLO/Source/alfas_functions_lhapdf.f	2013-12-20 05:39:24 +0000
+++ Template/NLO/Source/alfas_functions_lhapdf.f	2017-08-16 21:26:30 +0000
@@ -89,7 +89,7 @@
       ALPHAS=alphasPDF(Q)
       call cpu_time(tAfter)
       
-      tPDF = tPDF + (tAfter-tBefore)
+c      tPDF = tPDF + (tAfter-tBefore)
 
       RETURN
       END
 
=== added file 'Template/NLO/Source/extra_weights.f'
--- Template/NLO/Source/extra_weights.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/Source/extra_weights.f	2017-08-16 21:26:30 +0000
@@ -0,0 +1,43 @@
+      module extra_weights
+      
+         integer,parameter :: iwgtinfo=-5,maxscales=9,maxPDFs=200
+     $     ,maxPDFsets=25,maxdynscales=10
+         integer :: max_mom_str=1,max_mext=1,max_n_ctr=1
+         logical :: doreweight,lscalevar(maxdynscales)
+     $        ,lpdfvar(maxPDFsets)
+         integer :: iwgtnumpartn,jwgtinfo,mexternal
+     $        ,lhaPDFid(0:maxPDFsets),nmemPDF(maxPDFsets)
+     $        ,dyn_scale(0:maxdynscales),n_ctr_found,n_mom_conf
+         double precision :: wgtdegrem_xi,wgtdegrem_lxi,wgtdegrem_muF
+     $        ,wgtnstmp,wgtwnstmpmuf,wgtwnstmpmur,wgtnstmp_avgvirt
+     $        ,wgtref,scalevarR(0:maxscales),scalevarF(0:maxscales)
+     $        ,wgtxsecmu(maxscales,maxscales,maxdynscales)
+     $        ,wgtxsecPDF(0:maxPDFs,maxPDFsets),wgtbpower,wgtcpower
+     $        ,veto_multiplier,H1_factor_virt,veto_compensating_factor
+     $        ,born_wgt_veto
+         double precision,allocatable :: momenta_str(:,:,:)
+         character(len= 100) :: LHAPDFsetname(maxPDFsets) 
+         character(len=1024),allocatable :: n_ctr_str(:)
+        
+c input of cpower (checked against calculated value)
+         double precision,parameter :: cpowerinput=0d0
+c switch for running muR-dependent factor runfac=1(running)/0(fixed)
+         integer,parameter :: runfac=0
+c WARNING: If you set runfac=1 to include a muR-dependent factor
+c          make sure you modified the function rwgt_muR_dep_fac in
+c	   reweight_xsec.f and compute_cpower in fks_singular.f
+c          appropiately to include all muR dependent overall factors
+c          (except for alpha_s) in the calculation. This procedure 
+c          will be incorrect, if you miss one of the muR dependent
+c          factors or if there is a not factorizing muR dependent term.
+c	   You also have to set ren_group_coeff_in and cpowerinput
+c	   to the proper values.
+
+c first order coefficient of renormalization group equation of the
+c muR-dependent factor,
+c e.g. for masses:      ren_group_coeff = gamma0 = 3/2*C_F,
+c i.e. also for Yukawa: ren_group_coeff = gamma0
+         integer,parameter :: ren_group_coeff_in=0d0
+
+         save
+      end module extra_weights
 
=== removed file 'Template/NLO/Source/getissud.f'
--- Template/NLO/Source/getissud.f	2012-12-04 16:52:45 +0000
+++ Template/NLO/Source/getissud.f	1970-01-01 00:00:00 +0000
@@ -1,201 +0,0 @@
-C...GETISSUD performs an interpolation/extrapolation in 3 dimensions by
-C...fitting quadratic splines using 4 points in each dimension
-      double precision function getissud(ibeam,kfl,x1,x2,pt2)
-      implicit none
-
-      include 'sudgrid.inc'
-
-c   Arguments
-      integer ibeam,kfl
-      double precision x1,x2,pt2
-c   Storing values for the interpolation
-      double precision smallgrid(4,4),minigrid(4) ! pt2,x1
-c   Local variables
-      integer ipt2,ix1,ix2,ilo,ihi,i,j,k,kkfl,ipoints
-      double precision pt2i,x2i,x1i,minpoint,maxpoint,x(3)
-      integer nerr
-      data nerr/0/
-
-      getissud=0
-
-      x(1)=log(x2)
-      x(2)=x1
-      x(3)=log(pt2)
-
-      kkfl=kfl
-      if(ibeam.lt.0) kkfl=-kkfl
-      if(kkfl.lt.-2) kkfl=iabs(kfl)
-      if(iabs(kkfl).eq.21) kkfl=0
-      if(kkfl.eq.5) then
-        ipoints=2
-      else
-        ipoints=1
-      endif        
-      if(kkfl.gt.5) then
-        if(nerr.lt.10)
-     $     write(*,*)'GETISSUD Warning: flavor ',kfl,' not supported'
-        nerr=nerr+1
-        getissud=1
-        return
-      endif
-
-      if(x(1).lt.points(1,ipoints).or.
-     $   x(1).gt.points(nx2,ipoints).and.nerr.lt.10)
-     $   then
-        write(*,*) 'GETISSUD Warning: extrapolation in x2: ',x2
-        nerr=nerr+1
-      endif
-
-      if(x(2).lt.points(nx2+1,ipoints).or.
-     $   x(2).gt.points(nx2+nx1,ipoints)
-     $   .and.nerr.lt.10) then
-        write(*,*) 'GETISSUD Warning: extrapolation in x1: ',x1
-        nerr=nerr+1
-      endif
-
-      if(kkfl.eq.5.and.pt2.lt.22.3109)then
-        getissud=1d0
-        return
-      endif
-
-      if(kkfl.eq.5.and.x1.gt.0.6)then
-        getissud=0d0
-        return
-      endif
-
-      if(x(3).lt.points(nx2+nx1+1,ipoints)) then
-        write(*,*) 'GETISSUD Error! pt2 = ',exp(x(3)),' < ',
-     $     exp(points(nx2+nx1+1,ipoints)),' = min(pt2). Not allowed!'
-        write(*,*) 'You need to regenerate grid with new pt2min.'
-        stop
-      endif
-
-      if(x(3).lt.points(nx2+nx1+1,ipoints).or.
-     $   x(3).gt.points(nx2+nx1+npt2,ipoints)
-     $   .and.nerr.lt.10) then
-        write(*,*) 'GETISSUD Warning: extrapolation in pt2: ',pt2
-        nerr=nerr+1
-      endif
-
-
-c   Find nearest points by binary method
-c   x2
-      ilo=1
-      ihi=nx2
-      do while(ihi.gt.ilo+1)
-        ix2=ilo+(ihi-ilo)/2
-        if(x(1).gt.points(ix2,ipoints))then
-          ilo=ix2
-        else
-          ihi=ix2
-        endif
-      enddo
-      if(x(1).lt.points(ix2,ipoints))
-     $   ix2=ix2-1
-      ix2=max(2,min(ix2,nx2-2))
-
-c      print *,'x2: ',ix2,x(1),(points(i,ipoints),i=ix2-1,ix2+2)
-
-c   x1
-      ilo=1
-      ihi=nx1
-      do while(ihi.gt.ilo+1)
-        ix1=ilo+(ihi-ilo)/2
-        if(x(2).gt.points(nx2+ix1,ipoints))then
-          ilo=ix1
-        else
-          ihi=ix1
-        endif
-      enddo
-      if(x(2).lt.points(nx2+ix1,ipoints))
-     $   ix1=ix1-1
-      ix1=max(2,min(ix1,nx1-2))
-
-      do while(kkfl.eq.5.and.
-     $   points(nx2+ix1+2,ipoints).gt.0.6)
-        ix1=ix1-1
-      enddo
-
-c      print *,'x1: ',ix1,x(2),(points(nx2+i,ipoints),i=ix1-1,ix1+2)
-
-c   pt2
-      ilo=1
-      ihi=npt2
-      do while(ihi.gt.ilo+1)
-        ipt2=ilo+(ihi-ilo)/2
-        if(x(3).gt.points(nx2+nx1+ipt2,ipoints))then
-          ilo=ipt2
-        else
-          ihi=ipt2
-        endif
-      enddo
-      if(x(3).lt.points(nx2+nx1+ipt2,ipoints))
-     $   ipt2=ipt2-1
-      ipt2=max(2,min(ipt2,npt2-2))
-
-      do while(kkfl.eq.5.and.
-     $   exp(points(nx2+nx1+ipt2-1,ipoints)).lt.22.3109)
-        ipt2=ipt2+1
-      enddo
-c      print *,'pt2: ',ipt2,x(3),(points(nx2+nx1+i,ipoints),i=ipt2-1,ipt2+2)
-c      print *,'pt: ',ipt2,exp(x(3)/2),
-c     $   (exp(points(nx2+nx1+i,ipoints)/2),i=ipt2-1,ipt2+2)
-
-C   Now perform inter-/extra-polation
-
-C   Start with x2, which should have the flattest distribution
-C   Calculate sud(x2,ax1,apt2) for the 4x4 apt2 and ax1
-C   Then continue with pt2 and calculate sud(x2,ax1,pt2)
-C   for the 4 ax1
-C   Finally calculate sud(x2,x1,pt2)
-
-      do i=1,4
-        do j=1,4
-c          print *,'x1,pt:',points(nx2+ix1-2+i,ipoints),
-c     $       exp(points(nx2+nx1+ipt2-2+j,ipoints)/2)
-          call splint2(sudgrid(ix2-1,ix1-2+i,ipt2-2+j,kkfl),
-     $       points(ix2-1,ipoints),4,x(1),smallgrid(j,i))
-          smallgrid(j,i)=max(0d0,min(1d0,smallgrid(j,i)))
-        enddo
-      enddo
-
-      do i=1,4
-        call splint2(smallgrid(1,i),
-     $     points(nx2+nx1+ipt2-1,ipoints),4,x(3),minigrid(i))
-          minigrid(i)=max(0d0,min(1d0,minigrid(i)))
-      enddo
-
-      call splint2(minigrid,
-     $   points(nx2+ix1-1,ipoints),4,x(2),getissud)
-      getissud=max(0d0,min(1d0,getissud))
-      
-c      print *,'Result: ',getissud
-
-      return
-      end
-
-
-      subroutine splint2(ypoints,xpoints,npoints,x,ans)
-      implicit none
-
-C   arguments
-      integer npoints
-      double precision ypoints(npoints),xpoints(npoints)
-      double precision x,ans
-C   local variables
-      double precision a0,a1,a2,sd
-      integer ifail,i
-
-      CALL DLSQP2(npoints,xpoints,ypoints,a0,a1,a2,sd,ifail)
-
-c      print *,'Point, interpolation:'
-c      do i=1,npoints
-c        print *,exp(xpoints(i)),ypoints(i),
-c     $     a0+a1*xpoints(i)+a2*xpoints(i)**2
-c      enddo
-
-      ans=a0+a1*x+a2*x**2
-c      print *,x,ans
-
-      return
-      end
 
=== modified file 'Template/NLO/Source/make_opts.inc'
--- Template/NLO/Source/make_opts.inc	2016-06-28 14:44:52 +0000
+++ Template/NLO/Source/make_opts.inc	2017-08-16 21:26:30 +0000
@@ -48,11 +48,6 @@
   F2PY=$(DEFAULT_F2PY_COMPILER)
 endif
 
-# Set CXX unless it's defined by an environment variable
-ifeq ($(origin CXX),default)
-  CXX=$(DEFAULT_CPP_COMPILER)
-endif
-
 UNAME := $(shell uname -s)
 ifeq ($(origin LDFLAGS), undefined)
      LDFLAGS=$(STDLIB) $(MACFLAG)
@@ -94,10 +89,22 @@
   alfas_functions=alfas_functions_lhapdf
   llhapdf = $(shell $(lhapdf) --libs)
   reweight_xsec_events_pdf_dummy=
+  # check if we need to activate c++11 (for lhapdf6.2)
+  ifeq ($(origin CXX),default)
+    ifeq ($lhapdfversion$lhapdfsubversion,62)
+      CXX=$(DEFAULT_CPP_COMPILER) -std=c++11
+    else
+      CXX=$(DEFAULT_CPP_COMPILER)
+    endif
+  endif	     
 else
   alfas_functions=alfas_functions
   llhapdf=
   reweight_xsec_events_pdf_dummy=reweight_xsec_events_pdf_dummy.o
+  # Set CXX unless it's defined by an environment variable
+  ifeq ($(origin CXX),default)
+    CXX=$(DEFAULT_CPP_COMPILER)
+  endif
 endif
 
 # Option APPLGrid
 
=== modified file 'Template/NLO/Source/makefile'
--- Template/NLO/Source/makefile	2017-05-09 13:28:25 +0000
+++ Template/NLO/Source/makefile	2017-08-16 21:26:30 +0000
@@ -8,12 +8,8 @@
 IREGIDIR= ./IREGI/src/
 STDHEPDIR= ./StdHEP/
 
-PROCESS= matrix.o myamp.o 
-
-GENERIC = $(alfas_functions).o rw_routines.o kin_functions.o setrun.o	\
-	  run_printout.o dgauss.o ranmar.o getissud.o
-
-INCLUDEF = coupl.inc pmass.inc cluster.inc
+GENERIC = $(alfas_functions).o rw_routines.o kin_functions.o	\
+	  run_printout.o dgauss.o ranmar.o setrun.o
 
 .f.o: ; $(FC) $(FFLAGS) -c $*.f
 
@@ -21,19 +17,19 @@
     $(LIBDIR)libmodel.a $(LIBDIR)libcernlib.a param_card.inc
 	rm -f PDF/*.o
 
-
 $(LIBDIR)libdhelas.a: DHELAS
 	rm -f $(LIBDIR)libdhelas.a
 	cd DHELAS; make
 
 $(LIBDIR)libgeneric.a: $(GENERIC)
 	rm -f $(LIBDIR)libgeneric.a
-	ar cru libgeneric.a $(GENERIC)
+	ar cru libgeneric.a $(GENERIC) extra_weights.o
 	ranlib libgeneric.a
 	mv libgeneric.a $(LIBDIR)
+	cp -f extra_weights.mod $(LIBDIR)
 	rm -f $(alfas_functions).o 
 
-setrun.o: run_card.inc 
+setrun.o: run_card.inc extra_weights.o
 	$(FC) $(FFLAGS) -c -o setrun.o setrun.f
 
 $(LIBDIR)libpdf.a: PDF
 
=== modified file 'Template/NLO/Source/run.inc'
--- Template/NLO/Source/run.inc	2017-01-13 15:01:41 +0000
+++ Template/NLO/Source/run.inc	2017-08-16 21:26:30 +0000
@@ -54,11 +54,6 @@
       double precision bwcutoff
       common/to_bwcutoff/ bwcutoff
 c
-c     Sudakov grid file name
-c
-      character*130 issgridfile
-      common/to_sgridfile/issgridfile
-c
 c     kT/pT scheme for xqcut, clustering according to channel
 c
       integer ktscheme
 
=== removed file 'Template/NLO/Source/run_config.inc'
--- Template/NLO/Source/run_config.inc	2015-09-24 16:36:04 +0000
+++ Template/NLO/Source/run_config.inc	1970-01-01 00:00:00 +0000
@@ -1,54 +0,0 @@
-c*********************************************************************
-c     Parameters to configure running information for MadEvent
-c     The default values of these parameters should not need to be
-c     changed, unless there is a special need for optimization
-c*********************************************************************
-c
-c     The following parameters are used by symmetry.f in setting up the survey
-c
-      integer    icomp
-      parameter (icomp = 3)   !BW + Symmetry compression  0 == none	
-      integer    min_events_subprocess          !Minimum number of events
-      parameter (min_events_subprocess = 2000) !per iteration in each subprocess
-      integer    min_events_channel         !Minimum number of events
-      parameter (min_events_channel = 1000) !per iteration in each channel
-      integer    iter_survey                !Number of iterations for survey
-      parameter (iter_survey=4)
-      integer    nhel_survey                !Number of helicities to sum for each
-      parameter (nhel_survey=0)             !phase space point. (0 = sum_all)
-c
-c     The following parameters are used by gen_ximprove.f in running refine
-c
-      integer    min_events                 !Minimum number of events/iteration
-      parameter (min_events = 2000)         !to refine a channel
-      integer    max_events                 !Maximum number of events/iteration
-      parameter (max_events = 2000)        !to refine a channel
-      integer    max_iter                   !Maximum number of iterations
-      parameter (max_iter = 9)              !during refinement
-      integer    nhel_refine                !Number of helicities to sum for each
-      parameter (nhel_refine=0)             !phase space point. (0 = sum_all)
-c
-c     The following are used for parallel running
-c
-      character*(20) PBS_QUE
-      parameter     (PBS_QUE = 'madgraph')
-
-      integer    ChanPerJob
-      parameter (ChanPerJob=100000000)  !Number of channels / job for survey
-
-c      integer    max_np
-c      parameter (max_np=1)  !Number of channels / job for refine
-c
-c
-c
-	double precision trunc_max
-        parameter (trunc_max=0.01)
-c
-c	The following are used for grid type running
-c
-      double precision acc_wu
-      parameter       (acc_wu = 0.01)   !Desired accuracy from warmup run
-      integer      npoints_wu,        itmax_wu !warmup # points/iterations
-      parameter   (npoints_wu = 4000, itmax_wu = 8)
-      integer      min_gevents_wu
-      parameter   (min_gevents_wu=1000)   !Minumum # unweighted events to generate from channel
 
=== modified file 'Template/NLO/Source/setrun.f'
--- Template/NLO/Source/setrun.f	2017-05-09 11:54:11 +0000
+++ Template/NLO/Source/setrun.f	2017-08-16 21:26:30 +0000
@@ -6,47 +6,16 @@
 c 2. Collider parameters
 c 3. cuts
 c---------------------------------------------------------------------- 
+      use extra_weights
       implicit none
-c
-c     parameters
-c
-      integer maxpara
-      parameter (maxpara=1000)
-c
-c     local
-c     
-      integer npara
-      character*20 param(maxpara),value(maxpara)
-c
-c     include
-c
       include 'PDF/pdf.inc'
       include 'run.inc'
       include 'alfas.inc'
       include 'MODEL/coupl.inc'
-      include '../SubProcesses/reweight0.inc'
-
-      double precision D
-      common/to_dj/D
-c
-c     local
-c
-      character*20 ctemp
-      integer k,i,l1,l2
+      integer k,i
       character*132 buff
-C
-C     input cuts
-C
       include 'cuts.inc'
-C
-C     BEAM POLARIZATION
-C
-      REAL*8 POL(2)
-      common/to_polarization/ POL
-      data POL/1d0,1d0/
-c
-c     Les Houches init block (for the <init> info)
-c
+c Les Houches init block (for the <init> info)
       integer maxpup
       parameter(maxpup=100)
       integer idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
@@ -55,7 +24,6 @@
      &     idwtup,nprup,xsecup(maxpup),xerrup(maxpup),
      &     xmaxup(maxpup),lprup(maxpup)
 c
-      logical gridrun,gridpack
       integer          iseed
       common /to_seed/ iseed
       integer nevents
@@ -64,17 +32,13 @@
       integer iappl
       common /for_applgrid/ iappl
       integer idum
-C      
-c
-c----------
-c     start
-c----------
+c jet-rate distance. To be set to 1 for FxFx
+      double precision D
+      common/to_dj/D
+c Include all the parameters set in the run_card.dat
       include 'run_card.inc'
-      
-c MZ add the possibility to have shower_MC input lowercase
+c Change shower_MC string to upper case
       call to_upper(shower_MC)
-C
-
 c Determine if there is a need to do scale and/or PDF reweighting
       do_rwgt_scale=.false.
       do i=1,dyn_scale(0)
@@ -90,69 +54,41 @@
             exit
          endif
       enddo
-
 c Default scale and PDF choice used for the actual run
       dynamical_scale_choice=dyn_scale(1)
       lhaid=lhaPDFid(1)
-
-c merging cuts
-      xqcut=0d0
-      xmtc=0d0
-      d=1d0
-      
-c*********************************************************************
-c     Random Number Seed                                             *
-c*********************************************************************
-
-        gridrun=.false.
-        gridpack=.false.
-
-c************************************************************************     
-c     Renormalization and factorization scales                          *
-c************************************************************************     
-c
-
 c For backward compatibility
       scale = muR_ref_fixed
       q2fact(1) = muF1_ref_fixed**2      ! fact scale**2 for pdf1
       q2fact(2) = muF2_ref_fixed**2      ! fact scale**2 for pdf2     
       scalefact=muR_over_ref
       ellissextonfact=QES_over_ref
-
 c check that the event normalization input is reasoble
       buff = event_norm 
       call case_trap2(buff) ! requires a string of length 20 at least
       event_norm=buff 
-      if (event_norm(1:7).ne.'average' .and. event_norm(1:3).ne.'sum'
-     $     .and. event_norm(1:5).ne.'unity')then
+      if ( event_norm(1:7).ne.'average' .and.
+     $     event_norm(1:3).ne.'sum' .and.
+     $     event_norm(1:5).ne.'unity'.and.
+     $     event_norm(1:4).ne.'bias')then
          write (*,*) 'Do not understand the event_norm parameter'/
      &        /' in the run_card.dat. Possible options are'/
-     &        /' "average", "sum" or "unity". Current input is: ',
-     &        event_norm
-         open(unit=26,file='../../error',status='unknown')
-         write (26,*) 'Do not understand the event_norm parameter'/
-     &        /' in the run_card.dat. Possible options are'/
-     &        /' "average", "sum" or "unity". Current input is: ',
-     &        event_norm
-         
+     &        /' "average", "sum", "unity" or "bias". '/
+     &        /'Current input is: ', event_norm
          stop 1
       endif
-
-c info for reweight
-
+c Check parameters for FxFx/UNLOPS/NNLL-veto
       if ( ickkw.ne.0 .and. ickkw.ne.4 .and. ickkw.ne.3 .and.
      &     ickkw.ne.-1) then
          write (*,*) 'ickkw parameter not known. ickkw=',ickkw
          stop
       endif
-c$$$      ickkw=0
       chcluster=.false.
       ktscheme=1
-
-c !!! Default behavior changed (MH, Aug. 07) !!!
-c If no pdf, read the param_card and use the value from there and
-c order of alfas running = 2
-
+      xqcut=0d0
+      xmtc=0d0
+      D=1d0
+c Set alphaS(mZ)      
       if(lpp(1).ne.0.or.lpp(2).ne.0) then
          write(*,*) 'A PDF is used, so alpha_s(MZ)'/
      &        /' is going to be modified'
@@ -172,33 +108,29 @@
           write(*,*) 'The default order of alpha_s running is fixed to '
      &         ,nloop
       endif
-c !!! end of modification !!!
-
-C       Fill common block for Les Houches init info
+C Fill common block for Les Houches init info
       do i=1,2
-        if(lpp(i).eq.1.or.lpp(i).eq.2) then
-          idbmup(i)=2212
-        elseif(lpp(i).eq.-1.or.lpp(i).eq.-2) then
-          idbmup(i)=-2212
-        elseif(lpp(i).eq.3) then
-          idbmup(i)=11
-        elseif(lpp(i).eq.-3) then
-          idbmup(i)=-11
-        elseif(lpp(i).eq.0) then
-           open (unit=71,status='old',file='initial_states_map.dat')
-           read (71,*,err=100)idum,idum,idbmup(1),idbmup(2)
-           close (71)
-c          idbmup(i)=idup(i,1)
-        else
-          idbmup(i)=lpp(i)
-        endif
+         if(lpp(i).eq.1.or.lpp(i).eq.2) then
+            idbmup(i)=2212
+         elseif(lpp(i).eq.-1.or.lpp(i).eq.-2) then
+            idbmup(i)=-2212
+         elseif(lpp(i).eq.3) then
+            idbmup(i)=11
+         elseif(lpp(i).eq.-3) then
+            idbmup(i)=-11
+         elseif(lpp(i).eq.0) then
+            open (unit=71,status='old',file='initial_states_map.dat')
+            read (71,*,err=100)idum,idum,idbmup(1),idbmup(2)
+            close (71)
+         else
+            idbmup(i)=lpp(i)
+         endif
         ebmup(i)=ebeam(i)
       enddo
       call get_pdfup(pdlabel,pdfgup,pdfsup,lhaid)
-
+c Fill the nmemPDF(i) array with the number of PDF error set. This we
+c get from LHAPDF.
       if (lpdfvar(1) .and. (lpp(1).ne.0.or.lpp(2).ne.0) ) then
-c fill the nmemPDF(i) array with the number of PDF error set. This we
-c get from LHAPDF.
          call numberPDFm(1,nmemPDF(1))
          if (nmemPDF(1).eq.1) then
             nmemPDF(1)=0
@@ -207,100 +139,51 @@
       else
          nmemPDF(1)=0
       endif
-
-      return
- 99   write(*,*) 'error in reading'
       return
  100  write(*,*) '"initial_states_map.dat" not found (or incorrect'/
      $     /' format) by "Source/setrun"'
       stop 1
       end
 
-C-------------------------------------------------
-C   GET_PDFUP
-C   Convert MadEvent pdf name to LHAPDF number
-C-------------------------------------------------
-
       subroutine get_pdfup(pdfin,pdfgup,pdfsup,lhaid)
+C     Convert internal pdf name to LHAPDF number
       implicit none
-
       character*(*) pdfin
       integer mpdf
       integer npdfs,i,pdfgup(2),pdfsup(2),lhaid
-
       parameter (npdfs=16)
       character*7 pdflabs(npdfs)
-      data pdflabs/
-     $   'none',
-     $   'mrs02nl',
-     $   'mrs02nn',
-     $   'cteq4_m',
-     $   'cteq4_l',
-     $   'cteq4_d',
-     $   'cteq5_m',
-     $   'cteq5_d',
-     $   'cteq5_l',
-     $   'cteq5m1',
-     $   'cteq6_m',
-     $   'cteq6_l',
-     $   'cteq6l1',     
-     $   'nn23lo',
-     $   'nn23lo1',
-     $   'nn23nlo'/
+      data pdflabs/ 'none', 'mrs02nl', 'mrs02nn', 'cteq4_m', 'cteq4_l',
+     $     'cteq4_d', 'cteq5_m', 'cteq5_d', 'cteq5_l', 'cteq5m1',
+     $     'cteq6_m', 'cteq6_l', 'cteq6l1', 'nn23lo', 'nn23lo1',
+     $     'nn23nlo'/
       integer numspdf(npdfs)
-      data numspdf/
-     $   00000,
-     $   20250,
-     $   20270,
-     $   19150,
-     $   19170,
-     $   19160,
-     $   19050,
-     $   19060,
-     $   19070,
-     $   19051,
-     $   10000,
-     $   10041,
-     $   10042,
-     $   246800,
-     $   247000,
-     $   244800/
-
-
+      data numspdf/ 00000, 20250, 20270, 19150, 19170, 19160, 19050,
+     $     19060, 19070, 19051, 10000, 10041, 10042, 246800, 247000,
+     $     244800/
       if(pdfin.eq."lhapdf") then
-        write(*,*)'using LHAPDF'
-        do i=1,2
-           pdfgup(i)=-1
-           pdfsup(i)=lhaid
-        enddo
-        return
+         write(*,*)'using LHAPDF'
+         do i=1,2
+            pdfgup(i)=-1
+            pdfsup(i)=lhaid
+         enddo
+         return
       endif
-
-      
       mpdf=-1
       do i=1,npdfs
-        if(pdfin(1:len_trim(pdfin)) .eq. pdflabs(i))then
-          mpdf=numspdf(i)
-        endif
+         if(pdfin(1:len_trim(pdfin)) .eq. pdflabs(i))then
+            mpdf=numspdf(i)
+         endif
       enddo
-
       if(mpdf.eq.-1) then
-        write(*,*)'ERROR: pdf ',pdfin,' not implemented in get_pdfup.'
-        write(*,*)'known pdfs are'
-        write(*,*) pdflabs
-        open(unit=26,file='../../error',status='unknown')
-        write(26,*)'ERROR: pdf ',pdfin,' not implemented in get_pdfup.'
-        write(26,*)'known pdfs are'
-        write(26,*) pdflabs
-        stop 1
-c$$$        write(*,*)'using ',pdflabs(12)
-c$$$        mpdf=numspdf(12)
+         write(*,*)'ERROR: pdf ',pdfin,' not implemented in get_pdfup.'
+         write(*,*)'known pdfs are'
+         write(*,*) pdflabs
+         stop 1
       endif
-
       do i=1,2
-        pdfgup(i)=-1
-        pdfsup(i)=mpdf
+         pdfgup(i)=-1
+         pdfsup(i)=mpdf
       enddo
-
       return
       end
 
=== removed file 'Template/NLO/Source/sudgrid.inc'
--- Template/NLO/Source/sudgrid.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/sudgrid.inc	1970-01-01 00:00:00 +0000
@@ -1,4 +0,0 @@
-      integer npt2,nx1,nx2
-      parameter(npt2=40,nx1=80,nx2=20)
-      double precision points(nx2+nx1+npt2,2),sudgrid(nx2,nx1,npt2,-2:5)
-      common/sudgrid/points,sudgrid
 
=== modified file 'Template/NLO/SubProcesses/add_write_info.f'
--- Template/NLO/SubProcesses/add_write_info.f	2017-05-16 08:15:12 +0000
+++ Template/NLO/SubProcesses/add_write_info.f	2017-08-16 21:26:30 +0000
@@ -7,7 +7,6 @@
       include "nexternal.inc"
       include "born_nhel.inc"
       include "coloramps.inc"
-      include "reweight0.inc"
       include "nFKSconfigs.inc"
       include "leshouche_decl.inc"
       include "run.inc"
 
=== modified file 'Template/NLO/SubProcesses/analysis_lhe.f'
--- Template/NLO/SubProcesses/analysis_lhe.f	2017-02-14 15:57:59 +0000
+++ Template/NLO/SubProcesses/analysis_lhe.f	2017-08-16 21:26:30 +0000
@@ -30,29 +30,30 @@
       integer                                   npoints
       double precision            cross_section
       common /for_FixedOrder_lhe/ cross_section,npoints
+      character*10 MonteCarlo
       inquire(41,OPENED=lopen)  ! safety
       if (lopen) then
          backspace(41)          ! overwrite the final </eventgroup> tag
          write (41,*) nevents,sum_of_wgts,cross_section
          close(41)
          open(41, file='header.txt')
-         call write_lhef_header(41, 0, 'FO')
+         MonteCarlo='FO'
+         call write_lhef_header(41, 0, MonteCarlo)
          close(41)
       endif
       end
 
       subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+      use extra_weights
       implicit none
       include 'nexternal.inc'
-      include 'reweight.inc'
       include 'run.inc'
-      integer nwgt,max_weight
-      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
+      integer nwgt
 c      include 'genps.inc'
       integer istatus(nexternal)
       integer iPDG(nexternal)
       double precision p(0:4,nexternal)
-      double precision wgts(max_weight)
+      double precision wgts(*)
       integer ibody
       integer nwgts,nevents
       double precision sum_of_wgts
 
=== removed file 'Template/NLO/SubProcesses/c_weight.inc'
--- Template/NLO/SubProcesses/c_weight.inc	2016-03-21 08:45:42 +0000
+++ Template/NLO/SubProcesses/c_weight.inc	1970-01-01 00:00:00 +0000
@@ -1,24 +0,0 @@
-*     -*-fortran-*-
-
-      integer max_contr,max_wgt,max_iproc
-      parameter (max_contr=128,max_wgt=1024,max_iproc=100)
-      logical H_event(max_contr)
-      integer nFKS(max_contr),itype(max_contr),QCDpower(max_contr)
-     $     ,pdg(nexternal,0:max_contr),pdg_uborn(nexternal,0:max_contr)
-     $     ,parton_pdg_uborn(nexternal,max_iproc,max_contr)
-     $     ,parton_pdg(nexternal,max_iproc,max_contr),icontr,iwgt
-     $     ,plot_id(max_contr),niproc(max_contr),parton_pdf(nexternal
-     $     ,max_iproc,max_contr),icontr_sum(0:max_iproc,max_contr)
-     $     ,icontr_picked,iproc_picked
-      double precision momenta(0:3,nexternal,max_contr),momenta_m(0:3
-     $     ,nexternal,2,max_contr),wgt(3,max_contr),wgt_ME_tree(2
-     $     ,max_contr),bjx(2,max_contr),scales2(3,max_contr)
-     $     ,g_strong(max_contr),wgts(max_wgt,max_contr)
-     $     ,parton_iproc(max_iproc,max_contr),y_bst(max_contr)
-     $     ,plot_wgts(max_wgt,max_contr),shower_scale(max_contr)
-     $     ,unwgt(max_iproc,max_contr)
-      common /c_weights/momenta,momenta_m,wgt,wgt_ME_tree,wgts,plot_wgts
-     $     ,parton_iproc,bjx,scales2,g_strong,y_bst,shower_scale,unwgt
-     $     ,parton_pdg,nFKS,itype,QCDpower,pdg,pdg_uborn
-     $     ,parton_pdg_uborn,plot_id,niproc,icontr,iwgt,icontr_sum
-     $     ,icontr_picked,iproc_picked,H_event
 
=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2017-04-13 09:51:06 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2017-08-16 21:26:30 +0000
@@ -183,6 +183,7 @@
             enddo
           enddo
 
+          CALL UPDATE_AS_PARAM()
           call sborn(p_born, born)
           call sloopmatrix_thres(p_born,virt_wgts,tolerance,
      1 accuracies,return_code)
 
=== modified file 'Template/NLO/SubProcesses/collect_events.f'
--- Template/NLO/SubProcesses/collect_events.f	2016-02-24 12:45:15 +0000
+++ Template/NLO/SubProcesses/collect_events.f	2017-08-16 21:26:30 +0000
@@ -1,6 +1,6 @@
       program collect_events
       implicit none
-      character*120 string120,eventfile
+      character*512 string512,eventfile
       character*19 basicfile,nextbasicfile
       character*15 outputfile
       integer istep,i,numoffiles,nbunches,nevents,ievents,junit(80)
@@ -65,9 +65,9 @@
       nevents=0
       xtotal=0.d0
       do while (.true.)
-         read(10,'(120a)',err=2,end=2) string120
-         eventfile=string120(2:index(string120,'   '))
-         read(string120(index(string120,'   '):120),*)
+         read(10,'(512a)',err=2,end=2) string512
+         eventfile=string512(2:index(string512,'   '))
+         read(string512(index(string512,'   '):512),*)
      $    ievents,absxsec,xsecfrac
          if (ievents.eq.0) cycle
          nevents=nevents+ievents
@@ -209,8 +209,9 @@
       end
 
 
-      subroutine collect_all_evfiles(ioutput,numoffiles,junit,
-     #                               imaxevt,evwgt)
+      subroutine collect_all_evfiles(ioutput,numoffiles,junit,imaxevt
+     $     ,evwgt)
+      use extra_weights
       implicit none
       integer i_orig
       common /c_i_orig/i_orig
@@ -237,8 +238,6 @@
       data iseed/1/
       double precision rnd,fk88random
       external fk88random
-      logical debug
-      parameter (debug=.false.)
       integer nevents_file(80),proc_id(80)
       common /to_nevents_file/nevents_file,proc_id
       double precision xsecfrac_all(80)
@@ -249,18 +248,12 @@
       double precision xsecup_l(100),xerrup_l(100)
       integer lprup_l(100),nproc_l
       logical found_proc
-      include 'reweight_all.inc'
       include 'run.inc'
       integer proc_id_tot(0:100)
       double precision xsec(100),xsecABS,xerr(100)
       logical get_xsec_from_res1
       common/total_xsec/xsec,xerr,xsecABS,proc_id_tot,get_xsec_from_res1
 c
-      if(debug) then
-         write (*,*) ioutput,numoffiles,(junit(ii),ii=1,numoffiles)
-         write(ioutput,*)'test test test'
-         return
-      endif
       maxevt=0
       if (.not. get_xsec_from_res1) then
          do i=1,100
 
=== modified file 'Template/NLO/SubProcesses/combine_root.C'
--- Template/NLO/SubProcesses/combine_root.C	2013-12-02 16:53:51 +0000
+++ Template/NLO/SubProcesses/combine_root.C	2017-08-16 21:26:30 +0000
@@ -44,11 +44,13 @@
   //can be easily lifted
   int maxfiles=1000;
 
+  int numoffiles;
+  string filenames[maxfiles];
+  string currentdir;
+
   if(fileoffiles.is_open())
     {
-      string filenames[maxfiles];
 
-      int numoffiles;
       fileoffiles >> numoffiles;
 
       if(numoffiles>maxfiles){
@@ -56,7 +58,6 @@
 	return;
       }
 
-      string currentdir;
       fileoffiles >> currentdir;
 
       for(int i = 0; i < numoffiles; ++i)
@@ -122,8 +123,8 @@
       int llength = HistoContents.size();
       if(myverbose)cout << "Total number of histograms: " << llength << endl;
 
+      int llength0 = llength;
       if(i==0){
-	int llength0 = llength;
 	for(int j = 0; j <  llength; j++){
 	  TString histname = HistoContents[j]->GetName();
 	  if(myverbose){
 
=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2017-05-09 12:37:07 +0000
+++ Template/NLO/SubProcesses/cuts.f	2017-08-16 21:26:30 +0000
@@ -898,32 +898,47 @@
       end
 
 
-      subroutine unweight_function(p_born,unwgtfun)
-c This is a user-defined function to which to unweight the events
-c A non-flat distribution will generate events with a certain
-c weight. This is particularly useful to generate more events
-c (with smaller weight) in tails of distributions.
-c It computes the unwgt factor from the momenta and multiplies
-c the weight that goes into MINT (or vegas) with this factor.
-c Before writing out the events (or making the plots), this factor
-c is again divided out.
-c This function should be called with the Born momenta to be sure
-c that it stays the same for the events, counter-events, etc.
-c A value different from 1 makes that MINT (or vegas) does not list
-c the correct cross section.
+      subroutine bias_weight_function(p,ipdg,bias_wgt)
+c This is a user-defined function to which to bias the event generation.
+c A non-flat distribution will generate events with a certain weight
+c inversely proportinal to the bias_wgt. This is particularly useful to
+c generate more events (with smaller weight) in tails of distributions.
+c It computes the bias_wgt factor from the momenta and multiplies the
+c weight that goes into MINT (or vegas) with this factor.  Before
+c writing out the events (or making the plots), this factor is again
+c divided out. A value different from 1 makes that MINT (or vegas) does
+c not list the correct cross section, but the cross section can still be
+c computed from summing all the weights of the events (and dividing by
+c the number of events). Since the weights of the events are no longer
+c identical for all events, the statistical uncertainty on this total
+c cross section can be much larger than without including the bias.
+c
+c The 'bias_wgt' should be a IR-safe function of the momenta.
+c      
+c For this to be used, the 'event_norm' option in the run_card should be
+c set to
+c      'bias' = event_norm      
+c
       implicit none
       include 'nexternal.inc'
-      double precision unwgtfun,p_born(0:3,nexternal-1),shat,sumdot
-      external sumdot
-
-      unwgtfun=1d0
-
-c How to enhance the tails is very process dependent. But, it is
-c probably easiest to enhance the tails using shat, e.g.:
-c      shat=sumdot(p_born(0,1),p_born(0,2),1d0)
-c      unwgtfun=max(100d0**2,shat)/100d0**2
-c      unwgtfun=unwgtfun**2
-
+      double precision bias_wgt,p(0:3,nexternal),H_T
+      integer ipdg(nexternal),i
+
+      bias_wgt=1d0
+
+c How to enhance the tails is very process dependent. For example for
+c top quark production one could use:
+c      do i=1,nexternal
+c         if (ipdg(i).eq.6) then
+c            bias_wgt=sqrt(p(1,i)**2+p(2,i)**2)**3
+c         endif
+c      enddo
+c Or to use H_T^2 one does     
+c      H_T=0d0
+c      do i=3,nexternal
+c         H_T=H_T+sqrt(max(0d0,(p(0,i)+p(3,i))*(p(0,i)-p(3,i))))
+c      enddo
+c      bias_wgt=H_T**2
       return
       end
 
 
=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2017-05-09 11:54:11 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2017-08-16 21:26:30 +0000
@@ -2,6 +2,7 @@
 c**************************************************************************
 c     This is the driver for the whole calculation
 c**************************************************************************
+      use extra_weights
       implicit none
 C
 C     CONSTANTS
@@ -10,14 +11,13 @@
       parameter       (ZERO = 0d0)
       include 'nexternal.inc'
       include 'genps.inc'
-      include 'reweight.inc'
       INTEGER    ITMAX,   NCALL
 
       common/citmax/itmax,ncall
 C
 C     LOCAL
 C
-      integer i,j,l,l1,l2,ndim
+      integer i,j,l,l1,l2
       character*130 buf
 c
 c     Global
@@ -26,7 +26,8 @@
       include 'run.inc'
       include 'coupl.inc'
       
-c Vegas stuff
+c     Vegas stuff
+      integer         ndim
       common/tosigint/ndim
 
       real*8 sigint
@@ -175,6 +176,8 @@
          if(imode.eq.0)then
 c Don't safe the reweight information when just setting up the grids.
             doreweight=.false.
+            do_rwgt_scale=.false.
+            do_rwgt_pdf=.false.
             do kchan=1,nchans
                do i=1,ndimmax
                   do j=0,nintervals
@@ -210,10 +213,10 @@
          endif
 c
 c Setup for parton-level NLO reweighting
-         if(do_rwgt_scale.or.do_rwgt_pdf) call setup_fill_rwgt_NLOplot()
          call mint(sigint,ndim,ncall,itmax,imode,xgrid,ymax
      $        ,ymax_virt,ans,unc,chi2,nhits_in_grids)
          call topout
+         call deallocate_weight_lines
          write(*,*)'Final result [ABS]:',ans(1,0),' +/-',unc(1,0)
          write(*,*)'Final result:',ans(2,0),' +/-',unc(2,0)
          write(*,*)'chi**2 per D.o.F.:',chi2(1,0)
@@ -320,10 +323,6 @@
 c timing statistics
       include "timing_variables.inc"
       data tOLP/0.0/
-      data tFastJet/0.0/
-      data tPDF/0.0/
-      data tDSigI/0.0/
-      data tDSigR/0.0/
       data tGenPS/0.0/
       data tBorn/0.0/
       data tIS/0.0/
@@ -345,12 +344,12 @@
 
 
       double precision function sigint(xx,vegas_wgt,ifl,f)
+      use weight_lines
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'mint.inc'
       include 'nFKSconfigs.inc'
-      include 'c_weight.inc'
-      include 'reweight.inc'
       include 'run.inc'
       double precision xx(ndimmax),vegas_wgt,f(nintegrals),jac,p(0:3
      $     ,nexternal),rwgt,vol,sig,x(99),MC_int_wgt
@@ -505,6 +504,7 @@
 
       subroutine update_fks_dir(nFKS)
       implicit none
+      include 'run.inc'
       integer nFKS
       integer              nFKSprocess
       common/c_nFKSprocess/nFKSprocess
@@ -513,6 +513,7 @@
       call leshouche_inc_chooser()
       call setcuts
       call setfksfactor(.false.)
+      if (ickkw.eq.3) call configs_and_props_inc_chooser()
       return
       end
       
 
=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2017-05-09 12:41:29 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2017-08-16 21:26:30 +0000
@@ -2,6 +2,7 @@
 c**************************************************************************
 c     This is the driver for the whole calculation
 c**************************************************************************
+      use extra_weights
       implicit none
 C
 C     CONSTANTS
@@ -10,7 +11,6 @@
       parameter       (ZERO = 0d0)
       include 'nexternal.inc'
       include 'genps.inc'
-      include 'reweight.inc'
       INTEGER    ITMAX,   NCALL
 
       common/citmax/itmax,ncall
@@ -50,7 +50,9 @@
       double precision average_virtual(maxchannels),virtual_fraction(maxchannels)
       common/c_avg_virt/average_virtual,virtual_fraction
 
-      double precision weight
+      double precision weight,event_weight,inv_bias
+      character*7 event_norm
+      common /event_normalisation/event_norm
 c For MINT:
       real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
      $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
@@ -142,6 +144,8 @@
          doreweight=do_rwgt_scale.or.do_rwgt_pdf
       else
          doreweight=.false.
+         do_rwgt_scale=.false.
+         do_rwgt_pdf=.false.
       endif
       if (abrv(1:4).eq.'virt') then
          only_virt=.true.
@@ -201,6 +205,7 @@
          write (*,*) 'imode is ',imode
          call mint(sigintF,ndim,ncall,itmax,imode,xgrid,ymax,ymax_virt
      $        ,ans,unc,chi2,nhits_in_grids)
+         call deallocate_weight_lines
          open(unit=58,file='res_0',status='unknown')
          write(58,*)'Final result [ABS]:',ans(1,1),' +/-',unc(1,1)
          write(58,*)'Final result:',ans(2,1),' +/-',unc(2,1)
@@ -257,6 +262,7 @@
          write (*,*) 'imode is ',imode
          call mint(sigintF,ndim,ncall,itmax,imode,xgrid,ymax,ymax_virt
      $        ,ans,unc,chi2,nhits_in_grids)
+         call deallocate_weight_lines
          
 c If integrating the virtuals alone, we include the virtuals in
 c ans(1). Therefore, no need to have them in ans(5) and we have to set
@@ -307,7 +313,11 @@
          if (ickkw.eq.-1) putonshell=.false.
          unwgt=.true.
          open (unit=99,file='nevts',status='old',err=999)
-         read (99,*) nevts
+         if (event_norm(1:4).ne.'bias') then
+            read (99,*) nevts
+         else
+            read (99,*) nevts,event_weight
+         endif
          close(99)
          write(*,*) 'Generating ', nevts, ' events'
          if(nevts.eq.0) then
@@ -350,7 +360,11 @@
          absint=ans(1,1)+ans(5,1)
          uncer=unc(2,1)
 
-         weight=(ans(1,1)+ans(5,1))/ncall
+         if (event_norm(1:4).ne.'bias') then
+            weight=(ans(1,1)+ans(5,1))/ncall
+         else
+            weight=event_weight
+         endif
 
          if (abrv(1:3).ne.'all' .and. abrv(1:4).ne.'born' .and.
      $        abrv(1:4).ne.'virt') then
@@ -385,8 +399,13 @@
             call pick_unweight_contr(iFKS_picked)
             call update_fks_dir(iFKS_picked)
             call fill_rwgt_lines
+            if (event_norm(1:4).eq.'bias') then
+               call include_inverse_bias_wgt(inv_bias)
+               weight=event_weight*inv_bias
+            endif
             call finalize_event(x,weight,lunlhe,putonshell)
          enddo
+         call deallocate_weight_lines
          vn=-1
          call gen(sigintF,ndim,xgrid,ymax,ymax_virt,3,x,vn)
          write (*,*) 'Generation efficiencies:',x(1),x(4)
@@ -482,10 +501,6 @@
 c timing statistics
       include "timing_variables.inc"
       data tOLP/0.0/
-      data tFastJet/0.0/
-      data tPDF/0.0/
-      data tDSigI/0.0/
-      data tDSigR/0.0/
       data tGenPS/0.0/
       data tBorn/0.0/
       data tIS/0.0/
@@ -719,14 +734,7 @@
 c
 c     Here I want to set up with B.W. we map and which we don't
 c
-      dconfig = dconfig-iconfig
-      if (dconfig .eq. 0) then
-         write(*,*) 'Not subdividing B.W.'
-         lbw(0)=0
-      else
-         write(*,*) 'Error BW setting: not supported at NLO'
-         stop 1
-      endif
+      lbw(0)=0
  10   format( a)
  12   format( a,i4)
       end
@@ -739,12 +747,11 @@
 
 
       function sigintF(xx,vegas_wgt,ifl,f)
-c From dsample_fks
+      use weight_lines
       implicit none
       include 'mint.inc'
       include 'nexternal.inc'
       include 'nFKSconfigs.inc'
-      include 'c_weight.inc'
       include 'run.inc'
       logical firsttime,passcuts,passcuts_nbody,passcuts_n1body
       integer i,ifl,proc_map(0:fks_configs,0:fks_configs)
@@ -952,9 +959,11 @@
             call include_shape_in_shower_scale(p,iFKS)
          enddo
  12      continue
-
+         
 c Include PDFs and alpha_S and reweight to include the uncertainties
          call include_PDF_and_alphas
+c Include the weight from the bias_function
+         call include_bias_wgt
 c Sum the contributions that can be summed before taking the ABS value
          call sum_identical_contributions
 c Update the shower starting scale for the S-events after we have
@@ -979,7 +988,6 @@
       include 'nexternal.inc'
       include 'run.inc'
       include 'genps.inc'
-      include 'reweight_all.inc'
       include 'nFKSconfigs.inc'
       double precision lum,dlum
       external dlum
@@ -1003,19 +1011,12 @@
          write (*,*)'Using ickkw=4, include only 1 FKS dir per'/
      $        /' Born PS point (sum=0)'
       endif
-      maxproc_save=0
       do nFKSprocess=1,fks_configs
          call fks_inc_chooser()
 c Set Bjorken x's to some random value before calling the dlum() function
          xbk(1)=0.5d0
          xbk(2)=0.5d0
          lum=dlum()  ! updates IPROC
-         maxproc_save=max(maxproc_save,IPROC)
-         if (doreweight) then
-            call reweight_settozero()
-            call reweight_settozero_all(nFKSprocess*2,.true.)
-            call reweight_settozero_all(nFKSprocess*2-1,.true.)
-         endif
       enddo
       write (*,*) 'Total number of FKS directories is', fks_configs
 c For sum over identical FKS pairs, need to find the identical structures
@@ -1207,6 +1208,7 @@
 
       subroutine update_fks_dir(nFKS)
       implicit none
+      include 'run.inc'
       integer nFKS
       integer              nFKSprocess
       common/c_nFKSprocess/nFKSprocess
@@ -1215,6 +1217,7 @@
       call leshouche_inc_chooser()
       call setcuts
       call setfksfactor(.true.)
+      if (ickkw.eq.3) call configs_and_props_inc_chooser()
       return
       end
 
 
=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2017-05-09 13:42:19 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2017-08-16 21:26:30 +0000
@@ -1,9 +1,9 @@
       subroutine compute_born
 c This subroutine computes the Born matrix elements and adds its value
 c to the list of weights using the add_wgt subroutine
+      use extra_weights
       implicit none
       include 'nexternal.inc'
-      include 'reweight0.inc'
       include 'coupl.inc'
       include 'timing_variables.inc'
       double complex wgt_c(2)
@@ -31,9 +31,9 @@
       subroutine compute_nbody_noborn
 c This subroutine computes the soft-virtual matrix elements and adds its
 c value to the list of weights using the add_wgt subroutine
+      use extra_weights
       implicit none
       include 'nexternal.inc'
-      include 'reweight.inc'
       include 'coupl.inc'
       include 'run.inc'
       include 'timing_variables.inc'
@@ -94,10 +94,10 @@
       subroutine compute_real_emission(p,sudakov_damp)
 c This subroutine computes the real-emission matrix elements and adds
 c its value to the list of weights using the add_wgt subroutine
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'coupl.inc'
-      include 'reweight0.inc'
       include 'timing_variables.inc'
       double precision x,dot,f_damp,ffact,s_ev,fks_Sij,p(0:3,nexternal)
      $     ,wgt1,fx_ev,sudakov_damp
@@ -134,10 +134,10 @@
       subroutine compute_soft_counter_term(replace_MC_subt)
 c This subroutine computes the soft counter term and adds its value to
 c the list of weights using the add_wgt subroutine
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'coupl.inc'
-      include 'reweight0.inc'
       include 'timing_variables.inc'
       double precision wgt1,s_s,fks_Sij,fx_s,zero,replace_MC_subt,g22
       parameter (zero=0d0)
@@ -185,11 +185,11 @@
       subroutine compute_collinear_counter_term(replace_MC_subt)
 c This subroutine computes the collinear counter term and adds its value
 c to the list of weights using the add_wgt subroutine
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'coupl.inc'
       include 'fks_powers.inc'
-      include 'reweight.inc'
       include 'timing_variables.inc'
       double precision zero,one,s_c,fks_Sij,fx_c,deg_xi_c,deg_lxi_c,wgt1
      &     ,wgt3,g22,replace_MC_subt
@@ -248,10 +248,10 @@
       subroutine compute_soft_collinear_counter_term(replace_MC_subt)
 c This subroutine computes the soft-collinear counter term and adds its
 c value to the list of weights using the add_wgt subroutine
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'coupl.inc'
-      include 'reweight.inc'
       include 'fks_powers.inc'
       include 'timing_variables.inc'
       double precision zero,one,s_sc,fks_Sij,fx_sc,wgt1,wgt3,deg_xi_sc
@@ -314,6 +314,7 @@
       end
 
       subroutine compute_MC_subt_term(p,gfactsf,gfactcl,probne)
+      use extra_weights
       implicit none
 c This subroutine computes the MonteCarlo subtraction terms and adds
 c their values to the list of weights using the add_wgt subroutine. It
@@ -324,7 +325,6 @@
       include 'nexternal.inc'
       include 'madfks_mcatnlo.inc'
       include 'timing_variables.inc'
-      include 'reweight.inc'
       include 'coupl.inc'
       integer nofpartners,i
       double precision p(0:3,nexternal),gfactsf,gfactcl,probne,x,dot
@@ -613,13 +613,13 @@
 c Compute all the relevant prefactors for the Born and the soft-virtual,
 c i.e. all the nbody contributions. Also initialises the plots and
 c bpower.
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'genps.inc'
-      include 'reweight0.inc'
       include 'timing_variables.inc'
-      double precision pi,unwgtfun,vegas_wgt,enhance,xnoborn_cnt,xtot
+      double precision pi,vegas_wgt,enhance,xnoborn_cnt,xtot
      $     ,bpower,cpower,tiny
       data xnoborn_cnt /0d0/
       integer inoborn_cnt,i
@@ -694,7 +694,7 @@
          if(wgtcpower.ne.cpowerinput.and.dabs(cpower+1d0).gt.tiny)then
            write(*,*)'Inconsistency in the computation of cpower',
      #               wgtcpower,cpowerinput
-           write(*,*)'Check value in reweight0.inc'
+           write(*,*)'Check value in extra_weights.f'
            stop
          endif
          firsttime=.false.
@@ -732,11 +732,10 @@
             enhance=0d0
          endif
       endif
-      call unweight_function(p_born,unwgtfun)
       call set_cms_stuff(0)
 c f_* multiplication factors for Born and nbody
       f_b=jac_cnt(0)*xinorm_ev/(min(xiimax_ev,xiBSVcut_used)*shat/(16
-     $     *pi**2))*enhance*unwgtfun *fkssymmetryfactorBorn*vegas_wgt
+     $     *pi**2))*enhance*fkssymmetryfactorBorn*vegas_wgt
       f_nb=f_b
       call cpu_time(tAfter)
       tf_nb=tf_nb+(tAfter-tBefore)
@@ -753,7 +752,7 @@
       include 'fks_powers.inc'
       include 'coupl.inc'
       include 'timing_variables.inc'
-      double precision unwgtfun,vegas_wgt,enhance,xnoborn_cnt,xtot
+      double precision vegas_wgt,enhance,xnoborn_cnt,xtot
      &     ,prefact,prefact_cnt_ssc,prefact_deg,prefact_c,prefact_coll
      &     ,jac_ev,pi,prefact_cnt_ssc_c,prefact_coll_c,prefact_deg_slxi
      &     ,prefact_deg_sxi,zero
@@ -845,11 +844,10 @@
             enhance=0d0
          endif
       endif
-      call unweight_function(p_born,unwgtfun)
       prefact=xinorm_ev/xi_i_fks_ev/(1-y_ij_fks_ev)
 
 c f_* multiplication factors for real-emission, soft counter, ... etc.       
-      f_r=prefact*jac_ev*enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+      f_r=prefact*jac_ev*enhance*fkssymmetryfactor*vegas_wgt
       f_MC_S=f_r
       f_MC_H=f_r
       if (.not.nocntevents) then
@@ -857,9 +855,9 @@
      &        log(xicut_used/min(xiimax_ev,xiScut_used))/(1
      &        -y_ij_fks_ev)
          f_s=(prefact+prefact_cnt_ssc)*jac_cnt(0)*enhance
-     $        *unwgtfun*fkssymmetryfactor*vegas_wgt
+     $        *fkssymmetryfactor*vegas_wgt
          f_s_MC_S=prefact*jac_cnt(0)*enhance
-     $        *unwgtfun*fkssymmetryfactor*vegas_wgt
+     $        *fkssymmetryfactor*vegas_wgt
          f_s_MC_H=f_s_MC_S
 
          if (pmass(j_fks).eq.0d0) then
@@ -867,9 +865,9 @@
             prefact_coll=xinorm_cnt(1)/xi_i_fks_cnt(1)*log(delta_used
      $           /deltaS)/deltaS
             f_c=(prefact_c+prefact_coll)*jac_cnt(1)
-     $           *enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+     $           *enhance*fkssymmetryfactor*vegas_wgt
             f_c_MC_S=prefact_c*jac_cnt(1)
-     $           *enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+     $           *enhance*fkssymmetryfactor*vegas_wgt
             f_c_MC_H=f_c_MC_S
 
             call set_cms_stuff(1)
@@ -881,12 +879,12 @@
      $           *log(xicut_used/min(xiimax_cnt(1),xiScut_used))
      $           *log(delta_used/deltaS)/deltaS
             f_dc=jac_cnt(1)*prefact_deg/(shat/(32*pi**2))*enhance
-     $           *unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+     $           *fkssymmetryfactorDeg*vegas_wgt
             f_sc=(prefact_c+prefact_coll+prefact_cnt_ssc_c
-     &           +prefact_coll_c)*jac_cnt(2)*enhance*unwgtfun
+     &           +prefact_coll_c)*jac_cnt(2)*enhance
      &           *fkssymmetryfactorDeg*vegas_wgt
             f_sc_MC_S=prefact_c*jac_cnt(2)
-     $           *enhance*unwgtfun*fkssymmetryfactor*vegas_wgt
+     $           *enhance*fkssymmetryfactor*vegas_wgt
             f_sc_MC_H=f_sc_MC_S
 
             call set_cms_stuff(2)
@@ -898,13 +896,13 @@
      &           -log(min(xiimax_cnt(1),xiScut_used))**2 )*1/(2.d0
      &           *deltaS)
             f_dsc(1)=prefact_deg*jac_cnt(2)/(shat/(32*pi**2))*enhance
-     &           *unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+     &           *fkssymmetryfactorDeg*vegas_wgt
             f_dsc(2)=prefact_deg_sxi*jac_cnt(2)/(shat/(32*pi**2))
-     &           *enhance*unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+     &           *enhance*fkssymmetryfactorDeg*vegas_wgt
             f_dsc(3)=prefact_deg_slxi*jac_cnt(2)/(shat/(32*pi**2))
-     &           *enhance*unwgtfun*fkssymmetryfactorDeg*vegas_wgt
+     &           *enhance*fkssymmetryfactorDeg*vegas_wgt
             f_dsc(4)=( prefact_deg+prefact_deg_sxi )*jac_cnt(2)/(shat
-     &           /(32*pi**2))*enhance*unwgtfun*fkssymmetryfactorDeg
+     &           /(32*pi**2))*enhance*fkssymmetryfactorDeg
      &           *vegas_wgt
          else
             f_c=0d0
@@ -940,7 +938,7 @@
 
       
       subroutine add_wgt(type,wgt1,wgt2,wgt3)
-c Adds a contribution to the list in c_weight.inc. 'type' sets the type
+c Adds a contribution to the list in weight_lines. 'type' sets the type
 c of the contribution and wgt1..wgt3 are the coefficients multiplying
 c the logs. The arguments are:
 c     type=1 : real-emission
@@ -997,7 +995,7 @@
 c        corresponding to this contribution: wgt_ME_tree. This weight does
 c        include the 'ngluon' correction factor for the Born.
 c
-c Not set in this subroutine, but included in the c_weights common block
+c Not set in this subroutine, but included in the weight_lines module
 c are the
 c     wgts(iwgt,icontr) : weights including scale/PDFs/logs. These are
 c        normalised so that they can be used directly to compute cross
@@ -1016,16 +1014,16 @@
 c     parton_iproc(iproc,icontr) : value of the PDF for the iproc
 c        contribution
 c     parton_pdg(nexternal,iproc,icontr) : value of the PDG codes for
-c        the iproc contribution
+c     the iproc contribution
+      use weight_lines
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'genps.inc'
       include 'coupl.inc'
       include 'fks_info.inc'
-      include 'c_weight.inc'
       include 'q_es.inc'
-      include 'reweight0.inc'
       integer type,i,j
       double precision wgt1,wgt2,wgt3
       integer              nFKSprocess
@@ -1054,11 +1052,7 @@
       if (wgt2.ne.wgt2) return
       if (wgt3.ne.wgt3) return
       icontr=icontr+1
-      if (icontr.gt.max_contr) then
-         write (*,*) 'ERROR in add_wgt: too many contributions'
-     &        ,max_contr
-         stop 1
-      endif
+      call weight_lines_allocated(nexternal,icontr,max_wgt,max_iproc)
       itype(icontr)=type
       wgt(1,icontr)=wgt1
       wgt(2,icontr)=wgt2
@@ -1133,12 +1127,12 @@
       end
 
       subroutine include_veto_multiplier
+      use weight_lines
+      use extra_weights
       implicit none
 c Multiply all the weights by the NNLL-NLO jet veto Sudakov factors,
 c i.e., the term on the 2nd line of Eq.(20) of arXiv:1412.8408.
       include 'nexternal.inc'
-      include 'c_weight.inc'
-      include 'reweight.inc'
       integer i,j
       if (H1_factor_virt.ne.0d0) then
          call compute_veto_multiplier(H1_factor_virt,1d0,1d0
@@ -1158,10 +1152,10 @@
 c dependence and saves the weights in the wgts() array. The weights in
 c this array are now correctly normalised to compute the cross section
 c or to fill histograms.
+      use weight_lines
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
-      include 'c_weight.inc'
       include 'coupl.inc'
       include 'timing_variables.inc'
       include 'genps.inc'
@@ -1191,14 +1185,11 @@
          q2fact(2)=mu2_f
 c call the PDFs
          xlum = dlum()
-         if (iproc.gt.max_iproc) then
-            write (*,*) 'ERROR iproc too large',iproc,max_iproc
-            stop 1
-         endif
+         iwgt=1
+         call weight_lines_allocated(nexternal,max_contr,iwgt,iproc)
 c set_pdg_codes fills the niproc, parton_iproc, parton_pdg and parton_pdg_uborn
          call set_pdg_codes(iproc,pd,nFKSprocess,i)
 c iwgt=1 is the central value (i.e. no scale/PDF reweighting).
-         iwgt=1
          wgt_wo_pdf=(wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q) + wgt(3,i)
      &        *log(mu2_f/mu2_q))*g_strong(i)**QCDpower(i)
      &        *rwgt_muR_dep_fac(sqrt(mu2_r),sqrt(mu2_r))
@@ -1221,12 +1212,123 @@
       return
       end
 
+      subroutine include_bias_wgt
+c Include the weight from the bias_wgt_function to all the contributions
+c in icontr. This only changes the weight of the central value (after
+c inclusion of alphaS and parton luminosity). Both for 'wgts(1,icontr)'
+c as well as the the 'parton_iproc(1:niproc(icontr),icontr)', since
+c these are the ones used in MINT as well as for unweighting. Also the
+c 'virt_wgt_mint' and 'born_wgt_mint' are updated. Furthermore, to
+c include the weight also in the 'wgt' array that contain the
+c coefficients for PDF and scale computations. 
+      use weight_lines
+      implicit none
+      integer i,j
+      double precision bias
+      character*7 event_norm
+      common /event_normalisation/event_norm
+      double precision           virt_wgt_mint,born_wgt_mint
+      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
+c Set the bias_wgt to 1 in case we do not have to do any biassing
+      if (event_norm(1:4).ne.'bias') then
+         do i=1,icontr
+            bias_wgt(i)=1d0
+         enddo
+         return
+      endif
+c loop over all contributions
+      do i=1,icontr
+         if (itype(i).eq.1) then
+            ! use (n+1)-body momenta for the real emission. Pick the
+            ! first IPROC for parton PDGs.
+            call bias_weight_function(momenta_m(0,1,2,i),parton_pdg(1,1
+     $           ,i),bias)
+         else
+            ! use n-body momenta for all the other contributions. Pick
+            ! the first IPROC for parton PDGs.
+            call bias_weight_function(momenta_m(0,1,1,i),parton_pdg(1,1
+     $           ,i),bias)
+         endif
+         bias_wgt(i)=bias
+c Update the weights:
+         wgts(1,i)=wgts(1,i)*bias_wgt(i)
+         do j=1,niproc(i)
+            parton_iproc(j,i)=parton_iproc(j,i)*bias_wgt(i)
+         enddo
+         do j=1,3
+            wgt(j,i)=wgt(j,i)*bias_wgt(i)
+         enddo
+         if (itype(i).eq.14) then
+            virt_wgt_mint=virt_wgt_mint*bias_wgt(i)
+            born_wgt_mint=born_wgt_mint*bias_wgt(i)
+         endif
+      enddo
+      return
+      end
+
+      subroutine include_inverse_bias_wgt(inv_bias)
+c Update the inverse of the bias in the event weight. All information in
+c the rwgt_lines is NOT updated.
+      use weight_lines
+      use extra_weights
+      implicit none
+      include 'genps.inc'
+      include 'nFKSconfigs.inc'
+      integer i,ict,ipr,ii
+      double precision wgt_num,wgt_denom,inv_bias
+      character*7 event_norm
+      common /event_normalisation/event_norm
+      integer iproc_save(fks_configs),eto(maxproc,fks_configs)
+     $     ,etoi(maxproc,fks_configs),maxproc_found
+      common/cproc_combination/iproc_save,eto,etoi,maxproc_found
+      logical         Hevents
+      common/SHevents/Hevents
+      if (event_norm(1:4).ne.'bias') then
+         inv_bias=1d0
+         return
+      endif
+      wgt_num=0d0
+      wgt_denom=0d0
+      do i=1,icontr_sum(0,icontr_picked)
+         ict=icontr_sum(i,icontr_picked)
+         if (bias_wgt(ict).eq.0d0) then
+            write (*,*) "ERROR in include_inverse_bias_wgt: "/
+     $           /"bias_wgt is equal to zero",ict,bias_wgt
+            stop 1
+         endif
+c for all the rwgt_lines, remove the bias-wgt contribution from the
+c weights there. Note that the wgtref (also written in the event file)
+c keeps its contribution from the bias_wgt.
+         if (.not. Hevents) then
+            ipr=eto(etoi(iproc_picked,nFKS(ict)),nFKS(ict))
+            do ii=1,iproc_save(nFKS(ict))
+               if (eto(ii,nFKS(ict)).ne.ipr) cycle
+               wgt_denom=wgt_denom+parton_iproc(ii,ict)
+               wgt_num=wgt_num+parton_iproc(ii,ict)/bias_wgt(ict)
+            enddo
+         else
+            ipr=iproc_picked
+            wgt_denom=wgt_denom+parton_iproc(ipr,ict)
+            wgt_num=wgt_num+parton_iproc(ipr,ict)/bias_wgt(ict)
+         endif
+      enddo
+      if (abs((wgtref-wgt_denom)/(wgtref+wgt_denom)).gt.1d-10) then
+         write (*,*) "ERROR in include_inverse_bias_wgt: "/
+     $        /"reference weight not equal to recomputed weight",wgtref
+     $        ,wgt_denom
+         stop 1
+      endif
+c update the event weight to be written in the file
+      inv_bias=wgt_num/wgt_denom
+      return
+      end
+      
 
       subroutine set_pdg_codes(iproc,pd,iFKS,ict)
+      use weight_lines
       implicit none
       include 'nexternal.inc'
       include 'genps.inc'
-      include 'c_weight.inc'
       include 'fks_info.inc'
       integer j,k,iproc,ict,iFKS
       double precision  pd(0:maxproc),conv
@@ -1274,14 +1376,13 @@
       
       
       subroutine reweight_scale
-c Use the saved c_weight info to perform scale reweighting. Extends the
+c Use the saved weight_lines info to perform scale reweighting. Extends the
 c wgts() array to include the weights.
+      use weight_lines
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
-      include 'c_weight.inc'
-      include 'reweight.inc'
-      include 'reweightNLO.inc'
       include 'timing_variables.inc'
       integer i,kr,kf,iwgt_save,dd
       double precision xlum(maxscales),dlum,pi,mu2_r(maxscales),c_mu2_r
@@ -1298,7 +1399,7 @@
       if (icontr.eq.0) return
 c currently we have 'iwgt' weights in the wgts() array.
       iwgt_save=iwgt
-c loop over all the contributions in the c_weights common block
+c loop over all the contributions in the weight lines module
       do i=1,icontr
          iwgt=iwgt_save
          nFKSprocess=nFKS(i)
@@ -1328,11 +1429,8 @@
                do kr=1,nint(scalevarR(0))
                   if ((.not. lscalevar(dd)) .and. kr.ne.1) exit
                   iwgt=iwgt+1   ! increment the iwgt for the wgts() array
-                  if (iwgt.gt.max_wgt) then
-                     write (*,*) 'ERROR too many weights in '/
-     $                    /'reweight_scale',iwgt,max_wgt
-                     stop 1
-                  endif
+                  call weight_lines_allocated(nexternal,max_contr,iwgt
+     $                 ,max_iproc)
 c add the weights to the array
                   wgts(iwgt,i)=xlum(kf) * (wgt(1,i)+wgt(2,i)
      $                 *log(mu2_r(kr)/mu2_q)+wgt(3,i)*log(mu2_f(kf)
@@ -1349,15 +1447,14 @@
       end
 
       subroutine reweight_scale_NNLL
-c Use the saved c_weight info to perform scale reweighting. Extends the
+c Use the saved weight lines info to perform scale reweighting. Extends the
 c wgts() array to include the weights. Special for the NNLL+NLO jet-veto
 c computations (ickkw.eq.-1).
+      use weight_lines
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
-      include 'c_weight.inc'
-      include 'reweight.inc'
-      include 'reweightNLO.inc'
       include 'timing_variables.inc'
       integer i,ks,kh,iwgt_save
       double precision xlum(maxscales),dlum,pi,mu2_r(maxscales)
@@ -1394,7 +1491,7 @@
             endif
          enddo
       enddo
-c loop over all the contributions in the c_weights common block
+c loop over all the contributions in the weight lines module
       do i=1,icontr
          iwgt=iwgt_save
          nFKSprocess=nFKS(i)
@@ -1414,11 +1511,8 @@
                q2fact(2)=mu2_f(ks)
                xlum(ks) = dlum()
                iwgt=iwgt+1      ! increment the iwgt for the wgts() array
-               if (iwgt.gt.max_wgt) then
-                  write (*,*) 'ERROR too many weights in reweight_scale'
-     &                 ,iwgt,max_wgt
-                  stop 1
-               endif
+               call weight_lines_allocated(nexternal,max_contr,iwgt
+     $              ,max_iproc)
 c add the weights to the array
                if (itype(i).ne.7) then
                   wgts(iwgt,i)=xlum(ks) * (wgt(1,i)+wgt(2,i)
@@ -1445,14 +1539,13 @@
       end
 
       subroutine reweight_pdf
-c Use the saved c_weight info to perform PDF reweighting. Extends the
+c Use the saved weight_lines info to perform PDF reweighting. Extends the
 c wgts() array to include the weights.
+      use weight_lines
+      use extra_weights
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
-      include 'c_weight.inc'
-      include 'reweight.inc'
-      include 'reweightNLO.inc'
       include 'timing_variables.inc'
       integer n,i,nn
       double precision xlum,dlum,pi,mu2_r,mu2_f,mu2_q,rwgt_muR_dep_fac,g
@@ -1470,11 +1563,8 @@
 c allows for better caching of the PDFs
          do n=0,nmemPDF(nn)
             iwgt=iwgt+1
-            if (iwgt.gt.max_wgt) then
-               write (*,*) 'ERROR too many weights in reweight_pdf',iwgt
-     &              ,max_wgt
-               stop 1
-            endif
+            call weight_lines_allocated(nexternal,max_contr,iwgt
+     $           ,max_iproc)
             call InitPDFm(nn,n)
             do i=1,icontr
                nFKSprocess=nFKS(i)
@@ -1508,9 +1598,9 @@
 c that there is an unique PS configuration: at most one Born, one real
 c and one set of counter terms. Among other things, this means that one
 c must do MC over FKS directories.
+      use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
       include 'appl_common.inc'
       include 'nFKSconfigs.inc'
       include 'genps.inc'
@@ -1606,9 +1696,9 @@
 c the pdg_uborn (the PDG codes for the underlying Born process) the PDG
 c codes of i_fks and j_fks are combined to give the PDG code of the
 c mother and the extra (n+1) parton is given the PDG code of the gluon.
+      use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
       include 'fks_info.inc'
       include 'genps.inc'
       integer k,ict,iFKS
@@ -1651,9 +1741,9 @@
       subroutine get_wgt_nbody(sig)
 c Sums all the central weights that contribution to the nbody cross
 c section
+      use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
       double precision sig
       integer i
       sig=0d0
@@ -1670,9 +1760,9 @@
       subroutine get_wgt_no_nbody(sig)
 c Sums all the central weights that contribution to the cross section
 c excluding the nbody contributions.
+      use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
       double precision sig
       integer i
       sig=0d0
@@ -1688,20 +1778,20 @@
 
       subroutine fill_plots
 c Calls the analysis routine (which fill plots) for all the
-c contributions in the c_weight common block. Instead of really calling
+c contributions in the weight_lines module. Instead of really calling
 c it for all, it first checks if weights can be summed (i.e. they have
 c the same PDG codes and the same momenta) before calling the analysis
 c to greatly reduce the calls to the analysis routines.
+      use weight_lines
+      use extra_weights
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
-      include 'reweight0.inc'
       include 'timing_variables.inc'
       integer i,ii,j,max_weight
       logical momenta_equal,pdg_equal
       external momenta_equal,pdg_equal
-      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
-      double precision www(max_weight)
+      double precision,allocatable :: www(:)
+      save max_weight
       call cpu_time(tBefore)
       if (icontr.eq.0) return
 c fill the plots_wgts. Check if we can sum weights together before
@@ -1739,10 +1829,14 @@
       enddo
       do i=1,icontr
          if (plot_wgts(1,i).ne.0d0) then
-            if (iwgt.gt.max_weight) then
-               write (*,*) 'ERROR too many weights in fill_plots',iwgt
-     &              ,max_weight
-               stop 1
+            if (.not.allocated(www)) then
+               allocate(www(iwgt))
+               max_weight=iwgt
+            elseif(iwgt.ne.max_weight) then
+               write (*,*) 'Error in fill_plots (fks_singular.f): '/
+     $              /'number of weights should not vary between PS '/
+     $              /'points',iwgt,max_weight
+               stop
             endif
             do j=1,iwgt
                www(j)=plot_wgts(j,i)
@@ -1758,19 +1852,32 @@
 
       subroutine fill_mint_function(f)
 c Fills the function that is returned to the MINT integrator
+      use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
       include 'mint.inc'
       integer i
-      double precision f(nintegrals),sigint
+      double precision f(nintegrals),sigint,bias
       double precision           virt_wgt_mint,born_wgt_mint
       common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
       double precision virtual_over_born
       common /c_vob/   virtual_over_born
+      character*7 event_norm
+      common /event_normalisation/event_norm
       sigint=0d0
       do i=1,icontr
-         sigint=sigint+wgts(1,i)
+         if (event_norm(1:4).ne.'bias') then
+            sigint=sigint+wgts(1,i)
+         else
+            if (itype(i).eq.1) then
+               call bias_weight_function(momenta_m(0,1,2,i),parton_pdg(1
+     $              ,1,i),bias)
+            else
+               call bias_weight_function(momenta_m(0,1,1,i),parton_pdg(1
+     $              ,1,i),bias)
+            endif
+            sigint=sigint+wgts(1,i)*bias
+         endif
       enddo
       f(1)=abs(sigint)
       f(2)=sigint
@@ -1786,10 +1893,10 @@
 c Includes the shape function from the MC counter terms in the shower
 c starting scale. This function needs to be called (at least) once per
 c FKS configuration that is included in the current PS point.
+      use weight_lines
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
-      include 'c_weight.inc'
       include 'nFKSconfigs.inc'
       integer i,iFKS,Hevents,izero,mohdr
       double precision ddum(6),p(0:3,nexternal)
@@ -1847,9 +1954,9 @@
 c the MC counter terms for the H-events FKS configuration by FKS
 c configuration, while for the S-events also contributions from the
 c various FKS configurations can be summed together.
+      use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
       include 'genps.inc'
       include 'nFKSconfigs.inc'
       include 'fks_info.inc'
@@ -1953,9 +2060,9 @@
 c necessarily the same for all of these summed FKS configurations). Take
 c the weighted average over the FKS configurations as the shower scale
 c for the summed contribution.
+      use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
       include 'nFKSconfigs.inc'
       integer i,j,ict
       double precision tmp_wgt(fks_configs),showerscale(fks_configs)
@@ -2009,9 +2116,9 @@
       subroutine fill_mint_function_NLOPS(f,n1body_wgt)
 c Fills the function that is returned to the MINT integrator. Depending
 c on the imode we should or should not include the virtual corrections.
+      use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
       include 'mint.inc'
       integer i,j,ict
       double precision f(nintegrals),sigint,sigint1,sigint_ABS
@@ -2099,9 +2206,9 @@
       subroutine pick_unweight_contr(iFKS_picked)
 c Randomly pick (weighted by the ABS values) the contribution to a given
 c PS point that should be written in the event file.
+      use weight_lines
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
       include 'genps.inc'
       include 'nFKSconfigs.inc'
       include 'fks_info.inc'
@@ -2179,19 +2286,21 @@
 c Fills the lines, n_ctr_str, to be written in an event file with the
 c (internal) information to perform scale and/or PDF reweighting. All
 c information is available in each line to do the reweighting, apart
-c from the momenta: these are put in the momenta_str_l() array, and a
+c from the momenta: these are put in the momenta_str() array, and a
 c label in each of the n_ctr_str refers to a corresponding set of
-c momenta in the momenta_str_l() array.
+c momenta in the momenta_str() array.
+      use weight_lines
+      use extra_weights
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
-      include 'reweight0.inc'
       include 'genps.inc'
       include 'nFKSconfigs.inc'
       include 'fks_info.inc'
       integer k,i,ii,j,jj,ict,ipr,momenta_conf(2)
       logical momenta_equal,found
-      double precision conv,momenta_str_l(0:3,nexternal,max_n_ctr)
+      double precision conv
+      double precision,allocatable :: temp3(:,:,:)
+      character(len=1024),allocatable :: ctemp(:)
       external momenta_equal
       character*512 procid,str_temp
       parameter (conv=389379660d0) ! conversion to picobarns
@@ -2200,15 +2309,17 @@
       common/cproc_combination/iproc_save,eto,etoi,maxproc_found
       logical         Hevents
       common/SHevents/Hevents
+      if (.not.allocated(momenta_str)) allocate(momenta_str(0:3
+     $     ,max_mext,max_mom_str))
       wgtref=unwgt(iproc_picked,icontr_picked)
       n_ctr_found=0
       n_mom_conf=0
 c Loop over all the contributions in the picked contribution (the latter
-c is chosen in the pick_unweight_cont() subroutine)
+c is chosen in the pick_unweight_contr() subroutine)
       do i=1,icontr_sum(0,icontr_picked)
          ict=icontr_sum(i,icontr_picked)
 c Check if the current set of momenta are already available in the
-c momenta_str_l array. If not, add it.
+c momenta_str array. If not, add it.
          found=.false.
          do k=1,2
             do j=1,n_mom_conf
@@ -2216,7 +2327,7 @@
                   momenta_conf(k)=0
                   cycle
                endif
-               if (momenta_equal(momenta_str_l(0,1,j),
+               if (momenta_equal(momenta_str(0,1,j),
      &                           momenta_m(0,1,k,ict))) then
                   momenta_conf(k)=j
                   found=.true.
@@ -2225,12 +2336,20 @@
             enddo
             if (.not. found) then
                n_mom_conf=n_mom_conf+1
+               if (n_mom_conf.gt.max_mom_str .or. nexternal.gt.max_mext)
+     $              then
+                  allocate(temp3(0:3,max(nexternal,max_mext)
+     $                              ,max(n_mom_conf,max_mom_str)))
+                  temp3(0:3,1:min(nexternal,max_mext)
+     $                     ,1:min(max_mom_str,n_mom_conf))=momenta_str
+                  call move_alloc(temp3,momenta_str)
+                  max_mom_str=max(n_mom_conf,max_mom_str)
+                  max_mext=max(nexternal,max_mext)
+               endif
                do ii=1,nexternal
                   do jj=0,3
                      momenta_str(jj,ii,n_mom_conf)=
      &                                      momenta_m(jj,ii,k,ict)
-                     momenta_str_l(jj,ii,n_mom_conf)=
-     &                                      momenta_m(jj,ii,k,ict)
                   enddo
                enddo
                momenta_conf(k)=n_mom_conf
@@ -2244,6 +2363,15 @@
                if (eto(ii,nFKS(ict)).ne.ipr) cycle
                n_ctr_found=n_ctr_found+1
 
+               if (.not.allocated(n_ctr_str))
+     $              allocate(n_ctr_str(max_n_ctr))
+               if (n_ctr_found.gt.max_n_ctr) then
+                  allocate(ctemp(n_ctr_found))
+                  ctemp(1:max_n_ctr)=n_ctr_str
+                  call move_alloc(ctemp,n_ctr_str)
+                  max_n_ctr=n_ctr_found
+               endif
+               
                if (nincoming.eq.2) then
                   write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
      &                 (wgt(j,ict)*conv,j=1,3),(wgt_me_tree(j,ict),j=1,2),
@@ -2264,8 +2392,7 @@
      &              trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &              //trim(adjustl(procid))
 
-               write (str_temp,
-     &                    '(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12)')
+               write (str_temp,30)
      &              QCDpower(ict),
      &              (bjx(j,ict),j=1,2),
      &              (scales2(j,ict),j=1,3),
@@ -2276,7 +2403,8 @@
      &              fks_i_d(nFKS(ict)),
      &              fks_j_d(nFKS(ict)),
      &              parton_pdg_uborn(fks_j_d(nFKS(ict)),ii,ict),
-     &              parton_iproc(ii,ict)
+     &              parton_iproc(ii,ict),
+     &              bias_wgt(ict)
                n_ctr_str(n_ctr_found) =
      &              trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &              //trim(adjustl(str_temp))
@@ -2286,6 +2414,15 @@
             ipr=iproc_picked
             n_ctr_found=n_ctr_found+1
 
+            if (.not.allocated(n_ctr_str))
+     $           allocate(n_ctr_str(max_n_ctr))
+            if (n_ctr_found.gt.max_n_ctr) then
+               allocate(ctemp(n_ctr_found))
+               ctemp(1:max_n_ctr)=n_ctr_str
+               call move_alloc(ctemp,n_ctr_str)
+               max_n_ctr=n_ctr_found
+            endif
+
             if (nincoming.eq.2) then
                write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
      &              (wgt(j,ict)*conv,j=1,3),(wgt_me_tree(j,ict),j=1,2),
@@ -2306,7 +2443,7 @@
      &           trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &           //trim(adjustl(procid))
 
-            write (str_temp,'(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12)')
+            write (str_temp,30)
      &           QCDpower(ict),
      &           (bjx(j,ict),j=1,2),
      &           (scales2(j,ict),j=1,3),
@@ -2317,18 +2454,17 @@
      &           fks_i_d(nFKS(ict)),
      &           fks_j_d(nFKS(ict)),
      &           parton_pdg_uborn(fks_j_d(nFKS(ict)),ipr,ict),
-     &           parton_iproc(ipr,ict)
+     &           parton_iproc(ipr,ict),
+     &           bias_wgt(ict)
             n_ctr_str(n_ctr_found) =
      &           trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &           //trim(adjustl(str_temp))
 
 
          endif
-         if (n_ctr_found.ge.max_n_ctr) then
-            write (*,*) 'ERROR: too many contributions in <rwgt>'
-            stop1
-         endif
       enddo
+      return
+ 30   format(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12,1x,d18.12)
       end
       
       
@@ -2878,9 +3014,6 @@
       double complex xij_aor
       common/cxij_aor/xij_aor
 
-      logical rotategranny
-      common/crotategranny/rotategranny
-
       double precision cthbe,sthbe,cphibe,sphibe
       common/cbeangles/cthbe,sthbe,cphibe,sphibe
 
@@ -2918,20 +3051,7 @@
       E_i_fks = p(0,i_fks)
       z = 1d0 - E_i_fks/(E_i_fks+E_j_fks)
       t = z * shat/4d0
-      if(rotategranny .and. nexternal-1.ne.3 .and. nincoming.eq.2)then
-c Exclude 2->1 (at the Born level) processes: matrix elements are
-c independent of the PS point, but non-zero helicity configurations
-c might flip when rotating the momenta.
-        do i=1,nexternal-1
-          call trp_rotate_invar(p_born(0,i),p_born_rot(0,i),
-     #                          cthbe,sthbe,cphibe,sphibe)
-        enddo
-        CalculatedBorn=.false.
-        call sborn(p_born_rot,wgt1)
-        CalculatedBorn=.false.
-      else
-        call sborn(p_born,wgt1)
-      endif
+      call sborn(p_born,wgt1)
       call AP_reduced(j_type,i_type,t,z,ap)
       if (abs(j_type).eq.3 .and. i_type.eq.8) then
          Q=0d0
@@ -2945,10 +3065,6 @@
                pi(i)=p_i_fks_ev(i)
                pj(i)=p(i,j_fks)
             enddo
-            if(rotategranny)then
-              call trp_rotate_invar(pi,pi,cthbe,sthbe,cphibe,sphibe)
-              call trp_rotate_invar(pj,pj,cthbe,sthbe,cphibe,sphibe)
-            endif
             CALL IXXXSO(pi ,ZERO ,+1,+1,W1)        
             CALL OXXXSO(pj ,ZERO ,-1,+1,W2)        
             CALL IXXXSO(pi ,ZERO ,-1,+1,W3)        
@@ -2963,13 +3079,8 @@
          endif
 c Insert the extra factor due to Madgraph convention for polarization vectors
          imother_fks=min(i_fks,j_fks)
-         if(rotategranny)then
-           call getaziangles(p_born_rot(0,imother_fks),
-     #                       cphi_mother,sphi_mother)
-         else
-           call getaziangles(p_born(0,imother_fks),
-     #                       cphi_mother,sphi_mother)
-         endif
+         call getaziangles(p_born(0,imother_fks),
+     #                     cphi_mother,sphi_mother)
          wgt1(2) = -(cphi_mother-ximag*sphi_mother)**2 *
      #             wgt1(2) * azifact
          call Qterms_reduced_timelike(j_type, i_type, t, z, Q)
@@ -3040,24 +3151,7 @@
 c Thus, an extra factor z (implicit in the flux of the reduced Born
 c in FKS) has to be inserted here
       t = z*shat/4d0
-      if(j_fks.eq.2 .and. nexternal-1.ne.3 .and. nincoming.eq.2)then
-c Rotation according to innerpin.m. Use rotate_invar() if a more 
-c general rotation is needed.
-c Exclude 2->1 (at the Born level) processes: matrix elements are
-c independent of the PS point, but non-zero helicity configurations
-c might flip when rotating the momenta.
-        do i=1,nexternal-1
-          p_born_rot(0,i)=p_born(0,i)
-          p_born_rot(1,i)=-p_born(1,i)
-          p_born_rot(2,i)=p_born(2,i)
-          p_born_rot(3,i)=-p_born(3,i)
-        enddo
-        CalculatedBorn=.false.
-        call sborn(p_born_rot,wgt1)
-        CalculatedBorn=.false.
-      else
-        call sborn(p_born,wgt1)
-      endif
+      call sborn(p_born,wgt1)
       call AP_reduced(m_type,i_type,t,z,ap)
       if (abs(m_type).eq.3) then
          Q=0d0
@@ -3071,14 +3165,6 @@
                pi(i)=p_i_fks_ev(i)
                pj(i)=p(i,j_fks)
             enddo
-            if(j_fks.eq.2 .and. nincoming.eq.2)then
-c Rotation according to innerpin.m. Use rotate_invar() if a more 
-c general rotation is needed
-               pi(1)=-pi(1)
-               pi(3)=-pi(3)
-               pj(1)=-pj(1)
-               pj(3)=-pj(3)
-            endif
             CALL IXXXSO(pi ,ZERO ,+1,+1,W1)        
             CALL OXXXSO(pj ,ZERO ,-1,+1,W2)        
             CALL IXXXSO(pi ,ZERO ,-1,+1,W3)        
@@ -3092,13 +3178,8 @@
             azifact=Wij_angle/Wij_recta
          endif
 c Insert the extra factor due to Madgraph convention for polarization vectors
-         if(j_fks.eq.2 .and. nincoming.eq.2)then
-           cphi_mother=-1.d0
-           sphi_mother=0.d0
-         else
-           cphi_mother=1.d0
-           sphi_mother=0.d0
-         endif
+         cphi_mother=1.d0
+         sphi_mother=0.d0
          wgt1(2) = -(cphi_mother+ximag*sphi_mother)**2 *
      #             wgt1(2) * dconjg(azifact)
          call Qterms_reduced_spacelike(m_type, i_type, t, z, Q)
@@ -3508,13 +3589,13 @@
 
       subroutine sreal_deg(p,xi_i_fks,y_ij_fks,
      #                     collrem_xi,collrem_lxi)
+      use extra_weights
       implicit none
       include "genps.inc"
       include 'nexternal.inc'
       include "coupl.inc"
       include 'q_es.inc'
       include "run.inc"
-      include 'reweight.inc'
 
       double precision p(0:3,nexternal),collrem_xi,collrem_lxi
       double precision xi_i_fks,y_ij_fks
@@ -3679,8 +3760,6 @@
       include "run.inc"
       include "nexternal.inc"
       integer j_fks
-      double precision dlum
-      external dlum
       double precision zhw_used,xi_i_fks,xlum_mc_fact
       double precision xbjrk_ev(2),xbjrk_cnt(2,-2:2)
       common/cbjorkenx/xbjrk_ev,xbjrk_cnt
@@ -3808,15 +3887,6 @@
       endif
       if(.not.pass)i_momcmp_count=i_momcmp_count +1
 c
-      if(jac_cnt(0).gt.0.d0.and.jac.gt.0.d0)
-     #  call p_ev_vs_cnt(izero,i_fks,j_fks,p,p1_cnt,
-     #                   p_i_fks_ev,p_i_fks_cnt,
-     #                   xi_i_fks_ev,y_ij_fks_ev)
-      if(jac_cnt(1).gt.0.d0.and.jac.gt.0.d0)
-     #  call p_ev_vs_cnt(ione,i_fks,j_fks,p,p1_cnt,
-     #                   p_i_fks_ev,p_i_fks_cnt,
-     #                   xi_i_fks_ev,y_ij_fks_ev)
-c
       return
       end
 
@@ -3915,72 +3985,6 @@
       end
 
 
-      subroutine xmcompare_fsr(verbose,inum,iden,i_fks,j_fks,p,p1_cnt)
-      implicit none
-      include 'genps.inc'
-      include 'nexternal.inc'
-      logical verbose
-      integer inum,iden,i_fks,j_fks,iunit,ipart,i
-      double precision tiny,xnum,xden,xrat
-      double precision p(0:3,-max_branch:max_particles)
-      double precision p1_cnt(0:3,nexternal,-2:2)
-      parameter (iunit=6)
-      parameter (tiny=1.d-4)
-c
-      do ipart=1,nexternal
-        do i=0,3
-          xnum=p1_cnt(i,ipart,inum)
-          xden=p1_cnt(i,ipart,iden)
-          if(verbose)then
-            if(i.eq.0)then
-              write(iunit,*)' '
-              write(iunit,*)'part=',ipart
-            endif
-            call xprintout(iunit,xnum,xden)
-          else
-            if(ipart.ne.i_fks.and.ipart.ne.j_fks)then
-              if(xden.ne.0.d0)then
-                xrat=abs(1-xnum/xden)
-              else
-                xrat=abs(xnum)
-              endif
-              if(xrat.gt.tiny)then
-                write(*,*)'Kinematics of counterevents'
-                write(*,*)inum,iden
-                write(*,*)'is different. Particle:',ipart
-                stop
-              endif
-            endif
-          endif
-        enddo
-      enddo
-      do i=0,3
-        xnum=p1_cnt(i,i_fks,inum)+p1_cnt(i,j_fks,inum)
-        xden=p1_cnt(i,i_fks,iden)+p1_cnt(i,j_fks,iden)
-        if(verbose)then
-          if(i.eq.0)then
-            write(iunit,*)' '
-            write(iunit,*)'part=i+j'
-          endif
-          call xprintout(iunit,xnum,xden)
-        else
-          if(xden.ne.0.d0)then
-            xrat=abs(1-xnum/xden)
-          else
-            xrat=abs(xnum)
-          endif
-          if(xrat.gt.tiny)then
-            write(*,*)'Kinematics of counterevents'
-            write(*,*)inum,iden
-            write(*,*)'is different. Particle i+j'
-            stop
-          endif
-        endif
-      enddo
-      return
-      end
-
-
       subroutine xprintout(iunit,xv,xlim)
       implicit real*8(a-h,o-z)
 c
@@ -3993,70 +3997,6 @@
       end
 
 
-      subroutine p_ev_vs_cnt(icnt,i_fks,j_fks,p,p1_cnt,
-     #                       p_i_fks_ev,p_i_fks_cnt,
-     #                       xi_i_fks_ev,y_ij_fks_ev)
-      implicit none
-      include 'genps.inc'
-      include 'nexternal.inc'
-      integer icnt,i_fks,j_fks,ipart,i
-      double precision p(0:3,-max_branch:max_particles)
-      double precision p1_cnt(0:3,nexternal,-2:2)
-      double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
-      double precision xi_i_fks_ev,y_ij_fks_ev,tiny
-      double precision rat(0:3,nexternal+3),den(0:3,nexternal+3)
-      integer maxrat
-c
-c This routine is obsolete; the convergence checks are done elsewhere
-      return
-
-      do ipart=1,nexternal
-        do i=0,3
-          den(i,ipart)=p1_cnt(i,ipart,icnt)
-          if(den(i,ipart).ne.0.d0)then
-            rat(i,ipart)=p(i,ipart)/den(i,ipart)
-          else
-            rat(i,ipart)=p(i,ipart)
-          endif
-        enddo
-      enddo
-c
-      do i=0,3
-        den(i,nexternal+1)=p1_cnt(i,i_fks,icnt)+p1_cnt(i,j_fks,icnt)
-        if(den(i,nexternal+1).ne.0.d0)then
-          rat(i,nexternal+1)=(p(i,i_fks)+p(i,j_fks))/den(i,nexternal+1)
-        else
-          rat(i,nexternal+1)=p(i,i_fks)+p(i,j_fks)
-        endif
-      enddo
-c
-      if(icnt.eq.0)then
-        tiny=4*xi_i_fks_ev
-        maxrat=nexternal+3
-        do i=0,3
-          den(i,nexternal+2)=p_i_fks_cnt(i,0)
-          if(den(i,nexternal+2).ne.0.d0)then
-            rat(i,nexternal+2)=p_i_fks_ev(i)/den(i,nexternal+2)
-          else
-            rat(i,nexternal+2)=p_i_fks_ev(i)
-          endif
-        enddo
-        do i=0,3
-          den(i,nexternal+3)=p_i_fks_cnt(i,0)
-          if(den(i,nexternal+3).ne.0.d0)then
-            rat(i,nexternal+3)=p(i,i_fks)/den(i,nexternal+3)
-          else
-            rat(i,nexternal+3)=p(i,i_fks)
-          endif
-        enddo
-      else
-        tiny=2*sqrt(1-y_ij_fks_ev)
-        maxrat=nexternal+1
-      endif
-c
-      return
-      end
-
 
 c The following has been derived with minor modifications from the
 c analogous routine written for VBF
@@ -4213,183 +4153,8 @@
 
 
 
-
-      subroutine checksij(xsijvc,xsijlvc,xsijlim,
-     #                    xsumvc,xsumlvc,xsumlim,
-     #                    check,checkl,tolerance,
-     #                    iflag,imax,iev,ki,kk,ll,
-     #                    i_fks,j_fks,ilim,iret)
-c Analogous to checkres. Relevant to S functions
-      implicit none
-      real*8 xsijvc(15),xsijlvc,xsumvc(15),xsumlvc,check(15),checkl
-      real*8 xsijlim,xsumlim,tolerance
-      real*8 xsecvc(15),xseclvc
-      real*8 ckc(15),rckc(15),rat
-      logical found
-      integer iflag,imax,iev,ki,kk,ll,i_fks,j_fks,ilim,iret,ithrs,
-     # istop,iwrite,i,imin,icount,itype
-      parameter (ithrs=3)
-      parameter (istop=0)
-      parameter (iwrite=1)
-c
-      if(imax.gt.15)then
-        write(6,*)'Error in checksij: imax is too large',imax
-        stop
-      endif
-      itype=1
-      iret=0
- 100  continue
-      if(itype.eq.1)then
-        do i=1,imax
-          xsecvc(i)=xsijvc(i)
-        enddo
-        xseclvc=xsijlvc
-      elseif(itype.eq.2)then
-        do i=1,imax
-          xsecvc(i)=xsumvc(i)
-        enddo
-        xseclvc=xsumlvc
-      else
-        write(6,*)'Error in checksij: itype=',itype
-        stop
-      endif
-      do i=1,imax
-        if(xseclvc.eq.0.d0)then
-          ckc(i)=abs(xsecvc(i))
-        else
-          ckc(i)=abs(xsecvc(i)/xseclvc-1.d0)
-        endif
-      enddo
-      if(iflag.eq.0)then
-        rat=8.d0
-      elseif(iflag.eq.1)then
-        rat=2.d0
-      else
-        write(6,*)'Error in checksij: iflag=',iflag
-        write(6,*)' Must be 0 for soft, 1 for collinear'
-        stop
-      endif
-c
-      i=1
-      dowhile(ckc(i).gt.0.1d0)
-        i=i+1
-      enddo
-      imin=i
-      do i=imin,imax-1
-        if(ckc(i+1).gt.1.d-8)then
-c If this condition is replaced by .eq.0, the test will fail if the series
-c is made of elements all equal to the limit
-          rckc(i)=ckc(i)/ckc(i+1)
-        else
-c Element #i+1 of series equal to the limit, so it must pass the test
-          rckc(i)=rat*1.1d0
-        endif
-      enddo
-      icount=0
-      i=imin
-      dowhile(icount.lt.ithrs.and.i.lt.imax)
-        if(rckc(i).gt.rat)then
-          icount=icount+1
-        else
-          icount=0
-        endif
-        i=i+1
-      enddo
-c
-      if(icount.ne.ithrs)then
-        iret=iret+itype
-        if(istop.eq.1)then
-          write(6,*)'Test failed',iflag
-          write(6,*)'Event #',iev
-          stop
-        endif
-      endif
-      if(itype.eq.1.and.ki.eq.1.and.iflag.eq.0)then
-        itype=2
-        goto 100
-      endif
-c
-      if(ki.eq.1.and.ilim.eq.1)then
-        found=.false.
-        i=0
-        do while ((.not.found).and.i.lt.imax)
-          i=i+1
-          if(abs(check(i)-1.d0).gt.tolerance)then
-            found=.true.
-            itype=4
-          endif
-        enddo
-        if(.not.found)then
-          if(abs(checkl-1.d0).gt.tolerance)itype=4
-        endif
-        if(itype.eq.4)iret=iret+itype
-      endif
-c
-      if( iwrite.eq.1 .and.
-     #    iret.eq.1 .or.(iret.gt.1.and.ki.eq.1) )then
-        if(iret.gt.7)then
-          write(6,*)'Error in checksij: iret=',iret
-          stop
-        endif
-        write(77,*)'    '
-        if(iflag.eq.0)then
-          write(77,*)'Soft #',iev
-        elseif(iflag.eq.1)then
-          write(77,*)'Collinear #',iev
-        endif
-        write(77,*)'iret:',iret
-        write(77,*)'i_fks,j_fks:',i_fks,j_fks
-        if(iret.eq.1.or.iret.eq.3.or.iret.eq.5.or.iret.eq.7)then
-          write(77,*)'S_kl'
-          write(77,*)'k,kk,ll',ki,kk,ll
-          do i=1,imax
-             call xprintout(77,xsijvc(i),xsijlvc)
-          enddo
-        endif
-        if(iret.eq.2.or.iret.eq.3.or.iret.eq.6.or.iret.eq.7)then
-          write(77,*)'sum of S'
-          do i=1,imax
-             call xprintout(77,xsumvc(i),xsumlvc)
-          enddo
-        endif
-        if(iret.eq.4.or.iret.eq.5.or.iret.eq.6.or.iret.eq.7)then
-          write(77,*)'check to one'
-          do i=1,imax
-             call xprintout(77,check(i),checkl)
-          enddo
-        endif
-      endif
-c
-      if(ilim.eq.1)then
-        if( abs(xsijlvc-xsijlim).gt.1.d-6 .and. 
-     #    xsijlim.ne.-1.d0 )iret=iret+10
-        if( abs(xsumlvc-xsumlim).gt.1.d-6 .and.
-     #    xsumlim.ne.-1.d0 .and. iflag.eq.0)iret=iret+20
-        if(iwrite.eq.1.and.iret.ge.10)then
-          write(77,*)'    '
-          if(iflag.eq.0)then
-            write(77,*)'Soft #',iev
-          elseif(iflag.eq.1)then
-            write(77,*)'Collinear #',iev
-          endif
-          write(77,*)'iret:',iret
-          write(77,*)'i_fks,j_fks:',i_fks,j_fks
-          if((iret.ge.10.and.iret.lt.20).or.iret.ge.30)then
-            write(77,*)'limit of S_kl'
-            write(77,*)'k,kk,ll',ki,kk,ll
-            write(77,*)xsijlvc,xsijlim
-          endif
-          if(iret.ge.20)then
-            write(77,*)'limit of sum_j S_ij'
-            write(77,*)xsumlvc,xsumlim
-          endif
-        endif
-      endif
-      return
-      end
-
-
       subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt)
+      use extra_weights
       implicit none
       include "genps.inc"
       include 'nexternal.inc'
@@ -4401,7 +4166,6 @@
       common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
       include "run.inc"
       include "fks_powers.inc"
-      include 'reweight.inc'
       include "mint.inc"
       double precision p(0:3,nexternal),bsv_wgt,born_wgt,avv_wgt
       double precision pp(0:3,nexternal)
@@ -4755,10 +4519,6 @@
       double precision bpower,born_wgt
       double complex wgt1(2)
 
-      integer           isum_hel
-      logical                   multi_channel
-      common/to_matrix/isum_hel, multi_channel
-      integer isum_hel_orig
       integer i_fks,j_fks
       common/fks_indices/i_fks,j_fks
 
@@ -4767,13 +4527,6 @@
 
       double precision tiny
       parameter (tiny=1d-6)
-
-c Make sure that we sum over helicities (such that we do get a
-c non-zero Born)
-      isum_hel_orig = isum_hel
-      isum_hel=0
-      call get_helicity(i_fks,j_fks)
-
       calculatedBorn=.false.
       call sborn(p_born,wgt1)
 c Born contribution:
@@ -4815,7 +4568,6 @@
 c nothing funny happens later on
       g=g/10d0
       call update_as_param()
-      isum_hel=isum_hel_orig
       calculatedBorn=.false.
       call sborn(p_born,wgt1)
 
@@ -4825,25 +4577,21 @@
 c       This function computes the power of a muR-dependent factor which
 c       is stored in cpower. You need to modify it when you try to 
 c       reweight your cross section with a muR-dependent factor
-c       (runfac=1 in reweight0.inc)
+c       (runfac=1 in extra_weights.f)
 c Note: The implementation below only works for the Bottom Yukawa in
 c       the SM where "GC_33" contains the Yukawa, for other models
 c       or general muR-dependent factors you need to change GC_33
 c       to the corresponding coupling.
       subroutine compute_cpower(p_born,cpower)
+      use extra_weights
       implicit none
       include "nexternal.inc"
       include "coupl.inc"
-      include 'reweight.inc'
 
       double precision p_born(0:3,nexternal-1)
       double precision cpower,born_wgt
       double complex wgt1(2)
 
-      integer isum_hel
-      logical multi_channel
-      common/to_matrix/isum_hel, multi_channel
-      integer isum_hel_orig
       integer i_fks,j_fks
       common/fks_indices/i_fks,j_fks
 
@@ -4859,12 +4607,6 @@
 
 c   The following is relevant for a muR-dependent bottom-mass in Yukawa.
 c$$$
-c$$$c Make sure that we sum over helicities (such that we do get a
-c$$$c non-zero Born)
-c$$$      isum_hel_orig = isum_hel
-c$$$      isum_hel=0
-c$$$      call get_helicity(i_fks,j_fks)
-c$$$
 c$$$      calculatedBorn=.false.
 c$$$      call sborn(p_born,wgt1)
 c$$$c Born contribution:
@@ -4900,7 +4642,7 @@
 c$$$         cpower=dble(nint(cpower))
 c$$$         write(*,*)'cpower is', cpower
 c$$$c Check consistency with value used in reweighting
-c$$$c$$$         if( (doreweight.or.doNLOreweight) .and.
+c$$$c$$$         if( doreweight .and.
 c$$$c$$$     &        abs(cpower-wgtcpower).gt.tiny )then
 c$$$c$$$            write(*,*)'Error in compute_cpower'
 c$$$c$$$            write(*,*)'cpower(s) are:',cpower,wgtcpower
@@ -4911,7 +4653,6 @@
 c$$$c Change couplings back and recompute the Born to make sure that 
 c$$$c nothing funny happens later on
 c$$$      GC_33 = GC_33 / 10d0
-c$$$      isum_hel=isum_hel_orig
 c$$$      calculatedBorn=.false.
 c$$$      call sborn(p_born,wgt1)
 c$$$
@@ -5348,73 +5089,9 @@
       end
 
 
-      function m1l_finite_CDR(p,born)
-c Returns the finite part of virtual contribution, according to the
-c definitions given in (B.1) and (B.2). This function must include
-c the factor as/(2*pi)
-      implicit none
-      include "genps.inc"
-      include 'nexternal.inc'
-c      include "fks.inc"
-      integer fks_j_from_i(nexternal,0:nexternal)
-     &     ,particle_type(nexternal),pdg_type(nexternal)
-      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
-      include 'coupl.inc'
-      include 'q_es.inc'
-      double precision p(0:3,nexternal-1),m1l_finite_CDR,born
-      double precision CF,pi,aso2pi,shat,dot,xlgq2os
-      parameter (CF=4d0/3d0)
-      parameter (pi=3.1415926535897932385d0)
-c
-      aso2pi=g**2/(8*pi**2)
-c This is relevant to e+e- --> qqbar
-      shat=2d0*dot(p(0,1),p(0,2))
-      xlgq2os=log(QES2/shat)
-      m1l_finite_CDR=-aso2pi*CF*(xlgq2os**2+3*xlgq2os-pi**2+8.d0)*born
-      return
-      end
-
-
-      function m1l_W_finite_CDR(p,born)
-c Returns the finite part of virtual contribution, according to the
-c definitions given in (B.1) and (B.2). This function must include
-c the factor as/(2*pi)
-      implicit none
-      include "genps.inc"
-      include 'nexternal.inc'
-c      include "fks.inc"
-      integer fks_j_from_i(nexternal,0:nexternal)
-     &     ,particle_type(nexternal),pdg_type(nexternal)
-      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
-      include 'coupl.inc'
-      include 'q_es.inc'
-      double precision p(0:3,nexternal-1),m1l_W_finite_CDR,born
-      double precision CF,pi,aso2pi,shat,dot,xlgq2os
-      parameter (CF=4d0/3d0)
-      parameter (pi=3.1415926535897932385d0)
-c
-      aso2pi=g**2/(8*pi**2)
-      shat=2d0*dot(p(0,1),p(0,2))
-      xlgq2os=log(QES2/shat)
-
-c This is relevant to qqbar -> W 
-      m1l_W_finite_CDR=aso2pi*CF*(-xlgq2os**2-3d0*xlgq2os+pi**2-8d0)
-      m1l_W_finite_CDR=m1l_W_finite_CDR*born
-
-c This is relevant to gg -> H
-c$$$      m1l_W_finite_CDR=aso2pi*(-3d0*xlgq2os**2+11d0+3d0*pi**2)
-c$$$      m1l_W_finite_CDR=m1l_W_finite_CDR*born
-
-c This is relevant to bbbar -> H
-c$$$      m1l_W_finite_CDR=aso2pi
-c$$$     f     * (-4d0/3d0*xlgq2os**2
-c$$$     f        -8d0/3d0+(16d0/3d0+8d0/3d0)*pi**2/6d0)
-c$$$      m1l_W_finite_CDR=m1l_W_finite_CDR*born
-      return
-      end
-
-
       subroutine setfksfactor(match_to_shower)
+      use weight_lines
+      use extra_weights
       implicit none
 
       include 'mint.inc'
@@ -5449,11 +5126,9 @@
       include 'nexternal.inc'
       include 'fks_powers.inc'
       include 'nFKSconfigs.inc'
-      include 'c_weight.inc'
       integer fks_j_from_i(nexternal,0:nexternal)
      &     ,particle_type(nexternal),pdg_type(nexternal)
       common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
-      include 'reweight0.inc'
       include 'run.inc'
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
@@ -5473,8 +5148,6 @@
       common /cdelta_used/delta_used
       double precision xiScut_used,xiBSVcut_used
       common /cxiScut_used/xiScut_used,xiBSVcut_used
-      logical rotategranny
-      common/crotategranny/rotategranny
       double precision diagramsymmetryfactor_save(maxchannels)
       save diagramsymmetryfactor_save
       double precision diagramsymmetryfactor
@@ -5511,10 +5184,6 @@
       integer i_type,j_type,m_type
       common/cparticle_types/i_type,j_type,m_type
 
-c The value of rotategranny may be superseded later if phase space
-c parametrization allows it
-      rotategranny=.false.
-
       softtest=.false.
       colltest=.false.
       fold=0
@@ -5691,6 +5360,8 @@
 c Compute the identical particle symmetry factor that is in the
 c Born matrix elements.
          iden_born_FKS(nFKSprocess)=1
+         call weight_lines_allocated(nexternal,max_contr,max_wgt
+     $        ,max_iproc)
          call set_pdg(0,nFKSprocess)
          do i=1,nexternal
             iden(i)=1
@@ -5800,349 +5471,12 @@
       end
 
 
-      subroutine get_helicity(i_fks,j_fks)
-      implicit none
-      include "nexternal.inc"
-      include "born_nhel.inc"
-      include "madfks_mcatnlo.inc"
-      integer NHEL(nexternal,max_bhel*2),IHEL
-chel  include "helicities.inc"
-      include 'nFKSconfigs.inc'
-      double precision hel_fac
-      integer get_hel,skip(fks_configs)
-      common/cBorn/hel_fac,get_hel,skip
-      logical calculatedBorn
-      common/ccalculatedBorn/calculatedBorn
-      integer hel_wgt,hel_wgt_born,hel_wgt_real
-      integer nhelreal(nexternal,4),goodhelreal(4)
-      integer nhelrealall(nexternal,max_bhel*2)
-      common /c_nhelreal/ nhelreal,nhelrealall,goodhelreal,hel_wgt_real
-      integer nhelborn(nexternal-1,2),goodhelborn(2)
-      integer nhelbornall(nexternal-1,max_bhel)
-      common /c_nhelborn/ nhelborn,nhelbornall,goodhelborn,hel_wgt_born
-
-      integer           isum_hel
-      logical                   multi_channel
-      common/to_matrix/isum_hel, multi_channel
-
-      integer i,nexthel,j,i_fks,j_fks,ngood,k
-      data nexthel /0/
-      data ngood /0/
-      logical done,firsttime,all_set,chckr
-      data firsttime/.true./
-      integer goodhelr(0:4,max_bhel/2),goodhelb(0:2,max_bhel/2)
-      save goodhelr,goodhelb,all_set,chckr
-      double precision rnd,ran2
-      external ran2
-
-      character*4 abrv
-      common /to_abrv/ abrv
-      logical Hevents
-      common/SHevents/Hevents
-      logical usexinteg,mint
-      common/cusexinteg/usexinteg,mint
-
-c Do not change these two lines, because ./bin/compile_madfks.sh might
-c need to change them automatically
-      logical HelSum
-      parameter (HelSum=.true.)
-
-c************
-c goodhelr=2, real emission matrix element not yet calculated
-c             for this helicity
-c goodhelr=1, real emission matrix element calculated and non-zero
-c goodhelr=0, real emission matrix element calculated and zero,
-c             so can be skipped next time.
-c************
-      if (HelSum) return
-
-      if (isum_hel.ne.0) then ! MC over helicities
-c First, set the goodhelr and goodhelb to their starting values
-      if (firsttime) then
-         if ((mint .and. (.not.Hevents) .and. (abrv(1:2).eq.'vi' .or.
-     &        abrv.eq.'born' .or. abrv.eq.'grid' .or.
-     &        (.not.UseSudakov))) .or. (.not.mint .and. (abrv.eq.'born'
-     &        .or. abrv.eq.'grid' .or. abrv(1:2).eq.'vi'))) then
-c           if computing only the Born diagrams, should not
-c           consider real emission helicities            
-            chckr=.false.
-         else
-            chckr=.true.
-         endif
-         do i=1,fks_configs
-            skip(i)=1
-         enddo
-c read from file if possible
-         open(unit=65,file='goodhel.dat',status='old',err=532)
-         all_set=.true.
-         do j=0,4
-            read (65,*,err=532) (goodhelr(j,i),i=1,max_bhel/2)
-         enddo
-         do j=0,2
-            read (65,*,err=532) (goodhelb(j,i),i=1,max_bhel/2)
-         enddo
-         read(65,*,err=532) hel_wgt
-         hel_wgt_born=hel_wgt
-         hel_wgt_real=hel_wgt
-         do i=1,max_bhel/2
-            if ((chckr .and.
-     &           (goodhelb(0,i).eq.2 .or. goodhelr(0,i).eq.2)) .or.
-     &           (.not.chckr.and.goodhelb(0,i).eq.2)) all_set=.false.
-         enddo
-         close(65)
-         goto 533
-c if file does not exist or has wrong format, set all to 2
- 532     close(65)
-         write (*,*) 'Good helicities not found in file'
-         all_set=.false.
-         do j=0,4
-            do i=1,max_bhel/2
-               goodhelr(j,i)=2
-            enddo
-         enddo
-         do j=0,2
-            do i=1,max_bhel/2
-               goodhelb(j,i)=2
-            enddo
-         enddo
-         hel_wgt=max_bhel/2
-         hel_wgt_born=hel_wgt
-         hel_wgt_real=hel_wgt
- 533     continue
-         firsttime=.false.
-         goto 534 ! no previous event, so skip to the next helicity
-      endif
-
-c From previous event, check if there is an update
-      if (.not.all_set) then
-c real emission
-         if(goodhelr(0,ngood).eq.2) then
-            if ( goodhelreal(1).eq.0 .and.
-     &           goodhelreal(2).eq.0 .and.
-     &           goodhelreal(3).eq.0 .and.
-     &           goodhelreal(4).eq.0 ) then
-               do j=0,4
-                  goodhelr(j,ngood)=0
-               enddo
-            elseif( goodhelreal(1).le.1 .and.
-     &              goodhelreal(2).le.1 .and.
-     &              goodhelreal(3).le.1 .and.
-     &              goodhelreal(4).le.1 ) then
-               goodhelr(0,ngood)=1
-               do j=1,4
-                  goodhelr(j,ngood)=goodhelreal(j)
-               enddo
-            elseif (.not.(goodhelreal(1).eq.2 .and.
-     &                    goodhelreal(2).eq.2 .and.
-     &                    goodhelreal(2).eq.2 .and.
-     &                    goodhelreal(2).eq.2) ) then
-               write (*,*) 'Error #2 in get_helicities',
-     &              ngood,(goodhelr(j,ngood),j=0,4)
-               stop
-            endif
-         endif
-c Born and counter events
-         if(goodhelb(0,ngood).eq.2) then
-            if ( goodhelborn(1).eq.0 .and.
-     &           goodhelborn(2).eq.0 ) then
-               do j=0,2
-                  goodhelb(j,ngood)=0
-               enddo
-            elseif( goodhelborn(1).le.1 .and.
-     &              goodhelborn(2).le.1 ) then
-               goodhelb(0,ngood)=1
-               do j=1,2
-                  goodhelb(j,ngood)=goodhelborn(j)
-               enddo
-            elseif (.not.(goodhelborn(1).eq.2 .and.
-     &                    goodhelborn(2).eq.2) ) then
-               write (*,*) 'Error #3 in get_helicities',
-     &              nexthel,(goodhelb(j,ngood),j=0,2)
-               stop
-            endif
-         endif
-
-c Calculate new hel_wgt
-         hel_wgt=0
-         do i=1,max_bhel/2
-            if((chckr .and.
-     &           (goodhelb(0,i).ge.1.or.goodhelr(0,i).ge.1)) .or.
-     &           (.not.chckr .and. goodhelb(0,i).ge.1)) then
-               hel_wgt=hel_wgt+1
-            endif
-         enddo
-         hel_wgt_born=hel_wgt
-         hel_wgt_real=hel_wgt
-
-c check if all have been set, if so -> write to file
-         all_set=.true.
-         do i=1,max_bhel/2
-            if ((chckr .and.
-     &           (goodhelb(0,i).eq.2 .or. goodhelr(0,i).eq.2)) .or.
-     &           (.not.chckr.and.goodhelb(0,i).eq.2)) all_set=.false.
-         enddo
-         if (all_set) then
-            write (*,*) 'All good helicities have been found.',hel_wgt
-            open(unit=65,file='goodhel.dat',status='unknown')
-            do j=0,4
-               write (65,*) (goodhelr(j,i),i=1,max_bhel/2)
-            enddo
-            do j=0,2
-               write (65,*) (goodhelb(j,i),i=1,max_bhel/2)
-            enddo
-            write(65,*) hel_wgt
-            close(65)
-         endif
-      else
-         do i=1,4
-            if (goodhelr(i,ngood).ne.goodhelreal(i)) then
-               write (*,*)'Error #4 in get_helicities',i,ngood
-               stop
-            endif
-         enddo
-         do i=1,2
-            if (goodhelb(i,ngood).ne.goodhelborn(i)) then
-               write (*,*)'Error #5 in get_helicities',i,ngood
-               stop
-            endif
-         enddo
-      endif
-
-c Get the next helicity
- 534  continue
-      done=.false.
-      do while (.not.done)
-         if (nexthel.eq.max_bhel*2) nexthel=0
-         nexthel=nexthel+1
-         if(nhel(i_fks,nexthel).eq.1.and.nhel(j_fks,nexthel).eq.1) then
-            if (ngood.eq.max_bhel/2) ngood=0
-            ngood=ngood+1
-            if((chckr .and.
-     &           (goodhelr(0,ngood).ge.1.or.goodhelb(0,ngood).ge.1)).or.
-     &           (.not.chckr .and. goodhelb(0,ngood).ge.1)) then
-c Using random number to see if we have to go to the next.
-c Probably this is an overkill, but have to make sure that there is
-c no bias considering the *semi*-random numbers from VEGAS.
-               rnd=ran2()
-               if (rnd.le.1d0/dble(hel_wgt)) then
-                  done=.true.
-               endif
-            endif
-         endif
-      enddo
-
-      do i=1,nexternal
-         if (i.eq.i_fks) then
-            nhelreal(i,1)=1
-            nhelreal(i,2)=1
-            nhelreal(i,3)=-1
-            nhelreal(i,4)=-1
-         elseif (i.eq.j_fks) then
-            nhelreal(i,1)=1
-            nhelreal(i,2)=-1
-            nhelreal(i,3)=1
-            nhelreal(i,4)=-1
-         else
-            nhelreal(i,1)=nhel(i,nexthel)
-            nhelreal(i,2)=nhel(i,nexthel)
-            nhelreal(i,3)=nhel(i,nexthel)
-            nhelreal(i,4)=nhel(i,nexthel)
-         endif
-      enddo
-      do j=1,4
-         goodhelreal(j)=goodhelr(j,ngood)
-      enddo
-
-      do i=1,nexternal-1
-         if (i.eq.min(i_fks,j_fks)) then
-            nhelborn(i,1)=1
-            nhelborn(i,2)=-1
-         elseif(i.lt.max(i_fks,j_fks)) then
-            nhelborn(i,1)=nhel(i,nexthel)
-            nhelborn(i,2)=nhel(i,nexthel)
-         else
-            nhelborn(i,1)=nhel(i+1,nexthel)
-            nhelborn(i,2)=nhel(i+1,nexthel)
-         endif
-      enddo
-      do j=1,2
-         goodhelborn(j)=goodhelb(j,ngood)
-      enddo
-
-      else !isum_hel is zero, sum explicitly over helicities
-
-      do i=1,nexternal
-         do j=1,max_bhel*2
-            nhelrealall(i,j)=nhel(i,j)
-         enddo
-      enddo
-      do i=1,nexternal-1
-         k=0
-         do j=1,max_bhel*2
-            if (nhel(i_fks,j).eq.-1) then
-               k=k+1
-               if (i.lt.i_fks) then
-                  nhelbornall(i,k)=nhel(i,j)                  
-               elseif(i.ge.i_fks) then
-                  nhelbornall(i,k)=nhel(i+1,j)
-               endif
-            endif
-         enddo
-      enddo
-
-      endif
-      return
-      end
-
-      function get_ptrel(pp,i_fks,j_fks)
-      implicit none
-      include 'nexternal.inc'
-      double precision get_ptrel,pp(0:3,nexternal)
-      integer i_fks,j_fks
-      double precision tmp,psum(3)
-      integer i
-c
-      if(j_fks.le.2)then
-        tmp=sqrt(pp(1,i_fks)**2+pp(2,i_fks)**2)
-      else
-        do i=1,3
-          psum(i)=pp(i,i_fks)+pp(i,j_fks)
-        enddo
-        tmp=( pp(2,i_fks)*psum(1)-pp(1,i_fks)*psum(2) )**2+
-     #      ( pp(3,i_fks)*psum(1)-pp(1,i_fks)*psum(3) )**2+
-     #      ( pp(3,i_fks)*psum(2)-pp(2,i_fks)*psum(3) )**2
-        if(tmp.ne.0.d0)tmp=sqrt( tmp/
-     #       (psum(1)**2+psum(2)**2+psum(3)**2) )
-      endif
-      get_ptrel=tmp
-      return
-      end
-
-
-
-      FUNCTION FK88RANDOM(SEED)
-*     -----------------
-* Ref.: K. Park and K.W. Miller, Comm. of the ACM 31 (1988) p.1192
-* Use seed = 1 as first value.
-*
-      IMPLICIT INTEGER(A-Z)
-      REAL*8 MINV,FK88RANDOM
-      SAVE
-      PARAMETER(M=2147483647,A=16807,Q=127773,R=2836)
-      PARAMETER(MINV=0.46566128752458d-09)
-      HI = SEED/Q
-      LO = MOD(SEED,Q)
-      SEED = A*LO - R*HI
-      IF(SEED.LE.0) SEED = SEED + M
-      FK88RANDOM = SEED*MINV
-      END
-
 
       subroutine set_mu_central(ic,dd,c_mu2_r,c_mu2_f)
+      use weight_lines
+      use extra_weights
       implicit none
       include 'nexternal.inc'
-      include 'c_weight.inc'
-      include 'reweight0.inc'
       include 'run.inc'
       integer ic,dd,i,j
       double precision c_mu2_r,c_mu2_f,muR,muF,pp(0:3,nexternal)
 
=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
--- Template/NLO/SubProcesses/genps_fks.f	2017-05-16 14:24:36 +0000
+++ Template/NLO/SubProcesses/genps_fks.f	2017-08-16 21:26:30 +0000
@@ -2,8 +2,6 @@
       implicit none
       include 'genps.inc'
       include 'nexternal.inc'
-c     Timing profile statistics
-      include 'timing_variables.inc'
       integer ndim,iconfig
       double precision wgt,x(99),p(0:3,nexternal)
       integer iforest(2,-max_branch:-1,lmaxconfigs)
@@ -32,8 +30,6 @@
       include 'coupl.inc'
       include 'born_props.inc'
 c     
-      call cpu_time(tBefore)
-c
       this_config=iconfig
       iconf=iconfig
       iconfig0=iconfig
@@ -56,8 +52,6 @@
       enddo
       wgt=wgt*jac
 c
-      call cpu_time(tAfter)      
-      tGenPS = tGenPS + (tAfter-tBefore)
       return
       end 
       
@@ -313,8 +307,7 @@
 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_level,cBW_mass,cBW_width,s,x
-     $     ,xjac0,pass)
+     $     ,totmass,qwidth,qmass,cBW,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
@@ -1215,7 +1208,7 @@
       integer icountevts,i_fks,j_fks
       double precision xbjrk_born(2),tau_born,ycm_born,ycmhat,shat_born
      &     ,phi_i_fks,xpswgt,xjac,xiimax,xinorm,xp(0:3,nexternal),stot
-     &     ,x(2) ,y_ij_fks
+     &     ,x(2),y_ij_fks
       double precision shat,sqrtshat,tau,ycm,xbjrk(2),p_i_fks(0:3)
       logical pass
 c common blocks
@@ -1575,40 +1568,6 @@
       return
       end
 
-
-      function bwfunc(s,xm02,gah)
-c Returns the Breit Wigner function, normalized in such a way that
-c its integral in the range (-inf,inf) is one
-      implicit none
-      real*8 bwfunc,s,xm02,gah
-      real*8 pi,xm0
-      parameter (pi=3.1415926535897932d0)
-c
-      xm0=sqrt(xm02)
-      bwfunc=xm0*gah/(pi*((s-xm02)**2+xm02*gah**2))
-      return
-      end
-
-
-      function xbwmass3(t,xm02,ga,bwdelf,bwfmmn)
-c Returns the boson mass squared, given 0<t<1, the nominal mass (xm0),
-c and the mass range (implicit in bwdelf and bwfmmn). This function
-c is the inverse of F(M^2), where
-c   F(M^2)=\int_{xmlow2}^{M^2} ds BW(sqrt(s),M0,Ga)
-c   BW(M,M0,Ga)=M0 Ga/pi 1/((M^2-M0^2)^2+M0^2 Ga^2
-c and therefore eats up the Breit-Wigner when changing integration 
-c variable M^2 --> t
-      implicit none
-      real*8 xbwmass3,t,xm02,ga,bwdelf,bwfmmn
-      real*8 pi,xm0
-      parameter (pi=3.1415926535897932d0)
-c
-      xm0=sqrt(xm02)
-      xbwmass3=xm02+xm0*ga*tan(pi*bwdelf*t-bwfmmn)
-      return
-      end
-
-
       subroutine gentcms(pa,pb,t,phi,m1,m2,p1,pr,jac)
 c*************************************************************************
 c     Generates 4 momentum for particle 1, and remainder pr
@@ -1774,21 +1733,10 @@
       subroutine generate_tau_BW(stot,idim,x,mass,width,cBW,BWmass
      $     ,BWwidth,tau,jac)
       implicit none
-      real*8 pi
-      parameter (pi=3.1415926535897932d0)
-      integer nsamp
-      parameter (nsamp=1)
-      include 'run.inc'
-      include 'genps.inc'
-      integer cBW,icount,idim
-      data icount /0/
-      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
-      double precision tau_Born_lower_bound_save
-     $     ,tau_lower_bound_resonance_save,tau_lower_bound_save
+      integer cBW,idim
+      double precision stot,x,tau,jac,mass,width,BWmass(-1:1),BWwidth(
+     $     -1:1),s_mass,s
+      double precision smax,smin
       double precision tau_Born_lower_bound,tau_lower_bound_resonance
      &     ,tau_lower_bound
       common/ctau_lower_bound/tau_Born_lower_bound
@@ -1815,18 +1763,15 @@
 
 
       subroutine generate_tau(stot,idim,x,tau,jac)
-      double precision x,tau,jac
-      double precision roH,roHs,fract,ximax0,ximin0,tmp,fract1,fract2
-     &     ,roHj,stot,s,dum,dum3(-1:1),smin,smax,s_mass
-      integer nsamp,idim
-      parameter (nsamp=1)
+      implicit none
+      integer idim
+      double precision x,tau,jac,smin,smax,s_mass,s,tiny,dum,dum3(-1:1)
+     $     ,stot
+      parameter (tiny=1d-8)
       double precision tau_Born_lower_bound,tau_lower_bound_resonance
-     &     ,tau_lower_bound,tiny
-      parameter (tiny=1d-8)
+     $     ,tau_lower_bound
       common/ctau_lower_bound/tau_Born_lower_bound
-     &     ,tau_lower_bound_resonance,tau_lower_bound
-      character*4 abrv
-      common /to_abrv/ abrv
+     $     ,tau_lower_bound_resonance,tau_lower_bound
       smin=tau_born_lower_bound*stot
       smax=stot
       s_mass=tau_lower_bound_resonance*stot
@@ -1887,11 +1832,8 @@
 
       
       subroutine generate_inv_mass_sch(ns_channel,itree,m,sqrtshat_born
-     $     ,totmass,qwidth,qmass,cBW,cBW_level,cBW_mass,cBW_width,s,x
-     $     ,xjac0,pass)
+     $     ,totmass,qwidth,qmass,cBW,cBW_mass,cBW_width,s,x,xjac0,pass)
       implicit none
-      real*8 pi
-      parameter (pi=3.1415926535897932d0)
       include 'genps.inc'
       include 'nexternal.inc'
       integer ns_channel
@@ -1901,18 +1843,12 @@
       double precision sqrtshat_born,totmass,xjac0
       integer itree(2,-max_branch:-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
+      double precision smin,smax,totalmass
       logical pass
-      integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-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:nexternal),xi,fract,vol_new
-     $     ,shat,x_new
-      parameter (fract=0.1d0)
+      integer cBW(-nexternal:-1)
+      double precision cBW_mass(-1:1,-nexternal:-1),cBW_width(-1:1,
+     $     -nexternal:-1)
+      double precision s_mass(-nexternal:nexternal)
       common/to_phase_space_s_channel/s_mass
       pass=.true.
       totalmass=totmass
@@ -2233,9 +2169,8 @@
       include 'nexternal.inc'
       double precision p_born(0:3,nexternal-1),xmrec2,shat_born
       logical pass
-      integer imother
-      integer i
-      double precision recoilbar(0:3),dot,tmp
+      integer imother,i
+      double precision recoilbar(0:3),dot
       external dot
       pass=.true.
       do i=0,3
@@ -2280,15 +2215,10 @@
 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
+      double precision fract,A,B,C,D,E,F,G,bs(-1:1),maxi,mini
       integer j
 c
       if (itype.eq.1) then
 
=== modified file 'Template/NLO/SubProcesses/handling_lhe_events.f'
--- Template/NLO/SubProcesses/handling_lhe_events.f	2016-11-03 15:21:05 +0000
+++ Template/NLO/SubProcesses/handling_lhe_events.f	2017-08-16 21:26:30 +0000
@@ -1,14 +1,4 @@
 c Utility routines for LHEF. Originally taken from collect_events.f
-c
-c Note: the routines read_lhef_event and write_lhef_event use the common
-c blocks in reweight0.inc, relevant to reweight information. This is
-c independent of the process, and in particular of process-related
-c parameters such as nexternal, which is replaced here by (its supposed)
-c upper bound maxparticles. The arrays which have one dimension defined
-c by maxparticles may have a correspondence with process-specific ones,
-c and the dimensions of the latter are typically defined by nexternal.
-c Hence, one may need an explicit copy of one onto the other
-c
 
       block data
       integer event_id
@@ -20,9 +10,9 @@
       end
 
       subroutine write_lhef_header(ifile,nevents,MonteCarlo)
+      use extra_weights
       implicit none 
       include 'run.inc'
-      include 'reweight0.inc'
       integer idwgt,kk,ii,jj,nn,n
       integer ifile,nevents
       character*10 MonteCarlo
@@ -114,6 +104,7 @@
 
 
       subroutine write_lhef_header_banner(ifile,nevents,MonteCarlo,path)
+      use extra_weights
       implicit none 
       integer ifile, i, idwgt, nevents,iseed,ii,jj,kk,nn,n
       double precision mcmass(-16:21)
@@ -135,7 +126,6 @@
       character*150 buffer,buffer_lc,buffer2
       integer event_id
       common /c_event_id/ event_id
-      include 'reweight_all.inc'
 
 c     Set the event_id to 0. If 0 or positive, this value will be update
 c     in write_lhe_event. It is set to -99 through a block data
@@ -297,8 +287,8 @@
 
 
       subroutine read_lhef_header(ifile,nevents,MonteCarlo)
+      use extra_weights
       implicit none 
-      include 'reweight0.inc'
       include './run.inc'
       logical already_found
       integer ifile,nevents,i,ii,ii2,iistr,itemp
@@ -435,8 +425,8 @@
 c Avoid overloading read_lhef_header, meant to be used in utilities
       subroutine read_lhef_header_full(ifile,nevents,itempsc,itempPDF,
      #                                 MonteCarlo)
+      use extra_weights
       implicit none
-      include 'reweight0.inc'
       include 'run.inc'
       logical already_found
       integer ifile,nevents,i,ii,ii2,iistr,ipart,itempsc,itempPDF
@@ -652,6 +642,7 @@
       subroutine write_lhef_event(ifile,
      # NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
      # IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
+      use extra_weights
       implicit none
       INTEGER NUP,IDPRUP,IDUP(*),ISTUP(*),MOTHUP(2,*),ICOLUP(2,*)
       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,
@@ -669,7 +660,6 @@
       common/c_i_process/i_process
       integer nattr,npNLO,npLO
       common/event_attributes/nattr,npNLO,npLO
-      include 'reweight_all.inc'
       include './run.inc'
       include 'unlops.inc'
 c     if event_id is zero or positive (that means that there was a call
@@ -824,6 +814,7 @@
       subroutine read_lhef_event(ifile,
      # NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
      # IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
+      use extra_weights
       implicit none
       INTEGER NUP,IDPRUP,IDUP(*),ISTUP(*),MOTHUP(2,*),ICOLUP(2,*)
       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,
@@ -842,7 +833,6 @@
       common/c_i_process/i_process
       integer nattr,npNLO,npLO
       common/event_attributes/nattr,npNLO,npLO
-      include 'reweight_all.inc'
       include 'unlops.inc'
       include 'run.inc'
 c
@@ -876,11 +866,27 @@
          if(jwgtinfo.eq.-5 .or. jwgtinfo.eq.-9) then
             read(ifile,'(a)')string
             read(ifile,*) wgtref,n_ctr_found,n_mom_conf,wgtcpower
+            if (.not.allocated(momenta_str)) allocate(momenta_str(0:3
+     $           ,max_mext,max_mom_str))
+            if (n_mom_conf.gt.max_mom_str .or. mexternal.gt.max_mext)
+     $           then
+               deallocate(momenta_str)
+               max_mom_str=max(n_mom_conf,max_mom_str)
+               max_mext=max(mexternal,max_mext)
+               allocate(momenta_str(0:3,max_mext,max_mom_str))
+            endif
             do i=1,n_mom_conf
                do j=1,mexternal
                   read (ifile,*) (momenta_str(ii,j,i),ii=0,3)
                enddo
             enddo
+            if (.not.allocated(n_ctr_str))
+     $           allocate(n_ctr_str(max_n_ctr))
+            if (n_ctr_found.gt.max_n_ctr) then
+               deallocate(n_ctr_str)
+               max_n_ctr=n_ctr_found
+               allocate(n_ctr_str(max_n_ctr))
+            endif
             do i=1,n_ctr_found
                read (ifile,'(a)') n_ctr_str(i)
             enddo
@@ -955,6 +961,7 @@
       subroutine read_lhef_event_catch(ifile,
      # NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
      # IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
+      use extra_weights
       implicit none
       INTEGER NUP,IDPRUP,IDUP(*),ISTUP(*),MOTHUP(2,*),ICOLUP(2,*)
       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,
@@ -973,7 +980,6 @@
       common/c_i_process/i_process
       integer nattr,npNLO,npLO
       common/event_attributes/nattr,npNLO,npLO
-      include 'reweight_all.inc'
       include 'unlops.inc'
       include 'run.inc'
 c
@@ -1008,23 +1014,39 @@
       enddo
       read(ifile,'(a)')buff
       if(buff(1:1).eq.'#')then
-        read(buff,*)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
+         read(buff,*)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
      #                    fksfather_lhe,ipartner_lhe,
      #                    scale1_lhe,scale2_lhe,
      #                    jwgtinfo,mexternal,iwgtnumpartn,
      #         wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
-        if(jwgtinfo.eq.-5 .or. jwgtinfo.eq.-9) then
-           read(ifile,'(a)')string
-           read(ifile,*) wgtref,n_ctr_found,n_mom_conf,wgtcpower
-           do i=1,n_mom_conf
-              do j=1,mexternal
-                 read (ifile,*) (momenta_str(ii,j,i),ii=0,3)
-              enddo
-           enddo
-           do i=1,n_ctr_found
-              read (ifile,'(a)') n_ctr_str(i)
-           enddo
-           read(ifile,'(a)')string