Merge lp:~maddevelopers/mg5amcnlo/2.1.1 into lp:mg5amcnlo/lts

[Olivier Mattelaer]

I think that we are now in due time to release this version.
I would propose to do that on next friday (1 month since the latest so this is a good time).

So would be great to move forward on the branches that need to be merge here:
lp:~maddevelopers/mg5amcnlo/2.1.1_IOTests_MadeEasy ⇒ lp:~maddevelopers/mg5amcnlo/2.1.1
--> for me.

lp:~maddevelopers/mg5amcnlo/plugin_model ⇒ lp:~maddevelopers/mg5amcnlo/2.1.1
--> for Valentin

lp:~maddevelopers/mg5amcnlo/2.1.1_IOTests_MadeEasy ⇒ lp:~maddevelopers/mg5amcnlo/2.1.1
--> This is important that you look at this one, since it allows easy tests for the new features.
I tested it already and this is really nice. So please look at it.

Would be great to have all the branches merged soon such that we can make a proper validation of that version since it changes deep stuff (i.e. all the model name have changed)

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/2.1.1]

277. By Olivier Mattelaer on 2014-03-12

fix a problem a fix hardcoded vector size in MadWeight

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,

What about the branch to fix the 'local compilation' option for the experimentalists? Shouldn't we try to include this in the next release? It is not long to implement, and I can try to have it ready for thursday.

Cheers

[Olivier Mattelaer (olivier-mattelaer)]

Sure, this can (should actually) be include.
didn't mention it since it was not yet propose for merging.
But if things need to be done for that release, it can be include (if it is done on time obviously)

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/2.1.1]

278. By Olivier Mattelaer on 2014-03-12

fix a small glitch in Pythia8 output

279. By <email address hidden> on 2014-03-12

remove option BW_effects and zeromass_for_max_weight from MadSpin

280. By Rikkert Frederix on 2014-03-13

1. Fixed a bug for fixed order running introduced in r273.
2. Can now also do separate abrv=virt
3. Removed some obsolete options for 'sum' in driver_mintMC

281. By marco zaro on 2014-03-13

fixed bug for e+e- processes + lhapdf

282. By Rikkert Frederix on 2014-03-13

merge with understanding_cuts branch

283. By Rikkert Frederix on 2014-03-17

merge with NLO_BW_FIX

284. By Olivier Mattelaer on 2014-03-17

modify SysCalc according to the comment of Rikkert

285. By Rikkert Frederix on 2014-03-18

Merge with the running_yukawa branch

286. By Olivier Mattelaer on 2014-03-18

Multiple fixes
1) a problem with pbs cluster
2) a problem with 4l final state and mmll (bug still present, I have just added a warning)
3) add the possibility to put maxjetflavor to 6 (--forbidden with matching--)
4) Fix a problem for the MadWidth generation in cluster mode. (the TMP directory use was not accessible by the cluster-> crash)

287. By Olivier Mattelaer on 2014-03-18

Fix a problem with the selection of the number of node to use in multicore

288. By Valentin Hirschi on 2014-03-18

1. Merged with the branch implementing flexible loop output dependencies
   via the new MG5_aMC option 'output_dependencies'

289. By Olivier Mattelaer on 2014-03-19

fix a problem with PBS cluster (one more actually)

290. By Paolo Torrielli on 2014-03-19

merged the lepton_masses branch

291. By Olivier Mattelaer on 2014-03-19

a patch from Pavel for STDHEP

292. By Rikkert Frederix on 2014-03-19

Small fix in madfks_plot: during the grid setup of fixed order runs,
the doreweight=.false., therefore there was a list of error messages
written to the screen saying that some ratio is not defined
correctly. The code did give the correct results: the fix is simply to
remove the error messages.

293. By Valentin Hirschi on 2014-03-19

1. Merged with new IOTests

294. By Valentin Hirschi on 2014-03-19

1. Fixed issues in IOTests because of border effects. They are now
   completely gotten rid of by decorating the IOTests run functions
   with the 'set_global' decorator.

295. By Valentin Hirschi on 2014-03-21

1. Replaced avh_olo.f90 with the new version provided by Andreas which fixes the
   analytic continuation issue I spotted in the box scalar integral during my checks
   for the ML standalone processes for the paper
   All loop parallel tests still pass.
   I made sure the dilog warning are muted.

296. By Olivier Mattelaer on 2014-03-22

fix a bug related to SysCalc

297. By Paolo Torrielli on 2014-03-22

trivial change in amcatnlo_run_interface.py, in order not to break the PDFCODE
functionality in the shower if lhapdf is linked

298. By marco zaro on 2014-03-24

merged with 2.1.1-lhapdf

299. By Rikkert Frederix on 2014-03-25

Fixed a bug that the photon rapidity cut was only applied to the
*first* photon in the process. This bug was introduced in r282 (merge
with the understanding_cuts branch) and was never part of the official
release.

300. By marco zaro on 2014-03-25

fixed stupid bug in lhapdf functions

301. By marco zaro on 2014-03-25

added a better explaination about MSTW 90% pdf sets in lhapdf

302. By Olivier Mattelaer on 2014-03-25

merge with the plugin model

303. By Rikkert Frederix on 2014-03-26

merge with the scale_pdf_fixed_order_print branch

[Rikkert Frederix (frederix)]

I think most things that we want to include are there. I've been running many NLO processes locally, and never encountered a problem. So: approve.

review: Approve

[Olivier Mattelaer (olivier-mattelaer)]

Actually I'm waiting for the (very important branch)
2.1.1_dynscale.

This branch is intended very important bug in LO since it basically affects all SM processes with 2j (basically) They make the jobs crashing so at least this are not returning pointless results but this should be fixed.

This should be now in principle fixed (I should thanks Simon for the huge help on this), but we face some additional problem in the heft (We think due to the hgg vertex not handle correctly).
So we still have to work on it.

Cheers,

Olivier

review: Needs Fixing

[lp:~maddevelopers/mg5amcnlo/2.1.1]

304. By Rikkert Frederix on 2014-03-27

Small fix for fxfx merging

305. By marco zaro on 2014-03-27

run_mcatnlo function refactored a bit (in amcatnlo_run_interface)

306. By Rikkert Frederix on 2014-03-27

small fix in setcuts: for processes with t-channels some BW were
incorrectly flagged as 'can-never-go-on-shell'.

307. By Olivier Mattelaer on 2014-03-29

Fix Bug #1298515 (MadSpin with @ in the process)
Fix Bug #1298521 (Gridpack not returning the expected number of event)
Fix Question #246123 (crash on model if none particles as a zero mass/width)

308. By Olivier Mattelaer on 2014-03-29

merge with the branch 2.1.1_dynscale fixing problem with the QED/QCD classification of jet

309. By Olivier Mattelaer on 2014-03-29

remove a print statement

310. By Olivier Mattelaer on 2014-03-29

restore the previous version of load_para, looks like this function is still use in a lot of place.

311. By marco zaro on 2014-03-30

fixed warnings concerning PDFsets when doing e+e- processes

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2014-02-20 09:49:16 +0000
+++ MadSpin/decay.py	2014-03-30 11:08:12 +0000
@@ -1869,7 +1869,6 @@
         # input
         self.options = options
         #max_weight_arg = options['max_weight']  
-        #BW_effects = options['BW_effect']
         self.path_me = os.path.realpath(options['curr_dir']) 
         if options['ms_dir']:
             self.path_me = os.path.realpath(options['ms_dir'])
@@ -2012,7 +2011,6 @@
         """Running the full code""" 
     
         max_weight_arg = self.options['max_weight']  
-        BW_effects = self.options['BW_effect']
         decay_tools=decay_misc()
         
         #Next step: we need to determine which matrix elements are really necessary
@@ -2334,14 +2332,14 @@
             decay_mapping = self.get_process_identical_ratio(relation)
             return decay_mapping
         
-        BW_cut = self.options['BW_cut'] if self.options['BW_effect'] else 1e-6        
+        BW_cut = self.options['BW_cut']       
         
         #class the decay by class (nbody/pid)
         nbody_to_decay = collections.defaultdict(list)
         for decay in self.all_decay.values():
             id = decay['dc_branch']['tree'][-1]['label']
             id_final = decay['processes'][0].get_final_ids_after_decay()
-            cut = self.options['zeromass_for_max_weight']
+            cut = 0.0 
             mass_final = tuple([m if m> cut else 0 for m in map(self.pid2mass, id_final)])
             
             nbody_to_decay[(decay['nbody'], abs(id), mass_final)].append(decay)
@@ -2603,7 +2601,7 @@
             if '@' in proc:
                 proc, proc_nb = proc.split('@')
                 try:
-                    int(proc_nb)
+                    proc_nb = int(proc_nb)
                 except ValueError:
                     raise MadSpinError, 'MadSpin didn\'t allow order restriction after the @ comment: \"%s\" not valid' % proc_nb
                 proc_nb = '@ %i' % proc_nb 

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2014-02-16 01:00:13 +0000
+++ MadSpin/interface_madspin.py	2014-03-30 11:08:12 +0000
@@ -65,12 +65,11 @@
         self.decay = madspin.decay_misc()
         self.model = None
         
-        self.options = {'max_weight': -1, 'BW_effect': 1, 
+        self.options = {'max_weight': -1, 
                         'curr_dir': os.path.realpath(os.getcwd()),
                         'Nevents_for_max_weigth': 0,
                         'max_weight_ps_point': 400,
                         'BW_cut':-1,
-                        'zeromass_for_max_weight':5,
                         'nb_sigma':0,
                         'ms_dir':None,
                         'max_running_process':100}
@@ -509,6 +508,13 @@
         generate_all.mscmd = self 
         generate_all.pid2width = lambda pid: generate_all.banner.get('param_card', 'decay', abs(pid)).value
         generate_all.pid2mass = lambda pid: generate_all.banner.get('param_card', 'mass', abs(pid)).value
+        if generate_all.path_me != self.options['ms_dir']:
+            for decay in generate_all.all_ME.values():
+                decay['path'] = decay['path'].replace(generate_all.path_me, self.options['ms_dir'])
+                for decay2 in decay['decays']:
+                    decay2['path'] = decay2['path'].replace(generate_all.path_me, self.options['ms_dir'])
+            generate_all.path_me = self.options['ms_dir'] # directory can have been move
+            generate_all.ms_dir = generate_all.path_me
         
         if not hasattr(self.banner, 'param_card'):
             self.banner.charge_card('slha')

=== modified file 'MadSpin/src/driver.f'
--- MadSpin/src/driver.f	2014-02-18 16:35:27 +0000
+++ MadSpin/src/driver.f	2014-03-30 11:08:12 +0000
@@ -610,9 +610,9 @@
       integer  map_external2res(nexternal_prod) ! map (index in production) -> index in the full structure
       double precision p(0:3,-nexternal:nexternal)
  
-      integer idB, id1, index_p2
+      integer idB, id1, index_p2, which_initial, last_branch, d1,d2
       double precision pa(0:3), pb(0:3), p1(0:3), p2(0:3),pboost(0:3)
-      double precision pb_cms(0:3), p1_cms(0:3), p1_rot(0:3)
+      double precision pb_cms(0:3), p1_cms(0:3), p1_rot(0:3),temp
 
 c     common
       ! variables to keep track of the vegas numbers for the production part
@@ -641,6 +641,7 @@
       include  'configs_decay.inc'
 
       mapext2res=map_external2res
+      which_initial=-1
 c Determine number of s- and t-channel branches, at this point it
 c includes the s-channel p1+p2
 c      write(*,*) (itree(i,-1), i=1,2)
@@ -708,7 +709,7 @@
       enddo
 
       ! fill the new itree with the kinematics associated with the production
-      do i=-1,-(ns_channel+nt_channel)-1,-1
+      do i=-1,-(ns_channel+nt_channel)-1,-1  ! loop over invariants in the production
 c         write(*,*) 'i prod',i
 c         write(*,*) 'i full',i-ns_channel_decay
 c         write(*,*) 'd1 prod',itree(1,i)
@@ -717,12 +718,14 @@
 c            write(*,*) 'd1 full',itree(1,i)-ns_channel_decay
             itree_full(1,i-ns_channel_decay) = itree(1,i)-ns_channel_decay
          else 
+             if (itree(1,i).le.2) which_initial=itree(1,i)
              itree_full(1,i-ns_channel_decay) = map_external2res(itree(1,i))
          endif 
          if (itree(2,i).lt.0 ) then
 c            write(*,*) 'd2 full',itree(2,i)-ns_channel_decay
              itree_full(2,i-ns_channel_decay) = itree(2,i)-ns_channel_decay
          else 
+             if (itree(2,i).le.2) which_initial=itree(2,i)
              itree_full(2,i-ns_channel_decay) = map_external2res(itree(2,i))
          endif
          
@@ -742,10 +745,58 @@
            qmass_full(i-ns_channel_decay)=qmass(i)
            qwidth_full(i-ns_channel_decay)=qwidth(i)
          endif
+
+
       enddo 
 
-
-
+c
+c        MODIF March 5, 2014 (P.A.) 
+c        overwrite last t-channel invariant to avoid numerical unstabilities
+c        This is crucial if the last t-channel propa is connected to two massless particles
+c        (one inital + one final particles) which are colinear -> t ~ 0
+c        Indeed, in that case the extraction of t in the previous loop may return a positive 
+c        value: in that case, the code will always fail to generate momenta !!!
+
+         if (nt_channel.ge.2) then
+           last_branch=-ns_channel-nt_channel-1
+           if (itree(1, last_branch).ge.0.or.itree(2, last_branch).ge.0) then
+           ! only overwrite if last t-channel splitting is connected to initial particle + massless particle
+           if (itree(1, last_branch).ge.0) d1=itree(1, last_branch)
+           if (itree(2, last_branch).ge.0) d1=itree(2, last_branch)
+           if (which_initial.eq.1) then 
+              d2=2
+           elseif (which_initial.eq.2) then
+              d2=1
+           else
+              write(*,*) 'Problem with the last t-channel branching ',which_initial 
+              stop
+           endif 
+
+           temp=-2d0*dot(p_ext(0,d1),p_ext(0,d2))
+           if (map_external2res(d1).gt.0) then
+           temp=temp+m(map_external2res(d1))**2
+           else
+           temp=temp+qmass_full(map_external2res(d1))**2
+           endif
+           temp=temp+m(map_external2res(d2))**2
+
+           shat=2d0*dot(p_ext(0,1),p_ext(0,2))
+
+c           write(*,*) temp, fixedinv(-ns_channel-nt_channel-ns_channel_decay)
+
+c           pause
+           if (abs(temp-fixedinv(-ns_channel-nt_channel-ns_channel_decay)).lt.(shat/1d2)) then 
+           fixedinv(-ns_channel-nt_channel-ns_channel_decay)=temp
+           else
+c              write(*,*)  qmass_full(map_external2res(d1)), m(map_external2res(d2))
+c              write(*,*) d1, d2, itree(1, last_branch), itree(2, last_branch)
+              write(*,*) 'cannot extract last t-channel invariant', temp, fixedinv(-ns_channel-nt_channel-ns_channel_decay)
+              stop
+           endif
+           endif
+         endif
+
+c     END MODIF MARCH 5, 2014
  
       !write(*,*) -ns_channel-nt_channel-1
       !write(*,*) map_external2res(itree(2,-ns_channel-nt_channel-1))

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2013-12-14 22:59:27 +0000
+++ Template/LO/Cards/run_card.dat	2014-03-30 11:08:12 +0000
@@ -48,7 +48,8 @@
 #*********************************************************************
 # PDF CHOICE: this automatically fixes also alpha_s and its evol.    *
 #*********************************************************************
- 'cteq6l1'    = pdlabel     ! PDF set                                     
+ 'cteq6l1'    = pdlabel     ! PDF set (lhapdf for using lhapdf)
+  10042        = lhaid      ! PDF number used ONLY for LHAPDF
 #*********************************************************************
 # Renormalization and factorization scales                           *
 #*********************************************************************
@@ -115,7 +116,7 @@
  -1  = ptlmax    ! maximum pt for the charged leptons
  -1  = missetmax ! maximum missing Et (sum of neutrino's momenta)
 #*********************************************************************
-# Minimum and maximum E's (in the lab frame)                         *
+# Minimum and maximum E's (in the center of mass frame)              *
 #*********************************************************************
   0  = ej     ! minimum E for the jets 
   0  = eb     ! minimum E for the b 
@@ -162,6 +163,8 @@
  -1  = dralmax ! maxdistance between gamma and lepton
 #*********************************************************************
 # Minimum and maximum invariant mass for pairs                       *
+# WARNING: for four lepton final state mmll cut require to have      *
+#          different lepton masses for each flavor!                  *           
 #*********************************************************************
  0   = mmjj    ! min invariant mass of a jet pair 
  0   = mmbb    ! min invariant mass of a b pair 
@@ -263,12 +266,12 @@
 #*********************************************************************
    F  = use_syst      ! Enable systematics studies
 #
-#**************************************                                       
-# Parameter of the systematics study 
-#  will be use by SysCalc (if install)
-#**************************************                                     
+#**************************************
+# Parameter of the systematics study
+#  will be used by SysCalc (if installed)
+#**************************************                                  
 #
-0.5 1 2 = sys_scalefact  # Central scale factors
+0.5 1 2 = sys_scalefact  # factorization/renormalization scale factor
 0.5 1 2 = sys_alpsfact  # \alpha_s emission scale factors
 30 50 = sys_matchscale # variation of merging scale
 # PDF sets and number of members (0 or none for all members).

=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2013-11-24 16:38:02 +0000
+++ Template/LO/Source/dsample.f	2014-03-30 11:08:12 +0000
@@ -970,6 +970,7 @@
       endif
 
       im = ddum(j)
+      if (im.ge.ng) im = ng -1
       ip = im + 1
       ij = Minvar(j,ipole)
 c------

=== modified file 'Template/LO/Source/rw_routines.f'
--- Template/LO/Source/rw_routines.f	2012-06-23 03:07:49 +0000
+++ Template/LO/Source/rw_routines.f	2014-03-30 11:08:12 +0000
@@ -1,3 +1,107 @@
+      subroutine load_gridpack_para(npara,param,value)
+c----------------------------------------------------------------------
+c Read the params from the run_card.dat file
+c----------------------------------------------------------------------
+      implicit none
+c
+c     arguments
+c
+      character*20 param(*),value(*)
+      integer npara
+c
+c     local
+c
+      logical fopened,done
+      integer iunit
+      character*20 ctemp
+      integer k,i,l1,l2,iproc
+      character*132 buff
+      data iunit/21/
+c
+c     global
+c
+      integer ngroup
+      common/to_group/ngroup
+c
+c----------
+c     start
+c----------
+      npara=0
+      param(1)=' '
+      value(1)=' '
+c
+c     open file
+c
+      call open_file(iunit,'grid_card.dat',fopened)
+      if(fopened) then
+c
+c     first look for process-specific parameters
+c
+      done=.false.
+      do while(.not.done)
+         read(iunit,'(a132)',end=30,err=30) buff
+         if(buff(1:1).ne.'#' .and. index(buff,"=").gt.0
+     $        .and. index(buff,"@").gt.0) then
+            l1=index(buff,"@")
+            l2=index(buff,"!")
+            if(l2.eq.0) l2=l1+20  !maybe there is no comment...
+            read(buff(l1+1:l2),*,err=21) iproc
+            if(iproc.ne.ngroup) cycle
+
+            l1=index(buff,"=")
+            l2=index(buff,"@")
+            if(l2-l1.lt.0) cycle
+            npara=npara+1
+c
+             value(npara)=buff(1:l1-1)
+             ctemp=value(npara)
+             call case_trap2(ctemp)
+             value(npara)=ctemp
+c
+             param(npara)=" "//buff(l1+1:l2-1)
+             ctemp=param(npara)
+             call case_trap2(ctemp)
+             param(npara)=ctemp
+c
+ 21          cycle
+         endif
+       enddo
+ 30   rewind(iunit)
+c
+c     read in values
+c
+      done=.false.
+      do while(.not.done)
+         read(iunit,'(a132)',end=99,err=99) buff
+         if(buff(1:1).ne.'#' .and. index(buff,"=").gt.0
+     $        .and. index(buff,"@").le.0) then
+            l1=index(buff,"=")
+            l2=index(buff,"!")
+            if(l2.eq.0) l2=l1+20  !maybe there is no comment...
+            if(l2-l1.lt.0) cycle
+            npara=npara+1
+c
+             value(npara)=buff(1:l1-1)
+             ctemp=value(npara)
+             call case_trap2(ctemp)
+             value(npara)=ctemp
+c
+             param(npara)=" "//buff(l1+1:l2-1)
+c             write (*,*) param(npara),l1,l2
+             ctemp=param(npara)
+             call case_trap2(ctemp)
+             param(npara)=ctemp
+c             write(*,*) "New param:",param(npara)," = ", value(npara)
+c
+         endif
+      enddo
+ 99   close(iunit)
+      endif
+
+      return
+      end
+
+
       subroutine load_para(npara,param,value)
 c----------------------------------------------------------------------
 c Read the params from the run_card.dat file

=== modified file 'Template/LO/SubProcesses/cuts.f'
--- Template/LO/SubProcesses/cuts.f	2014-01-29 15:48:36 +0000
+++ Template/LO/SubProcesses/cuts.f	2014-03-30 11:08:12 +0000
@@ -26,11 +26,8 @@
 C     INPUT:
 C            P(0:3,1)           MOMENTUM OF INCOMING PARTON
 C            P(0:3,2)           MOMENTUM OF INCOMING PARTON
-C            P(0:3,3)           MOMENTUM OF d
-C            P(0:3,4)           MOMENTUM OF b
-C            P(0:3,5)           MOMENTUM OF bbar
-C            P(0:3,6)           MOMENTUM OF e+
-C            P(0:3,7)           MOMENTUM OF ve
+C            P(0:3,3)           MOMENTUM OF ...
+C            ALL MOMENTA ARE IN THE REST FRAME!!
 C            COMMON/JETCUTS/   CUTS ON JETS
 C     OUTPUT:
 C            TRUE IF EVENTS PASSES ALL CUTS LISTED

=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2014-02-05 23:28:12 +0000
+++ Template/LO/SubProcesses/reweight.f	2014-03-30 11:08:12 +0000
@@ -191,6 +191,14 @@
 c**************************************************
 c   Traces particle lines according to CKKW rules
 c**************************************************
+c     ipart gives the external particle number corresponding to the present
+c     quark or gluon line.
+c     For t-channel lines, ipart(1) contains the connected beam.
+c     For s-channel lines, it depends if it is quark or gluon line:
+c     For quark lines, ipart(2) is 0 and ipart(1) connects to the corresponding
+c     final-state quark. For gluons, if it splits into two gluons,
+c     it connects to the hardest gluon. If it splits into qqbar, it ipart(1) is
+c     the hardest and ipart(2) is the softest.
       implicit none
 
       include 'ncombs.inc'
@@ -202,6 +210,9 @@
       integer ipdg(n_max_cl),ipart(2,n_max_cl)
       logical isjet
       external isjet
+      integer iddgluon, iddother, idgluon, idother
+      logical isqcd
+      external isqcd
 
       idmo=ipdg(imo)
       idda1=ipdg(ida1)
@@ -273,26 +284,54 @@
           ipart(1,imo)=ipart(1,ida2)
           ipart(2,imo)=ipart(2,ida2)
         endif
-      else if(idmo.eq.21)then
-      write(*,*) '277',ida1, ida2,isjet(ida1),isjet(ida2)
+      else if(idmo.eq.21 .and. abs(idda1).le.6 .and.
+     $        abs(idda2).le.6) then
 c     gluon -> quark anti-quark: use both, but take hardest as 1
-c     BUT if only one is a jet
-        if(isjet(ida1).and..not.isjet(ida2)) then
-          ipart(1,imo)=ipart(1,ida1)
-          ipart(2,imo)=ipart(1,ida2)
-        else if(isjet(ida2).and..not.isjet(ida1)) then
-          ipart(1,imo)=ipart(1,ida2)
-          ipart(2,imo)=ipart(1,ida1)        
-        else if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
+        if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
      $     p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
-          write(*,*) '1 wins'
           ipart(1,imo)=ipart(1,ida1)
           ipart(2,imo)=ipart(1,ida2)
         else
-          write(*,*) '2 wins'
           ipart(1,imo)=ipart(1,ida2)
           ipart(2,imo)=ipart(1,ida1)
         endif
+      else if(idmo.eq.21.and.(idda1.eq.21.or.idda2.eq.21))then
+         if(idda1.eq.21) then
+            iddgluon = idda1
+            idgluon = ida1
+            iddother = idda2
+            idother = ida2
+         else
+            iddgluon = idda2
+            iddother = idda1
+            idgluon = ida2
+            idother = ida1
+         endif
+         if (isqcd(idother))then
+c        gluon -> gluon + scalar octet Choose hardest one
+            if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
+     $         p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
+               ipart(1,imo)=ipart(1,ida1)
+               ipart(2,imo)=ipart(2,ida1)
+            else
+               ipart(1,imo)=ipart(1,ida2)
+               ipart(2,imo)=ipart(2,ida2)
+            endif
+         else
+c        gluon -> gluon + Higgs use the gluon one
+               ipart(1,imo)=ipart(1,idgluon)
+               ipart(2,imo)=ipart(2,idgluon)
+         endif
+      else if(idmo.eq.21) then
+c     gluon > octet octet Choose hardest one
+            if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
+     $         p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
+               ipart(1,imo)=ipart(1,ida1)
+               ipart(2,imo)=ipart(2,ida1)
+            else
+               ipart(1,imo)=ipart(1,ida2)
+               ipart(2,imo)=ipart(2,ida2)
+            endif
       else if(idmo.eq.idda1.or.idmo.eq.idda1+sign(1,idda2))then
 c     quark -> quark-gluon or quark-Z or quark-h or quark-W
         ipart(1,imo)=ipart(1,ida1)
@@ -331,7 +370,6 @@
       idmo=ipdg(imo)
       idda1=ipdg(ida1)
       idda2=ipdg(ida2)
-
 c     Check QCD vertex
       if(islast.or..not.isqcd(idmo).or..not.isqcd(idda1).or.
      &     .not.isqcd(idda2)) then
@@ -353,13 +391,12 @@
       endif        
 
 c     FS clustering
-      if(isjet(idda1).and.(isjet(idmo).or.idmo.eq.idda2).or.
-     $   isjet(idda2).and.(isjet(idmo).or.idmo.eq.idda1)) then
+      if((isjet(idda1).and.(isjet(idmo).or.idmo.eq.idda2)).or.
+     $   (isjet(idda2).and.(isjet(idmo).or.idmo.eq.idda1))) then
          isjetvx=.true.
       else
          isjetvx=.false.
       endif
-      
       return
       end
 
@@ -493,6 +530,7 @@
       integer fsnum(2),ida(2),imo,jcode
       logical chclusold,fail,increasecode
       save chclusold
+      integer tmpindex
 
       logical isqcd,isjet,isparton,cluster,isjetvx
       integer ifsno
@@ -577,21 +615,24 @@
         call ipartupdate(p,imocl(n),idacl(n,1),idacl(n,2),
      $       ipdgcl(1,igraphs(1),iproc),ipart)
       enddo
-      
+
 c     Prepare beam related variables for scale and jet determination
       do i=1,2
          ibeam(i)=ishft(1,i-1)
 c        jfirst is first parton splitting on this side
          jfirst(i)=0
-c        jlast is last parton on this side
+c        jlast is last parton on this side This means
+c        the last cluster which is still QCD.
          jlast(i)=0
-c        jcentral is the central scale vertex on this side
+c        jcentral is the central scale vertex on this side. i.e it stops
+c        when the T channel particles is not colored anymore.
          jcentral(i)=0
 c        qcdline gives whether this IS line is QCD
          qcdline(i)=isqcd(ipdgcl(ibeam(i),igraphs(1),iproc))
 c        partonline gives whether this IS line is parton (start out true for any QCD)
          partonline(i)=qcdline(i)
 c        goodjet gives whether this cluster line is considered a jet
+c        i.e. if all related/previous clustering are jet
          goodjet(ibeam(i))=partonline(i)
       enddo
 
@@ -611,8 +652,8 @@
 c     increasecode gives whether we should increase jcode at next vertex
       increasecode=.false.
       do n=1,nexternal-2
-        do i=1,2
-          do j=1,2
+        do i=1,2 ! index of the child in the interaction
+          do j=1,2 ! j index of the beam
             if(idacl(n,i).eq.ibeam(j))then
 c             IS clustering
               ibeam(j)=imocl(n)
@@ -634,6 +675,9 @@
                  jlast(j)=n
                  partonline(j)=goodjet(ida(3-i)).and.
      $                isjet(ipdgcl(imo,igraphs(1),iproc))
+              else if (jfirst(j).eq.0) then
+                 jfirst(j) = n
+                 goodjet(imo)=.false.
               else
                  goodjet(imo)=.false.
               endif
@@ -653,7 +697,7 @@
 c             FS line is a jet line
               if(goodjet(ida(3-i))) then
                  if(partonline(j).or.
-     $              ipdgcl(ida(3-i),igraphs(1),iproc).eq.21)then
+     $ ipdgcl(ida(3-i),igraphs(1),iproc).eq.21)then
 c                   Need to include gluon to avoid soft singularity
                     iqjets(ipart(1,ida(3-i)))=1 ! 1 means for sure jet
                  else
@@ -674,20 +718,21 @@
            if(.not.isjetvx(imocl(n),idacl(n,1),idacl(n,2),
      $        ipdgcl(1,igraphs(1),iproc),ipart,n.eq.nexternal-2)) then
 c          Remove non-gluon jets that lead up to non-jet vertices
-              if(ipart(1,imocl(n)).gt.2)then ! ipart(1) set and not IS line
-c                The ishft gives the FS particle corresponding to imocl
-                 if(ipdgcl(ishft(1,ipart(1,imocl(n))-1),igraphs(1),iproc).ne.21)
+           if(ipart(1,imocl(n)).gt.2)then ! ipart(1) set and not IS line
+c          The ishft gives the FS particle corresponding to imocl
+              if(ipdgcl(ishft(1,ipart(1,imocl(n))-1),igraphs(1),iproc).ne.21)
      $              iqjets(ipart(1,imocl(n)))=0
-              endif
-              if(ipart(2,imocl(n)).gt.2)then ! ipart(1) set and not IS line
-c                The ishft gives the FS particle corresponding to imocl
-                 if(ipdgcl(ishft(1,ipart(2,imocl(n))-1),igraphs(1),iproc).ne.21)
-     $             iqjets(ipart(2,imocl(n)))=0
-              endif
+           endif
+           if(ipart(2,imocl(n)).gt.2)then ! ipart(1) set and not IS line
+c             The ishft gives the FS particle corresponding to imocl
+              if(ipdgcl(ishft(1,ipart(2,imocl(n))-1),igraphs(1),iproc).ne.21)
+     $              iqjets(ipart(2,imocl(n)))=0
+           endif
 c          Set goodjet to false for mother
               goodjet(imocl(n))=.false.
               cycle
            endif
+
 c          This is a jet vertex, so set jet flag for final-state jets
 c          ifsno gives leg number if daughter is FS particle, otherwise 0
            fsnum(1)=ifsno(idacl(n,1),ipart)
@@ -703,7 +748,7 @@
 c          Flag mother as good jet if PDG is jet and both daughters are jets
            goodjet(imocl(n))=
      $          (isjet(ipdgcl(imocl(n),igraphs(1),iproc)).and.
-     $          goodjet(idacl(n,1)).and.goodjet(idacl(n,1)))
+     $          goodjet(idacl(n,1)).and.goodjet(idacl(n,2)))
         endif
       enddo
 
@@ -782,7 +827,7 @@
      $           write(*,*) 'Bad clustering, jets fail. Reclustering ',
      $           iconfig
             chcluster=.true.
-            goto 100
+c            goto 100 ! not
          endif
       endif
       
@@ -936,13 +981,17 @@
             q2fact(1)=pt2ijcl(nexternal-2)
             q2fact(2)=q2fact(1)
          endif
+      elseif(jcentral(1).eq.0)then
+            q2fact(1) = pt2ijcl(jfirst(1))
+      elseif(jcentral(2).eq.0)then
+            q2fact(2) = pt2ijcl(jfirst(2))
       elseif(ickkw.eq.2.or.pdfwgt)then
 c     Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
 c     f1(x1,pt2E) is given by DSIG, just need to set scale.
 c     Use the minimum scale found for fact scale in ME
-         if(jlast(1).gt.0.and.jfirst(1).lt.jlast(1))
+         if(jlast(1).gt.0.and.jfirst(1).le.jlast(1))
      $        q2fact(1)=min(pt2ijcl(jfirst(1)),q2fact(1))
-         if(jlast(2).gt.0.and.jfirst(2).lt.jlast(2))
+         if(jlast(2).gt.0.and.jfirst(2).le.jlast(2))
      $        q2fact(2)=min(pt2ijcl(jfirst(2)),q2fact(2))
       endif
 
@@ -1118,6 +1167,7 @@
          write(*,*) 'Set process number ',ipsel
       endif
 
+      if (use_syst.and.igraphs(1).eq.0) igraphs(1) = 1 ! happens if use_syst=T BUT fix scale
 c     Set incoming particle identities
       ipdgcl(1,igraphs(1),iproc)=idup(1,ipsel,iproc)
       ipdgcl(2,igraphs(1),iproc)=idup(2,ipsel,iproc)

=== modified file 'Template/LO/SubProcesses/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2013-11-29 01:25:36 +0000
+++ Template/LO/SubProcesses/setcuts.f	2014-03-30 11:08:12 +0000
@@ -178,7 +178,7 @@
          if (abs(idup(i,1,1)).eq.14) do_cuts(i)=.false.  ! no cuts on vm vm~
          if (abs(idup(i,1,1)).eq.16) do_cuts(i)=.false.  ! no cuts on vt vt~
 c-flavor-jets
-         if (abs(idup(i,1,1)).le.min(maxjetflavor,5)) then
+         if (abs(idup(i,1,1)).le.min(maxjetflavor,6)) then
               is_a_j(i)=.true.
 c              write(*,*)'jet:ithe pdg is ',abs(idup(i,1,1)),' maxflavor=',maxjetflavor
          else if (abs(idup(i,1,1)).ge.maxjetflavor+1 .and. abs(idup(i,1,1)).le.5) then
@@ -458,7 +458,7 @@
       do i=nincoming+1,nexternal-1
       do j=nincoming+1,nexternal-1
          if(j.lt.i)then
-            s_min(i,j) = s_min(j,i)
+            s_min(i,j) = max(s_min(j,i),s_min(i,j))
          else
             smin=0.0d0**2
 

=== modified file 'Template/LO/bin/internal/gen_jpeg-pl'
--- Template/LO/bin/internal/gen_jpeg-pl	2012-02-24 23:14:05 +0000
+++ Template/LO/bin/internal/gen_jpeg-pl	2014-03-30 11:08:12 +0000
@@ -62,7 +62,7 @@
 #
 #  Now create jpeg file for card
 #
-  if (! -e "../../card.jpg") {
+  if (! -e "../../HTML/card.jpg") {
     system ("/bin/bash -c \"head -352 matrix$imatrix.ps >& junk.ps\" ");
     open(JUNK,">> junk.ps") || die "Error opening junk.ps";  
 
@@ -72,7 +72,7 @@
 
     system ("/bin/bash -c \"cat matrix$imatrix.ps | sed 1,352d >> junk.ps\" ");
 
-    system "/bin/bash -c \"nice gs  \-sDEVICE\=jpeg \-sOutputFile\=card.jpg \-q \-dNOPAUSE \-dBATCH \-g180x150 ./junk.ps; rm -f junk.ps; cp -p card.jpg ../../card.jpg\" ";
+    system "/bin/bash -c \"nice gs  \-sDEVICE\=jpeg \-sOutputFile\=card.jpg \-q \-dNOPAUSE \-dBATCH \-g180x150 ./junk.ps; rm -f junk.ps; cp -p card.jpg ../../HTML/card.jpg\" ";
   }
   if ($imatrix eq "") {$imatrix = 0;}
   $imatrix = $imatrix + 1;

=== modified file 'Template/LO/bin/internal/syscalc_template.dat'
--- Template/LO/bin/internal/syscalc_template.dat	2013-11-21 18:58:16 +0000
+++ Template/LO/bin/internal/syscalc_template.dat	2014-03-30 11:08:12 +0000
@@ -1,6 +1,9 @@
 # Central scale factors
 scalefact:
 %(sys_scalefact)s
+# choice of correlation scheme between muF and muR
+scalecorrelation:
+%(sys_scalecorrelation)s
 # \alpha_s emission scale factors
 alpsfact:
 %(sys_alpsfact)s

=== removed file 'Template/MadWeight/src/phasespace.inc'
--- Template/MadWeight/src/phasespace.inc	2014-01-14 10:43:25 +0000
+++ Template/MadWeight/src/phasespace.inc	1970-01-01 00:00:00 +0000
@@ -1,13 +0,0 @@
-	  include 'maxparticles.inc' 
-      integer max_branches
-      parameter (max_branches=max_particles-1)
-      integer max_configs
-      parameter (max_configs=100)
-
-c     channel position
-      integer config_pos,perm_pos
-      common /to_config/config_pos,perm_pos
-
-      integer max_channel
-      parameter (max_channel=10000)
-

=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2014-01-24 16:51:50 +0000
+++ Template/NLO/Cards/run_card.dat	2014-03-30 11:08:12 +0000
@@ -137,7 +137,8 @@
  1.0  = epsgamma  ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
  .true.  = isoEM  ! isolate photons from EM energy (photons and leptons)
 #***********************************************************************
-# maximal pdg code for quark to be considered as a jet                 *
+# Maximal PDG code for quark to be considered a jet when applying cuts.*
+# At least all massless quarks of the model should be included here.   *
 #***********************************************************************
- 5 = maxjetflavor
+ 4 = maxjetflavor
 #***********************************************************************

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_V.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_V.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_V.f	2014-03-30 11:08:12 +0000
@@ -21,7 +21,7 @@
       PARAMETER (PI=3.14159265358979312D0)
       integer j,kk,l,jpr
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_h.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_h.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_h.f	2014-03-30 11:08:12 +0000
@@ -26,39 +26,39 @@
       l=(kk-1)*40
       call mbook(l+1,'Higgs pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+2,'Higgs pT '//weights_info(kk)
+      call mbook(l+2,'Higgs pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+3,'Higgs log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+4,'Higgs pT,|y_H|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+5,'Higgs pT,|y_H|<2 '//weights_info(kk)
+      call mbook(l+5,'Higgs pT 1,|y_H|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+6,'Higgs log[pT],|y_H|<2 '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
 
       call mbook(l+7,'j1 pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+8,'j1 pT '//weights_info(kk)
+      call mbook(l+8,'j1 pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+9,'j1 log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+10,'j1 pT,|y_j1|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+11,'j1 pT,|y_j1|<2 '//weights_info(kk)
+      call mbook(l+11,'j1 pT 1,|y_j1|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+12,'j1 log[pT],|y_j1|<2 '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
 
       call mbook(l+13,'Inc j pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+14,'Inc j pT '//weights_info(kk)
+      call mbook(l+14,'Inc j pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+15,'Inc j log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+16,'Inc j pT,|y_Ij|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,2.d2)
-      call mbook(l+17,'Inc j pT,|y_Ij|<2 '//weights_info(kk)
+      call mbook(l+17,'Inc j pT 1,|y_Ij|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,5.d2)
       call mbook(l+18,'Inc j log[pT],|y_Ij|<2'//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_hjj.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_hjj.f	2014-02-18 17:20:34 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_hjj.f	2014-03-30 11:08:12 +0000
@@ -42,7 +42,7 @@
 
          call mbook(l+  2,'Higgs pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+  3,'Higgs pT '
+         call mbook(l+  3,'Higgs pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+  4,'Higgs logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -53,7 +53,7 @@
 
          call mbook(l+  7,'j1 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+  8,'j1 pT '
+         call mbook(l+  8,'j1 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+  9,'j1 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -64,7 +64,7 @@
 
          call mbook(l+ 12,'j2 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 13,'j2 pT '
+         call mbook(l+ 13,'j2 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 14,'j2 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -75,7 +75,7 @@
 
          call mbook(l+ 17,'j3 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 18,'j3 pT '
+         call mbook(l+ 18,'j3 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 19,'j3 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -86,7 +86,7 @@
 
          call mbook(l+ 22,'H+j1 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 23,'H+j1 pT '
+         call mbook(l+ 23,'H+j1 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 24,'H+j1 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -97,7 +97,7 @@
 
          call mbook(l+ 27,'j1+j2 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 28,'j1+j2 pT '
+         call mbook(l+ 28,'j1+j2 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 29,'j1+j2 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -108,7 +108,7 @@
 
          call mbook(l+ 32,'syst pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 33,'syst pT '
+         call mbook(l+ 33,'syst pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 34,'syst logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -151,7 +151,7 @@
          call mbook(l+ 50,'jveto pT '
      &       //cc(i)//weights_info(kk),
      1    ((vetomax-vetomin)/dble(nbinveto)),vetomin,vetomax)
-         call mbook(l+ 51,'jveto pT '
+         call mbook(l+ 51,'jveto pT 1 '
      &       //cc(i)//weights_info(kk),
      1    ((2d0*vetomax-vetomin)/dble(nbinveto)),
      1    vetomin,2d0*vetomax)

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_lplm.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_lplm.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_lplm.f	2014-03-30 11:08:12 +0000
@@ -36,7 +36,7 @@
       l=(kk-1)*42+(j-1)*21
       call mbook(l+ 1,'V pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+ 2,'V pt      '//cc(j)//weights_info(kk)
+      call mbook(l+ 2,'V pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+ 3,'V log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)
@@ -49,7 +49,7 @@
 c
       call mbook(l+ 7,'lm pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+ 8,'lm pt      '//cc(j)//weights_info(kk)
+      call mbook(l+ 8,'lm pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+ 9,'lm log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)
@@ -57,7 +57,7 @@
      &     ,0.2d0,-9.d0,9.d0)
       call mbook(l+11,'lp pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+12,'lp pt      '//cc(j)//weights_info(kk)
+      call mbook(l+12,'lp pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+13,'lp log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_lvl.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_lvl.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_lvl.f	2014-03-30 11:08:12 +0000
@@ -21,7 +21,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_taptam.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_taptam.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_taptam.f	2014-03-30 11:08:12 +0000
@@ -16,7 +16,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_tj.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_tj.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_tj.f	2014-03-30 11:08:12 +0000
@@ -18,7 +18,7 @@
       include 'reweight0.inc'
       integer j,kk,l
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
@@ -27,7 +27,7 @@
       common/cwgtsinfo/weights_info
       call inihist
       nwgt_analysis=nwgt
-      do j=1,2
+      do j=1,1
       do kk=1,nwgt_analysis
       l=(kk-1)*48+(j-1)*24
       call mbook(l+ 1,'t pt        '//cc(j)//weights_info(kk)
@@ -323,7 +323,7 @@
       enddo
       nbjet=count_bj
 c fill the histograms
-      do i=1,2
+      do i=1,1
          do kk=1,nwgt_analysis
             l=(kk-1)*48+(i-1)*24
             call mfill(l+1,pttop,www(kk))

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_ttx.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_ttx.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_ttx.f	2014-03-30 11:08:12 +0000
@@ -55,9 +55,9 @@
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
         call mbook(l+12,'delta y          '
      &       //cc(i)//weights_info(kk),0.2d0,-4.d0,4.d0)
-        call mbook(l+13,'tt azimt         '
+        call mbook(l+13,'tt azimt 1       '
      &       //cc(i)//weights_info(kk),pi/60.d0,2*pi/3,pi)
-        call mbook(l+14,'tt del R         '
+        call mbook(l+14,'tt del R 1       '
      &       //cc(i)//weights_info(kk),pi/60.d0,2*pi/3,4*pi/3)
         call mbook(l+15,'y_tb             '
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
@@ -65,11 +65,11 @@
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
         call mbook(l+17,'tt log[pi-azimt] '
      &       //cc(i)//weights_info(kk),0.05d0,-4.d0,0.1d0)
-        call mbook(l+18,'tt pt            '
+        call mbook(l+18,'tt pt 1          '
      &       //cc(i)//weights_info(kk),20.d0,80.d0,2000.d0)
-        call mbook(l+19,'tb pt            '
+        call mbook(l+19,'tb pt 1          '
      &       //cc(i)//weights_info(kk),20.d0,400.d0,2400.d0)
-        call mbook(l+20,'t pt             '
+        call mbook(l+20,'t pt 1           '
      &       //cc(i)//weights_info(kk),20.d0,400.d0,2400.d0)
       enddo
       enddo

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_ttx_v2.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_ttx_v2.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_pp_ttx_v2.f	2014-03-30 11:08:12 +0000
@@ -18,7 +18,7 @@
       integer j,kk,l,i
       PARAMETER (PI=3.14159265358979312D0)
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_rates.f'
--- Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_rates.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/mcatnlo_hwan_rates.f	2014-03-30 11:08:12 +0000
@@ -12,7 +12,7 @@
       include 'reweight0.inc'
       integer i,kk,l
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_V_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_V_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_V_hepmc.f	2014-03-30 11:08:12 +0000
@@ -21,7 +21,7 @@
       PARAMETER (PI=3.14159265358979312D0)
       integer j,kk,l,jpr
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_h_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_h_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_h_hepmc.f	2014-03-30 11:08:12 +0000
@@ -32,39 +32,39 @@
       l=(kk-1)*40
       call mbook(l+1,'Higgs pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+2,'Higgs pT '//weights_info(kk)
+      call mbook(l+2,'Higgs pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+3,'Higgs log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+4,'Higgs pT,|y_H|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+5,'Higgs pT,|y_H|<2 '//weights_info(kk)
+      call mbook(l+5,'Higgs pT 1,|y_H|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+6,'Higgs log[pT],|y_H|<2 '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
 
       call mbook(l+7,'j1 pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+8,'j1 pT '//weights_info(kk)
+      call mbook(l+8,'j1 pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+9,'j1 log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+10,'j1 pT,|y_j1|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+11,'j1 pT,|y_j1|<2 '//weights_info(kk)
+      call mbook(l+11,'j1 pT 1,|y_j1|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+12,'j1 log[pT],|y_j1|<2 '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
 
       call mbook(l+13,'Inc j pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+14,'Inc j pT '//weights_info(kk)
+      call mbook(l+14,'Inc j pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+15,'Inc j log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+16,'Inc j pT,|y_Ij|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,2.d2)
-      call mbook(l+17,'Inc j pT,|y_Ij|<2 '//weights_info(kk)
+      call mbook(l+17,'Inc j pT 1,|y_Ij|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,5.d2)
       call mbook(l+18,'Inc j log[pT],|y_Ij|<2'//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_hjj_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_hjj_hepmc.f	2014-02-18 17:20:34 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_hjj_hepmc.f	2014-03-30 11:08:12 +0000
@@ -47,7 +47,7 @@
 
          call mbook(l+  2,'Higgs pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+  3,'Higgs pT '
+         call mbook(l+  3,'Higgs pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+  4,'Higgs logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -58,7 +58,7 @@
 
          call mbook(l+  7,'j1 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+  8,'j1 pT '
+         call mbook(l+  8,'j1 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+  9,'j1 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -69,7 +69,7 @@
 
          call mbook(l+ 12,'j2 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 13,'j2 pT '
+         call mbook(l+ 13,'j2 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 14,'j2 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -80,7 +80,7 @@
 
          call mbook(l+ 17,'j3 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 18,'j3 pT '
+         call mbook(l+ 18,'j3 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 19,'j3 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -91,7 +91,7 @@
 
          call mbook(l+ 22,'H+j1 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 23,'H+j1 pT '
+         call mbook(l+ 23,'H+j1 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 24,'H+j1 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -102,7 +102,7 @@
 
          call mbook(l+ 27,'j1+j2 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 28,'j1+j2 pT '
+         call mbook(l+ 28,'j1+j2 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 29,'j1+j2 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -113,7 +113,7 @@
 
          call mbook(l+ 32,'syst pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 33,'syst pT '
+         call mbook(l+ 33,'syst pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 34,'syst logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -156,7 +156,7 @@
          call mbook(l+ 50,'jveto pT '
      &       //cc(i)//weights_info(kk),
      1    ((vetomax-vetomin)/dble(nbinveto)),vetomin,vetomax)
-         call mbook(l+ 51,'jveto pT '
+         call mbook(l+ 51,'jveto pT 1 '
      &       //cc(i)//weights_info(kk),
      1    ((2d0*vetomax-vetomin)/dble(nbinveto)),
      1    vetomin,2d0*vetomax)

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_lplm_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_lplm_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_lplm_hepmc.f	2014-03-30 11:08:12 +0000
@@ -42,7 +42,7 @@
       l=(kk-1)*42+(j-1)*21
       call mbook(l+ 1,'V pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+ 2,'V pt      '//cc(j)//weights_info(kk)
+      call mbook(l+ 2,'V pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+ 3,'V log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)
@@ -55,7 +55,7 @@
 c
       call mbook(l+ 7,'lm pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+ 8,'lm pt      '//cc(j)//weights_info(kk)
+      call mbook(l+ 8,'lm pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+ 9,'lm log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)
@@ -63,7 +63,7 @@
      &     ,0.2d0,-9.d0,9.d0)
       call mbook(l+11,'lp pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+12,'lp pt      '//cc(j)//weights_info(kk)
+      call mbook(l+12,'lp pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+13,'lp log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_lvl_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_lvl_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_lvl_hepmc.f	2014-03-30 11:08:12 +0000
@@ -21,7 +21,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_taptam_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_taptam_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_taptam_hepmc.f	2014-03-30 11:08:12 +0000
@@ -16,7 +16,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_tj_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_tj_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_tj_hepmc.f	2014-03-30 11:08:12 +0000
@@ -20,7 +20,7 @@
       PARAMETER (PI=3.14159265358979312D0)
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
@@ -35,7 +35,7 @@
       weights_info(nwgt)="central value  "
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       nwgt_analysis=nwgt
-      do j=1,2
+      do j=1,1
       do kk=1,nwgt_analysis
       l=(kk-1)*48+(j-1)*24
       call mbook(l+ 1,'t pt        '//cc(j)//weights_info(kk)
@@ -313,7 +313,7 @@
       enddo
       nbjet=count_bj
 c fill the histograms
-      do i=1,2
+      do i=1,1
          do kk=1,nwgt_analysis
             l=(kk-1)*48+(i-1)*24
             call mfill(l+1,pttop,www(kk))

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_ttx_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_ttx_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_ttx_hepmc.f	2014-03-30 11:08:12 +0000
@@ -60,9 +60,9 @@
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
         call mbook(l+12,'delta y          '
      &       //cc(i)//weights_info(kk),0.2d0,-4.d0,4.d0)
-        call mbook(l+13,'tt azimt         '
+        call mbook(l+13,'tt azimt 1       '
      &       //cc(i)//weights_info(kk),pi/60.d0,2*pi/3,pi)
-        call mbook(l+14,'tt del R         '
+        call mbook(l+14,'tt del R 1       '
      &       //cc(i)//weights_info(kk),pi/60.d0,2*pi/3,4*pi/3)
         call mbook(l+15,'y_tb             '
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
@@ -70,11 +70,11 @@
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
         call mbook(l+17,'tt log[pi-azimt] '
      &       //cc(i)//weights_info(kk),0.05d0,-4.d0,0.1d0)
-        call mbook(l+18,'tt pt            '
+        call mbook(l+18,'tt pt 1          '
      &       //cc(i)//weights_info(kk),20.d0,80.d0,2000.d0)
-        call mbook(l+19,'tb pt            '
+        call mbook(l+19,'tb pt 1          '
      &       //cc(i)//weights_info(kk),20.d0,400.d0,2400.d0)
-        call mbook(l+20,'t pt             '
+        call mbook(l+20,'t pt 1           '
      &       //cc(i)//weights_info(kk),20.d0,400.d0,2400.d0)
       enddo
       enddo

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_ttx_v2_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_ttx_v2_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_pp_ttx_v2_hepmc.f	2014-03-30 11:08:12 +0000
@@ -18,7 +18,7 @@
       PARAMETER (PI=3.14159265358979312D0)
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_rates_hepmc.f'
--- Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_rates_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/HWPPAnalyzer/mcatnlo_hwan_rates_hepmc.f	2014-03-30 11:08:12 +0000
@@ -13,7 +13,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_V_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_V_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_V_hepmc.f	2014-03-30 11:08:12 +0000
@@ -21,7 +21,7 @@
       PARAMETER (PI=3.14159265358979312D0)
       integer j,kk,l,jpr,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_h_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_h_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_h_hepmc.f	2014-03-30 11:08:12 +0000
@@ -32,39 +32,39 @@
       l=(kk-1)*40
       call mbook(l+1,'Higgs pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+2,'Higgs pT '//weights_info(kk)
+      call mbook(l+2,'Higgs pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+3,'Higgs log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+4,'Higgs pT,|y_H|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+5,'Higgs pT,|y_H|<2 '//weights_info(kk)
+      call mbook(l+5,'Higgs pT 1,|y_H|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+6,'Higgs log[pT],|y_H|<2 '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
 
       call mbook(l+7,'j1 pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+8,'j1 pT '//weights_info(kk)
+      call mbook(l+8,'j1 pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+9,'j1 log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+10,'j1 pT,|y_j1|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+11,'j1 pT,|y_j1|<2 '//weights_info(kk)
+      call mbook(l+11,'j1 pT 1,|y_j1|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+12,'j1 log[pT],|y_j1|<2 '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
 
       call mbook(l+13,'Inc j pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+14,'Inc j pT '//weights_info(kk)
+      call mbook(l+14,'Inc j pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+15,'Inc j log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+16,'Inc j pT,|y_Ij|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,2.d2)
-      call mbook(l+17,'Inc j pT,|y_Ij|<2 '//weights_info(kk)
+      call mbook(l+17,'Inc j pT 1,|y_Ij|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,5.d2)
       call mbook(l+18,'Inc j log[pT],|y_Ij|<2'//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_hjj_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_hjj_hepmc.f	2014-02-18 17:20:34 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_hjj_hepmc.f	2014-03-30 11:08:12 +0000
@@ -47,7 +47,7 @@
 
          call mbook(l+  2,'Higgs pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+  3,'Higgs pT '
+         call mbook(l+  3,'Higgs pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+  4,'Higgs logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -58,7 +58,7 @@
 
          call mbook(l+  7,'j1 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+  8,'j1 pT '
+         call mbook(l+  8,'j1 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+  9,'j1 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -69,7 +69,7 @@
 
          call mbook(l+ 12,'j2 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 13,'j2 pT '
+         call mbook(l+ 13,'j2 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 14,'j2 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -80,7 +80,7 @@
 
          call mbook(l+ 17,'j3 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 18,'j3 pT '
+         call mbook(l+ 18,'j3 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 19,'j3 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -91,7 +91,7 @@
 
          call mbook(l+ 22,'H+j1 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 23,'H+j1 pT '
+         call mbook(l+ 23,'H+j1 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 24,'H+j1 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -102,7 +102,7 @@
 
          call mbook(l+ 27,'j1+j2 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 28,'j1+j2 pT '
+         call mbook(l+ 28,'j1+j2 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 29,'j1+j2 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -113,7 +113,7 @@
 
          call mbook(l+ 32,'syst pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 33,'syst pT '
+         call mbook(l+ 33,'syst pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 34,'syst logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -156,7 +156,7 @@
          call mbook(l+ 50,'jveto pT '
      &       //cc(i)//weights_info(kk),
      1    ((vetomax-vetomin)/dble(nbinveto)),vetomin,vetomax)
-         call mbook(l+ 51,'jveto pT '
+         call mbook(l+ 51,'jveto pT 1 '
      &       //cc(i)//weights_info(kk),
      1    ((2d0*vetomax-vetomin)/dble(nbinveto)),
      1    vetomin,2d0*vetomax)

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_lplm_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_lplm_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_lplm_hepmc.f	2014-03-30 11:08:12 +0000
@@ -42,7 +42,7 @@
       l=(kk-1)*42+(j-1)*21
       call mbook(l+ 1,'V pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+ 2,'V pt      '//cc(j)//weights_info(kk)
+      call mbook(l+ 2,'V pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+ 3,'V log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)
@@ -55,7 +55,7 @@
 c
       call mbook(l+ 7,'lm pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+ 8,'lm pt      '//cc(j)//weights_info(kk)
+      call mbook(l+ 8,'lm pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+ 9,'lm log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)
@@ -63,7 +63,7 @@
      &     ,0.2d0,-9.d0,9.d0)
       call mbook(l+11,'lp pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+12,'lp pt      '//cc(j)//weights_info(kk)
+      call mbook(l+12,'lp pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+13,'lp log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_lvl_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_lvl_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_lvl_hepmc.f	2014-03-30 11:08:12 +0000
@@ -21,7 +21,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_taptam_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_taptam_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_taptam_hepmc.f	2014-03-30 11:08:12 +0000
@@ -16,7 +16,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_tj_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_tj_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_tj_hepmc.f	2014-03-30 11:08:12 +0000
@@ -20,7 +20,7 @@
       PARAMETER (PI=3.14159265358979312D0)
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
@@ -35,7 +35,7 @@
       enddo
       nwgt=nnn+1
       nwgt_analysis=nwgt
-      do j=1,2
+      do j=1,1
       do kk=1,nwgt_analysis
       l=(kk-1)*48+(j-1)*24
       call mbook(l+ 1,'t pt        '//cc(j)//weights_info(kk)
@@ -315,7 +315,7 @@
       enddo
       nbjet=count_bj
 c fill the histograms
-      do i=1,2
+      do i=1,1
          do kk=1,nwgt_analysis
             l=(kk-1)*48+(i-1)*24
             call mfill(l+1,pttop,www(kk))

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_ttx_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_ttx_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_ttx_hepmc.f	2014-03-30 11:08:12 +0000
@@ -60,9 +60,9 @@
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
         call mbook(l+12,'delta y          '
      &       //cc(i)//weights_info(kk),0.2d0,-4.d0,4.d0)
-        call mbook(l+13,'tt azimt         '
+        call mbook(l+13,'tt azimt 1       '
      &       //cc(i)//weights_info(kk),pi/60.d0,2*pi/3,pi)
-        call mbook(l+14,'tt del R         '
+        call mbook(l+14,'tt del R 1       '
      &       //cc(i)//weights_info(kk),pi/60.d0,2*pi/3,4*pi/3)
         call mbook(l+15,'y_tb             '
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
@@ -70,11 +70,11 @@
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
         call mbook(l+17,'tt log[pi-azimt] '
      &       //cc(i)//weights_info(kk),0.05d0,-4.d0,0.1d0)
-        call mbook(l+18,'tt pt            '
+        call mbook(l+18,'tt pt 1          '
      &       //cc(i)//weights_info(kk),20.d0,80.d0,2000.d0)
-        call mbook(l+19,'tb pt            '
+        call mbook(l+19,'tb pt 1          '
      &       //cc(i)//weights_info(kk),20.d0,400.d0,2400.d0)
-        call mbook(l+20,'t pt             '
+        call mbook(l+20,'t pt 1           '
      &       //cc(i)//weights_info(kk),20.d0,400.d0,2400.d0)
       enddo
       enddo

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_ttx_v2_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_ttx_v2_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_ttx_v2_hepmc.f	2014-03-30 11:08:12 +0000
@@ -18,7 +18,7 @@
       PARAMETER (PI=3.14159265358979312D0)
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_rates_hepmc.f'
--- Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_rates_hepmc.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_rates_hepmc.f	2014-03-30 11:08:12 +0000
@@ -13,7 +13,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i,nnn
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_V.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_V.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_V.f	2014-03-30 11:08:12 +0000
@@ -23,7 +23,7 @@
       PARAMETER (PI=3.14159265358979312D0)
       integer j,kk,l,jpr
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_h.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_h.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_h.f	2014-03-30 11:08:12 +0000
@@ -26,39 +26,39 @@
       l=(kk-1)*40
       call mbook(l+1,'Higgs pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+2,'Higgs pT '//weights_info(kk)
+      call mbook(l+2,'Higgs pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+3,'Higgs log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+4,'Higgs pT,|y_H|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+5,'Higgs pT,|y_H|<2 '//weights_info(kk)
+      call mbook(l+5,'Higgs pT 1,|y_H|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+6,'Higgs log[pT],|y_H|<2 '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
 
       call mbook(l+7,'j1 pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+8,'j1 pT '//weights_info(kk)
+      call mbook(l+8,'j1 pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+9,'j1 log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+10,'j1 pT,|y_j1|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+11,'j1 pT,|y_j1|<2 '//weights_info(kk)
+      call mbook(l+11,'j1 pT 1,|y_j1|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+12,'j1 log[pT],|y_j1|<2 '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
 
       call mbook(l+13,'Inc j pT '//weights_info(kk)
      $     ,2.d0,0.d0,200.d0)
-      call mbook(l+14,'Inc j pT '//weights_info(kk)
+      call mbook(l+14,'Inc j pT 1 '//weights_info(kk)
      $     ,5.d0,0.d0,500.d0)
       call mbook(l+15,'Inc j log[pT] '//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)
       call mbook(l+16,'Inc j pT,|y_Ij|<2 '//weights_info(kk)
      $     ,2.d0,0.d0,2.d2)
-      call mbook(l+17,'Inc j pT,|y_Ij|<2 '//weights_info(kk)
+      call mbook(l+17,'Inc j pT 1,|y_Ij|<2 '//weights_info(kk)
      $     ,5.d0,0.d0,5.d2)
       call mbook(l+18,'Inc j log[pT],|y_Ij|<2'//weights_info(kk)
      $     ,0.05d0,0.1d0,5.d0)

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_hjj.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_hjj.f	2014-02-18 17:20:34 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_hjj.f	2014-03-30 11:08:12 +0000
@@ -41,7 +41,7 @@
 
          call mbook(l+  2,'Higgs pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+  3,'Higgs pT '
+         call mbook(l+  3,'Higgs pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+  4,'Higgs logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -52,7 +52,7 @@
 
          call mbook(l+  7,'j1 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+  8,'j1 pT '
+         call mbook(l+  8,'j1 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+  9,'j1 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -63,7 +63,7 @@
 
          call mbook(l+ 12,'j2 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 13,'j2 pT '
+         call mbook(l+ 13,'j2 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 14,'j2 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -74,7 +74,7 @@
 
          call mbook(l+ 17,'j3 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 18,'j3 pT '
+         call mbook(l+ 18,'j3 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 19,'j3 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -85,7 +85,7 @@
 
          call mbook(l+ 22,'H+j1 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 23,'H+j1 pT '
+         call mbook(l+ 23,'H+j1 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 24,'H+j1 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -96,7 +96,7 @@
 
          call mbook(l+ 27,'j1+j2 pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 28,'j1+j2 pT '
+         call mbook(l+ 28,'j1+j2 pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 29,'j1+j2 logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -107,7 +107,7 @@
 
          call mbook(l+ 32,'syst pT '
      &       //cc(i)//weights_info(kk),8.0d0,0.d0,400.d0)
-         call mbook(l+ 33,'syst pT '
+         call mbook(l+ 33,'syst pT 1 '
      &       //cc(i)//weights_info(kk),16.0d0,0.d0,800.d0)
          call mbook(l+ 34,'syst logpT '
      &       //cc(i)//weights_info(kk),0.08d0,0.d0,4.d0)
@@ -150,7 +150,7 @@
          call mbook(l+ 50,'jveto pT '
      &       //cc(i)//weights_info(kk),
      1    ((vetomax-vetomin)/dble(nbinveto)),vetomin,vetomax)
-         call mbook(l+ 51,'jveto pT '
+         call mbook(l+ 51,'jveto pT 1 '
      &       //cc(i)//weights_info(kk),
      1    ((2d0*vetomax-vetomin)/dble(nbinveto)),
      1    vetomin,2d0*vetomax)

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_lplm.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_lplm.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_lplm.f	2014-03-30 11:08:12 +0000
@@ -36,7 +36,7 @@
       l=(kk-1)*42+(j-1)*21
       call mbook(l+ 1,'V pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+ 2,'V pt      '//cc(j)//weights_info(kk)
+      call mbook(l+ 2,'V pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+ 3,'V log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)
@@ -49,7 +49,7 @@
 c
       call mbook(l+ 7,'lm pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+ 8,'lm pt      '//cc(j)//weights_info(kk)
+      call mbook(l+ 8,'lm pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+ 9,'lm log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)
@@ -57,7 +57,7 @@
      &     ,0.2d0,-9.d0,9.d0)
       call mbook(l+11,'lp pt      '//cc(j)//weights_info(kk)
      &     ,2.d0,0.d0,200.d0)
-      call mbook(l+12,'lp pt      '//cc(j)//weights_info(kk)
+      call mbook(l+12,'lp pt 1    '//cc(j)//weights_info(kk)
      &     ,10.d0,0.d0,1000.d0)
       call mbook(l+13,'lp log[pt] '//cc(j)//weights_info(kk)
      &     ,0.05d0,0.1d0,5.d0)

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_lvl.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_lvl.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_lvl.f	2014-03-30 11:08:12 +0000
@@ -21,7 +21,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_taptam.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_taptam.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_taptam.f	2014-03-30 11:08:12 +0000
@@ -16,7 +16,7 @@
       include 'reweight0.inc'
       integer j,kk,l,i
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_tj.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_tj.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_tj.f	2014-03-30 11:08:12 +0000
@@ -17,7 +17,7 @@
       include 'reweight0.inc'
       integer j,kk,l
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis
@@ -26,7 +26,7 @@
       common/cwgtsinfo/weights_info
       call inihist
       nwgt_analysis=nwgt
-      do j=1,2
+      do j=1,1
       do kk=1,nwgt_analysis
       l=(kk-1)*48+(j-1)*24
       call mbook(l+ 1,'t pt        '//cc(j)//weights_info(kk)
@@ -336,7 +336,7 @@
       enddo
       nbjet=count_bj
 c fill the histograms
-      do i=1,2
+      do i=1,1
          do kk=1,nwgt_analysis
             l=(kk-1)*48+(i-1)*24
             call mfill(l+1,pttop,www(kk))

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_ttx.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_ttx.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_ttx.f	2014-03-30 11:08:12 +0000
@@ -55,9 +55,9 @@
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
         call mbook(l+12,'delta y          '
      &       //cc(i)//weights_info(kk),0.2d0,-4.d0,4.d0)
-        call mbook(l+13,'tt azimt         '
+        call mbook(l+13,'tt azimt 1       '
      &       //cc(i)//weights_info(kk),pi/60.d0,2*pi/3,pi)
-        call mbook(l+14,'tt del R         '
+        call mbook(l+14,'tt del R 1       '
      &       //cc(i)//weights_info(kk),pi/60.d0,2*pi/3,4*pi/3)
         call mbook(l+15,'y_tb             '
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
@@ -65,11 +65,11 @@
      &       //cc(i)//weights_info(kk),0.1d0,-4.d0,4.d0)
         call mbook(l+17,'tt log[pi-azimt] '
      &       //cc(i)//weights_info(kk),0.05d0,-4.d0,0.1d0)
-        call mbook(l+18,'tt pt            '
+        call mbook(l+18,'tt pt 1          '
      &       //cc(i)//weights_info(kk),20.d0,80.d0,2000.d0)
-        call mbook(l+19,'tb pt            '
+        call mbook(l+19,'tb pt 1          '
      &       //cc(i)//weights_info(kk),20.d0,400.d0,2400.d0)
-        call mbook(l+20,'t pt             '
+        call mbook(l+20,'t pt 1           '
      &       //cc(i)//weights_info(kk),20.d0,400.d0,2400.d0)
       enddo
       enddo

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_ttx_v2.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_ttx_v2.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_ttx_v2.f	2014-03-30 11:08:12 +0000
@@ -18,7 +18,7 @@
       integer j,kk,l,i
       PARAMETER (PI=3.14159265358979312D0)
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_rates.f'
--- Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_rates.f	2014-02-12 10:28:36 +0000
+++ Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_rates.f	2014-03-30 11:08:12 +0000
@@ -15,7 +15,7 @@
       include 'reweight0.inc'
       integer i,kk,l
       character*5 cc(2)
-      data cc/'     ','     '/
+      data cc/'     ','Born '/
       integer nwgt,max_weight,nwgt_analysis
       common/cnwgt/nwgt
       common/c_analysis/nwgt_analysis

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIG6.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIG6.Script	2013-12-13 06:40:45 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIG6.Script	2014-03-30 11:08:12 +0000
@@ -274,6 +274,9 @@
 CMASSHW=$CMASS
 if (( $(bc <<< "$B_MASS <  0.0") )); then BMASSHW=$BMASS; fi
 if (( $(bc <<< "$B_MASS >= 0.0") )); then BMASSHW=$B_MASS; fi
+EMASSHW=$EMASS
+MUMASSHW=$MUMASS
+TAUMASSHW=$TAUMASS
 GMASSHW=$GMASS
 
 if [ $BEAM1 == 1 ]; then PART1="P"; elif [ $BEAM1 == -1 ]; then PART1="PBAR";
@@ -339,6 +342,7 @@
   $TMASS $TWIDTH                  ! M_top, Gamma_top
   $HGGMASS $HGGWIDTH              ! M_H, Ga_H
   $UMASSHW $DMASSHW $SMASSHW $CMASSHW $BMASSHW $GMASSHW ! quark and gluon masses
+  $EMASSHW $MUMASSHW $TAUMASSHW   ! Lepton masses
   $MAXPR_HW                       ! Number of events to be printed out
   $ERR_FR_HW                      ! Tolerated fraction of errors
   $B_STABLE_HW                    ! .TRUE. if it is single-top production

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIGPP.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIGPP.Script	2013-12-13 06:40:45 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIGPP.Script	2014-03-30 11:08:12 +0000
@@ -321,6 +321,9 @@
 CMASSHW=$CMASS
 if (( $(bc <<< "$B_MASS <  0.0") )); then BMASSHW=$BMASS; fi
 if (( $(bc <<< "$B_MASS >= 0.0") )); then BMASSHW=$B_MASS; fi
+EMASSHW=$EMASS
+MUMASSHW=$MUMASS
+TAUMASSHW=$TAUMASS
 GMASSHW=$GMASS
 
 if [ $BEAM1 == 1 ]; then PART1="P"; elif [ $BEAM1 == -1 ]; then PART1="PBAR";
@@ -456,21 +459,27 @@
 
 # TODO set the correct incoming particles
 
-# set the masses and widths
+# set masses and widths
+set /Herwig/Particles/e-:NominalMass $EMASSHW
+set /Herwig/Particles/e+:NominalMass $EMASSHW
+set /Herwig/Particles/mu-:NominalMass $MUMASSHW
+set /Herwig/Particles/mu+:NominalMass $MUMASSHW
+set /Herwig/Particles/tau-:NominalMass $TAUMASSHW
+set /Herwig/Particles/tau+:NominalMass $TAUMASSHW
 set /Herwig/Particles/W+:NominalMass $WMASS
+set /Herwig/Particles/W-:NominalMass $WMASS
+set /Herwig/Particles/Z0:NominalMass $ZMASS
+set /Herwig/Particles/t:NominalMass $TMASS
+set /Herwig/Particles/tbar:NominalMass $TMASS
+set /Herwig/Particles/h0:NominalMass $HGGMASS
+set /Herwig/Particles/b:NominalMass $BMASSHW
+set /Herwig/Particles/bbar:NominalMass $BMASSHW
 set /Herwig/Particles/W+:Width $WWIDTH
-set /Herwig/Particles/W-:NominalMass $WMASS
 set /Herwig/Particles/W-:Width $WWIDTH
-set /Herwig/Particles/Z0:NominalMass $ZMASS
 set /Herwig/Particles/Z0:Width $ZWIDTH
-set /Herwig/Particles/t:NominalMass $TMASS
 set /Herwig/Particles/t:Width $TWIDTH
-set /Herwig/Particles/tbar:NominalMass $TMASS
 set /Herwig/Particles/tbar:Width $TWIDTH
-set /Herwig/Particles/h0:NominalMass $HGGMASS
 set /Herwig/Particles/h0:Width $HGGWIDTH
-#set /Herwig/Particles/b:NominalMass $BMASS
-#set /Herwig/Particles/bbar:NominalMass $BMASS
 # --------------------------------------------------------------------- #
 # 3.) ADDITIONAL SETTINGS
 

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6PT.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6PT.Script	2013-12-13 06:40:45 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6PT.Script	2014-03-30 11:08:12 +0000
@@ -276,6 +276,9 @@
 CMASSPY=$CMASS
 if (( $(bc <<< "$B_MASS <  0.0") )); then BMASSPY=$BMASS; fi
 if (( $(bc <<< "$B_MASS >= 0.0") )); then BMASSPY=$B_MASS; fi
+EMASSPY=$EMASS
+MUMASSPY=$MUMASS
+TAUMASSPY=$TAUMASS
 GMASSPY=$GMASS
 
 if [ $BEAM1 == 1 ]; then PART1="P"; elif [ $BEAM1 == -1 ]; then PART1="PBAR";
@@ -355,6 +358,7 @@
   $TMASS $TWIDTH                  ! M_top, Gamma_top
   $HGGMASS $HGGWIDTH              ! M_H, Ga_H
   $UMASSPY $DMASSPY $SMASSPY $CMASSPY $BMASSPY $GMASSPY ! quark and gluon masses
+  $EMASSPY $MUMASSPY $TAUMASSPY   ! Lepton masses
   $RNDEVSEED_PY                   ! Random seed (0 for default)
   $IS_4L_PY                       ! .TRUE. if it is four-lepton production
 EOF

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6Q.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6Q.Script	2013-12-13 06:40:45 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA6Q.Script	2014-03-30 11:08:12 +0000
@@ -276,6 +276,9 @@
 CMASSPY=$CMASS
 if (( $(bc <<< "$B_MASS <  0.0") )); then BMASSPY=$BMASS; fi
 if (( $(bc <<< "$B_MASS >= 0.0") )); then BMASSPY=$B_MASS; fi
+EMASSPY=$EMASS
+MUMASSPY=$MUMASS
+TAUMASSPY=$TAUMASS
 GMASSPY=$GMASS
 
 if [ $BEAM1 == 1 ]; then PART1="P"; elif [ $BEAM1 == -1 ]; then PART1="PBAR";
@@ -355,6 +358,7 @@
   $TMASS $TWIDTH                  ! M_top, Gamma_top
   $HGGMASS $HGGWIDTH              ! M_H, Ga_H
   $UMASSPY $DMASSPY $SMASSPY $CMASSPY $BMASSPY $GMASSPY ! quark and gluon masses
+  $EMASSPY $MUMASSPY $TAUMASSPY   ! Lepton masses
   $RNDEVSEED_PY                   ! Random seed (0 for default)
   $IS_4L_PY                       ! .TRUE. if it is four-lepton production
 EOF

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script	2014-02-07 10:16:08 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script	2014-03-30 11:08:12 +0000
@@ -271,6 +271,9 @@
 CMASSPY=$CMASS
 if (( $(bc <<< "$B_MASS <  0.0") )); then BMASSPY=$BMASS; fi
 if (( $(bc <<< "$B_MASS >= 0.0") )); then BMASSPY=$B_MASS; fi
+EMASSPY=$EMASS
+MUMASSPY=$MUMASS
+TAUMASSPY=$TAUMASS
 GMASSPY=$GMASS
 
 if [ $BEAM1 == 1 ]; then iA=2212; elif [ $BEAM1 == -1 ]; then iA=-2212;
@@ -433,9 +436,9 @@
 ! 5) Final-state shower.
 TimeShower:pTmaxMatch = 1          ! Use scalup (re-check)
 TimeShower:pTmaxFudge = 1.         ! Factor changing the max scale
-#TimeShower:alphaSvalue = 0.1383   ! Alpha_s(MZ) in final-state shower
-#TimeShower:alphaSorder = 2        ! Alpha_s running order in final-state shower
-#TimeShower:alphaorder = 0         ! Alpha_EM running order in final-state shower
+TimeShower:alphaSvalue = 0.118     ! Alpha_s(MZ) in final-state shower
+TimeShower:alphaSorder = 1         ! Alpha_s running order in final-state shower
+TimeShower:alphaEMorder = 0        ! Alpha_EM running order in final-state shower
 TimeShower:interleave = on         ! If on, FSR interleaved with ISR
 TimeShower:allowBeamRecoil = on    ! If off, no energy transfer to from ISR to FSR
 TimeShower:dampenBeamRecoil = off  ! Dampens the effect of beam recoil
@@ -458,9 +461,9 @@
 ! 6) Initial-state shower.
 SpaceShower:pTmaxMatch = 1         ! Use scalup (re-check)
 SpaceShower:pTmaxFudge = 1.        ! Factor changing the max scale
-#SpaceShower:alphaSvalue = 0.137   ! Alpha_s(MZ) in final-state shower
-#SpaceShower:alphaSorder = 2       ! Alpha_s running order in final-state shower
-#SpaceShower:alphaEMorder = 0      ! Alpha_EM running order in final-state shower
+SpaceShower:alphaSvalue = 0.118    ! Alpha_s(MZ) in final-state shower
+SpaceShower:alphaSorder = 1        ! Alpha_s running order in final-state shower
+SpaceShower:alphaEMorder = 0       ! Alpha_EM running order in final-state shower
 SpaceShower:QCDshower = on         ! QCD final-state shower
 SpaceShower:QEDshowerByQ = off     ! Prevent quarks from radiating photons
 SpaceShower:QEDshowerByL = off     ! Prevent leptons from radiating photons
@@ -481,6 +484,9 @@
 4:m0 = $CMASSPY                    ! charm mass
 5:m0 = $BMASSPY                    ! bottom mass
 6:m0 = $TMASS                      ! top mass
+11:m0 = $EMASSPY                   ! electron mass
+13:m0 = $MUMASSPY                  ! muon mass
+15:m0 = $TAUMASSPY                 ! tauon mass
 23:m0 = $ZMASS                     ! Z mass
 24:m0 = $WMASS                     ! W mass
 25:m0 = $HGGMASS                   ! Higgs mass

=== modified file 'Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f'
--- Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2013-12-13 13:30:05 +0000
+++ Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2014-03-30 11:08:12 +0000
@@ -13,7 +13,7 @@
       CHARACTER *7 NORM_EVENT
       INTEGER MQQ
       COMMON/cMQQ/MQQ
-      REAL*8 TMPLAM,GAMT0,ERR_FR
+      REAL*8 TMPLAM,GAMT0,ERR_FR,MME,MMM,MMT
       INTEGER IPDF
       CHARACTER * 70 LHAPDF
       LOGICAL LHACRTL,OLDFORM,PI_STABLE,WP_STABLE,WM_STABLE,Z_STABLE,
@@ -106,6 +106,12 @@
       DO I=1,5
          RMASS(I+6)=RMASS(I)
       ENDDO
+      WRITE(*,*)'Enter lepton (e,mu,tau) masses: should be zero'
+      READ(*,*)mme,mmm,mmt
+      if(mme.ne.0d0.or.mmm.ne.0d0.or.mmt.ne.0d0)then
+         write(*,*)'nonzero lepton masses'
+         stop
+      endif
 C Set electron and muon masses equal to zero to avoid rounding problems
       RMASS(121)=0.D0
       RMASS(123)=0.D0

=== modified file 'Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f'
--- Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f	2013-12-13 06:40:45 +0000
+++ Template/NLO/MCatNLO/srcPythia/madfks_pydriver.f	2014-03-30 11:08:12 +0000
@@ -14,10 +14,10 @@
       CHARACTER *100 TMPSTR
       INTEGER N,NPAD,K,MSTU,MSTJ,KCHG,MDCY,MDME,KFDP,MSEL,MSELPD,PYCOMP,
      &MSUB,KFIN,MSTP,MSTI,MODPDF(2),MAXEV,MAXPR,ID,NGENPD,NGEN
-      DOUBLE PRECISION M_Z,M_W,M_T,M_H,M_D,M_U,M_S,M_C,M_B,M_G,
+      DOUBLE PRECISION M_Z,M_W,M_T,M_H,M_D,M_U,M_S,M_C,M_B,M_G,M_E,M_MU,M_TAU,
      &G_Z,G_W,G_T,G_H,QCDLAM,PBEAM1,PBEAM2,EMMIN,EMMAX,GAMMAX,XSEC
       DOUBLE PRECISION P,V,PARJ,PMAS,PARF,VCKM,BRAT,CKIN,PARP,PARI,PARU
-      COMMON/CMASSES/M_Z,M_W,M_T,M_H,M_D,M_U,M_S,M_C,M_B,M_G
+      COMMON/CMASSES/M_Z,M_W,M_T,M_H,M_D,M_U,M_S,M_C,M_B,M_G,M_E,M_MU,M_TAU
       COMMON/CWIDTHS/G_Z,G_W,G_T,G_H
 C QQIN IS THE EVENT FILE
       CHARACTER*50 QQIN,PART1,PART2
@@ -66,6 +66,9 @@
       M_C=0D0
       M_B=0D0
       M_G=0D0
+      M_E=0D0
+      M_MU=0D0
+      M_TAU=0D0
       G_Z=0D0
       G_W=0D0
       G_T=0D0
@@ -225,6 +228,8 @@
       READ(*,*)M_H,G_H
       WRITE(*,*)'Enter quark (d,u,s,c,b) and gluon masses'
       READ(*,*)M_D,M_U,M_S,M_C,M_B,M_G
+      WRITE(*,*)'Enter lepton (e,mu,tau) masses'
+      READ(*,*)M_E,M_MU,M_TAU
 
 C---OPEN LH FILE
       OPEN(UNIT=MSTP(161),FILE=QQIN)
@@ -434,9 +439,9 @@
 C THESE QUANTITIES IS HEP-PH/0603175 AND IN PARENTHESES I INDICATE THE PAGE
       IMPLICIT NONE
       INTEGER ISHW_PY,JPR
-      DOUBLE PRECISION M_Z,M_W,M_T,M_H,M_D,M_U,M_S,M_C,M_B,M_G,
+      DOUBLE PRECISION M_Z,M_W,M_T,M_H,M_D,M_U,M_S,M_C,M_B,M_G,M_E,M_MU,M_TAU,
      &G_Z,G_W,G_T,G_H,QCDLAM
-      COMMON/CMASSES/M_Z,M_W,M_T,M_H,M_D,M_U,M_S,M_C,M_B,M_G
+      COMMON/CMASSES/M_Z,M_W,M_T,M_H,M_D,M_U,M_S,M_C,M_B,M_G,M_E,M_MU,M_TAU
       COMMON/CWIDTHS/G_Z,G_W,G_T,G_H
       COMMON/CQCD/QCDLAM
       COMMON/CISHW/ISHW_PY
@@ -460,6 +465,9 @@
       PMAS(4,1)=M_C
       PMAS(5,1)=M_B
       PMAS(6,1)=M_T
+      PMAS(11,1)=M_E
+      PMAS(13,1)=M_MU
+      PMAS(15,1)=M_TAU
       PMAS(21,1)=M_G
 C
       PMAS(6,2)=G_T

=== modified file 'Template/NLO/Source/cuts.inc'
--- Template/NLO/Source/cuts.inc	2013-11-21 14:36:04 +0000
+++ Template/NLO/Source/cuts.inc	2014-03-30 11:08:12 +0000
@@ -1,5 +1,5 @@
 C
-C     NEW CUTS
+C     DEFAULT CUTS FROM RUN_CARD>DAt
 C
       REAL*8 JETALGO,JETRADIUS,PTJ,ETAJ  ! jet cuts
       REAL*8 PTL,ETAL,DRLL,MLL           ! lepton cuts
@@ -10,74 +10,10 @@
       LOGICAL ISOEM                           ! photon cuts
       COMMON /TO_ISOGAMMA_CUTS/PTGMIN,R0GAMMA,XN,EPSGAMMA,ETAGAMMA
       COMMON /TO_ISOGAMMA_EM/ISOEM
-
-C
-C     AUTOMATIC CUTS
-C
-c      REAL*8 PTJ,PTB,PTA,PTL,PTONIUM  ! ptmin
-      REAL*8 PTB,PTA,PTONIUM  ! ptmin
-      REAL*8 EJ,EB,EA,EL,ETAONIUM      ! min energy
-c      REAL*8 ETAJ,ETAB,ETAA,ETAL ! min rapidity
-      REAL*8 ETAB,ETAA ! min rapidity
-c      REAL*8 DRJJ,DRBB,DRLL,DRAA,DRBJ,DRAJ,DRJL,DRAB,DRBL,DRAL ! min distance
-      REAL*8 DRJJ,DRBB,DRAA,DRBJ,DRAJ,DRJL,DRAB,DRBL,DRAL ! min distance
-      REAL*8 MMJJ,MMLL,MMAA,MMBB !min inv mass
-      COMMON /TO_AUTO_CUTS/
-c     &     PTJ ,PTB ,PTA ,PTL ,PTONIUM, ! pt
-     &     PTB ,PTA ,PTONIUM, ! pt
-     &     EJ  ,EB  ,EA  ,EL  , ! energy
-c     &     ETAJ,ETAB,ETAA,ETAL,ETAONIUM, ! rapidity 
-     &     ETAB,ETAA,ETAONIUM, ! rapidity 
-     &     DRJJ,DRAJ,DRBJ,DRJL, ! delta R
-     &     DRBB,DRAB,DRBL,
-     &     DRAA,DRAL,
-c     &     DRLL,
-     &     MMJJ,MMLL,MMAA,MMBB ! invariant masses of opposite sign
-
-C     
-C     MIN/MAX cuts
-C
-      REAL*8 PTJmax ,PTBmax ,PTAmax ,PTLmax    ! max pt
-      REAL*8 EJmax  ,EBmax  ,EAmax  ,ELmax     ! max rap
-      REAL*8 ETAJmin,ETABmin,ETAAmin,ETALmin   ! min rap
-      REAL*8 DRJJmax,DRBBmax,DRLLmax,DRAAmax
-      REAL*8 DRBJmax,DRAJmax,DRJLmax,DRABmax,DRBLmax,DRALmax
-      REAL*8 MMJJmax,MMLLmax,MMAAmax,MMBBmax !max inv mass 
-      REAL*8 MMNL,MMNLMAX ! invariant mass of all leptons
-      REAL*8 cutuse
-      REAL*8 ptj1min,ptj2min,ptj3min,ptj4min
-      REAL*8 ptj1max,ptj2max,ptj3max,ptj4max
-      REAL*8 ht2min,ht3min,ht4min
-      REAL*8 ht2max,ht3max,ht4max
-      REAL*8 htjmin,htjmax
-      REAL*8 misset,missetmax,ptheavy
       integer maxjetflavor
-	
-      COMMON /to_min_max_cuts/
-     &     PTJmax,PTBmax,PTAmax,PTLmax,
-     &     EJmax,EBmax,EAmax,ELmax,
-     &     ETAJmin,ETABmin,ETAAmin,ETALmin,
-     &     DRJJmax,DRBBmax,DRLLmax,DRAAmax,
-     &     DRBJmax,DRAJmax,DRJLmax,DRABmax,DRBLmax,DRALmax,
-     &     MMJJmax,MMLLmax,MMAAmax,MMBBmax,
-     &     MMNL,MMNLmax,
-     &     ptj1min,ptj1max,
-     &     ptj2min,ptj2max,
-     &     ptj3min,ptj3max,
-     &     ptj4min,ptj4max,
-     &     cutuse,
-     &     ht2min,ht3min,ht4min,
-     &     ht2max,ht3max,ht4max,
-     &     htjmin,htjmax,
-     &     misset,missetmax,ptheavy,
-     &     maxjetflavor
-
-C
-C     SPECIAL CUTS
-C
-C-- this are somewhat old kind of cuts that where used for specific
-c   purposes
-c
-      real*8 xptj,xptb,xpta,xptl,xmtc
-      real*8 xetamin,xqcut,deltaeta
-      common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
+      COMMON /to_min_max_cuts/maxjetflavor
+C
+C     SPECIAL CUTS (FOR MERGING)
+C
+      real*8 xmtc,xqcut
+      common /to_specxpt/xmtc,xqcut

=== modified file 'Template/NLO/Source/makefile'
--- Template/NLO/Source/makefile	2013-03-18 13:23:52 +0000
+++ Template/NLO/Source/makefile	2014-03-30 11:08:12 +0000
@@ -1,8 +1,11 @@
 include make_opts
 
 LIBDIR= ../lib/
+MCATNLOLIBDIR= ../MCatNLO/lib/
 BINDIR= ../bin/
 PDFDIR= ./PDF/
+CUTTOOLSDIR= ./CutTools/
+STDHEPDIR= ./StdHEP/
 
 PROCESS= hfill.o matrix.o myamp.o 
 
@@ -75,15 +78,37 @@
 	$(FC) $(FFLAGS) -o combine_events $(COMBINE) $(LIBDIR)/libmodel.a $(LIBDIR)/libpdf.a 
 	mv combine_events $(BINDIR)
 
+CutTools: $(LIBDIR)libcts.a
+
+$(LIBDIR)libcts.a: $(CUTTOOLSDIR)
+	cd $(CUTTOOLSDIR); make
+	ln -sf ../Source/CutTools/includects/libcts.a $(LIBDIR)libcts.a
+	ln -sf ../Source/CutTools/includects/mpmodule.mod $(LIBDIR)mpmodule.mod
+
+StdHEP: $(MCATNLOLIBDIR)libstdhep.a  
+
+$(MCATNLOLIBDIR)libstdhep.a: $(STDHEPDIR)
+	cd $(STDHEPDIR); make all
+	ln -sf ../../Source/StdHEP/lib/libstdhep.a $(MCATNLOLIBDIR)libstdhep.a
+	ln -sf ../../Source/StdHEP/lib/libFmcfio.a $(MCATNLOLIBDIR)libFmcfio.a
+
 clean4pdf:
 	rm -f ../lib/libpdf.a
 	rm -f ../lib/libgeneric.a
 
+cleanCT:
+	cd $(CUTTOOLSDIR); make clean; cd ..
+
+cleanStdHEP:
+	cd $(STDHEPDIR); make clean; cd ..
+
 clean:
 	rm -f *.o
 	rm -f param_card.inc run_card.inc
-	cd MODEL; make clean
+	cd MODEL; make clean; cd ..
 	cd CERNLIB; make clean; cd ..
-	cd PDF; make clean
-	cd DHELAS; make clean
+	cd PDF; make clean; cd ..
+	cd DHELAS; make clean; cd ..
+	cd $(CUTTOOLSDIR); make clean; cd ..
+	cd $(STDHEPDIR); make clean; cd ..
 	rm -f $(BINDIR)/combine_events $(BINDIR)/gen_ximprove

=== modified file 'Template/NLO/Source/setrun.f'
--- Template/NLO/Source/setrun.f	2013-10-24 14:27:07 +0000
+++ Template/NLO/Source/setrun.f	2014-03-30 11:08:12 +0000
@@ -76,149 +76,11 @@
       call to_upper(shower_MC)
 
 
-
-
-c*********************************************************************
-c     Minimum pt's                                                   *
-c*********************************************************************
-      ptb     = 0d0
-      pta     = 0d0
-      misset  = 0d0
-      ptonium = 0d0
-      
-c*********************************************************************
-c     Maximum pt's                                                   *
-c*********************************************************************
-      ptjmax=1d5
-      ptbmax=1d5
-      ptamax=1d5
-      ptlmax=1d5
-      missetmax=1d5
-
-c*********************************************************************
-c     Maximum rapidity (absolute value)                              *
-c*********************************************************************
-      etab=1d2
-      etaa=1d2
-      etaonium=1d2
-      etajmin=0d0
-      etabmin=0d0
-      etaamin=0d0
-      etalmin=0d0
-
-      ej=0d0
-      eb=0d0
-      ea=0d0
-      el=0d0
-
-c*********************************************************************
-c     Maximum E's                                                    *
-c*********************************************************************
-      ejmax=1d5
-      ebmax=1d5
-      eamax=1d5
-      elmax=1d5
-
-c*********************************************************************
-c     Minimum DeltaR distance                                        *
-c*********************************************************************
-      drjj=0d0
-      drbb=0d0
-      draa=0d0
-      drbj=0d0
-      draj=0d0
-      drjl=0d0
-      drab=0d0
-      drbl=0d0
-      dral=0d0
-
-c*********************************************************************
-c     Maximum DeltaR distance                                        *
-c*********************************************************************
-      drjjmax=1d2
-      drbbmax=1d2
-      drllmax=1d2
-      draamax=1d2
-      drbjmax=1d2
-      drajmax=1d2
-      drjlmax=1d2
-      drabmax=1d2
-      drblmax=1d2
-      dralmax=1d2
-
-c*********************************************************************
-c     Minimum invariant mass for pairs                               *
-c*********************************************************************
-      mmjj=0d0
-      mmbb=0d0
-      mmaa=0d0
-      mmll=0d0
-
-c*********************************************************************
-c     Maximum invariant mass for pairs                               *
-c*********************************************************************
-      mmjjmax=1d5
-      mmbbmax=1d5
-      mmaamax=1d5
-      mmllmax=1d5
-
-c*********************************************************************
-c     Min Maxi invariant mass for all leptons                        *
-c*********************************************************************
-      mmnl=0d0
-      mmnlmax=1d5
-
-c*********************************************************************
-c     Inclusive cuts                                                 *
-c*********************************************************************
-      xptj=0d0
-      xptb=0d0
-      xpta=0d0
-      xptl=0d0
+c merging cuts
+      xqcut=0d0
       xmtc=0d0
-
-c*********************************************************************
-c     WBF cuts                                                       *
-c*********************************************************************
-      xetamin=0d0
-      deltaeta=0d0
-
-c*********************************************************************
-c     Jet measure cuts                                               *
-c*********************************************************************
-      xqcut=0d0
       d=1d0
-
-c*********************************************************************
-c Set min pt of one heavy particle                                   *
-c*********************************************************************
-        ptheavy=0d0
-
-c*********************************************************************
-c Check   the pt's of the jets sorted by pt                          *
-c*********************************************************************
-        ptj1min=0d0
-        ptj1max=1d5
-        ptj2min=0d0
-        ptj2max=1d5
-        ptj3min=0d0
-        ptj3max=1d5
-        ptj4min=0d0
-        ptj4max=1d5
-        cutuse=0d0
-
-c*********************************************************************
-c Check  Ht                                                          *
-c*********************************************************************
-        ht2min=0d0
-        ht3min=0d0
-        ht4min=0d0
-        ht2max=1d5
-        ht3max=1d5
-        ht4max=1d5
-        htjmin=0d0
-        htjmax=1d5
-
+      
 c*********************************************************************
 c     Random Number Seed                                             *
 c*********************************************************************

=== modified file 'Template/NLO/SubProcesses/MCmasses_HERWIG6.inc'
--- Template/NLO/SubProcesses/MCmasses_HERWIG6.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/SubProcesses/MCmasses_HERWIG6.inc	2014-03-30 11:08:12 +0000
@@ -3,4 +3,10 @@
       mcmass(3)=0.50d0
       mcmass(4)=1.55d0
       mcmass(5)=4.95d0
+      mcmass(11)=0.d0
+      mcmass(12)=0.d0
+      mcmass(13)=0.d0
+      mcmass(14)=0.d0
+      mcmass(15)=0.d0
+      mcmass(16)=0.d0
       mcmass(21)=0.75d0

=== modified file 'Template/NLO/SubProcesses/MCmasses_HERWIGPP.inc'
--- Template/NLO/SubProcesses/MCmasses_HERWIGPP.inc	2012-10-29 06:09:23 +0000
+++ Template/NLO/SubProcesses/MCmasses_HERWIGPP.inc	2014-03-30 11:08:12 +0000
@@ -3,4 +3,10 @@
       mcmass(3)=0.d0
       mcmass(4)=0.d0
       mcmass(5)=0.d0
+      mcmass(11)=0.510998928d-3
+      mcmass(12)=0.d0
+      mcmass(13)=0.1056583715d0
+      mcmass(14)=0.d0
+      mcmass(15)=1.77682d0
+      mcmass(16)=0.d0
       mcmass(21)=0.d0

=== modified file 'Template/NLO/SubProcesses/MCmasses_PYTHIA6PT.inc'
--- Template/NLO/SubProcesses/MCmasses_PYTHIA6PT.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/SubProcesses/MCmasses_PYTHIA6PT.inc	2014-03-30 11:08:12 +0000
@@ -3,4 +3,10 @@
       mcmass(3)=0.50d0
       mcmass(4)=1.50d0
       mcmass(5)=4.80d0
+      mcmass(11)=0.510998928d-3
+      mcmass(12)=0.d0
+      mcmass(13)=0.1056583715d0
+      mcmass(14)=0.d0
+      mcmass(15)=1.77682d0
+      mcmass(16)=0.d0
       mcmass(21)=0.0d0

=== modified file 'Template/NLO/SubProcesses/MCmasses_PYTHIA6Q.inc'
--- Template/NLO/SubProcesses/MCmasses_PYTHIA6Q.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/SubProcesses/MCmasses_PYTHIA6Q.inc	2014-03-30 11:08:12 +0000
@@ -3,4 +3,10 @@
       mcmass(3)=0.50d0
       mcmass(4)=1.50d0
       mcmass(5)=4.80d0
+      mcmass(11)=0.510998928d-3
+      mcmass(12)=0.d0
+      mcmass(13)=0.1056583715d0
+      mcmass(14)=0.d0
+      mcmass(15)=1.77682d0
+      mcmass(16)=0.d0
       mcmass(21)=0.0d0

=== modified file 'Template/NLO/SubProcesses/MCmasses_PYTHIA8.inc'
--- Template/NLO/SubProcesses/MCmasses_PYTHIA8.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/SubProcesses/MCmasses_PYTHIA8.inc	2014-03-30 11:08:12 +0000
@@ -3,4 +3,10 @@
       mcmass(3)=0.50d0
       mcmass(4)=1.50d0
       mcmass(5)=4.80d0
+      mcmass(11)=0.510998928d-3
+      mcmass(12)=0.d0
+      mcmass(13)=0.1056583715d0
+      mcmass(14)=0.d0
+      mcmass(15)=1.77682d0
+      mcmass(16)=0.d0
       mcmass(21)=0.0d0

=== modified file 'Template/NLO/SubProcesses/add_write_info.f'
--- Template/NLO/SubProcesses/add_write_info.f	2013-08-29 18:50:43 +0000
+++ Template/NLO/SubProcesses/add_write_info.f	2014-03-30 11:08:12 +0000
@@ -1069,7 +1069,7 @@
       common/to_mass/  emass
 
 c Monte Carlo masses: use PDG conventions
-      double precision mcmass(-5:21)
+      double precision mcmass(-16:21)
       common/cmcmass/mcmass
 
 c Masses used by write_events_lhe
@@ -1096,6 +1096,8 @@
       xmj=-1.d0
       xm1=-1.d0
       xm2=-1.d0
+c WARNING: what follows will need to be reconsidered the case of 
+c QED corrections, for what is relevant to i_fks and j_fks
       do i=1,nexternal
         if(i.eq.i_fks)then
           if(pmass(i).ne.0.d0)then
@@ -1186,12 +1188,13 @@
         else
           idpart=jpart(1,i)
           if( idpart.eq.21 .or.
-     #        (abs(idpart).ge.1.and.abs(idpart).le.5) )then
+     #        (abs(idpart).ge.1.and.abs(idpart).le.5) .or.
+     #        (abs(idpart).ge.11.and.abs(idpart).le.16) )then
             if(pmass(i).eq.0.d0)then
               tmpmass=mcmass(idpart)
             else
-c If MadFKS has a non-zero mass for a "light" quark, one probably wants
-c to use that in the shower phase as well (ie bottom or charm production).
+c If MadFKS has a non-zero mass for a "light" quark or lepton, one probably 
+c wants to use that in the shower phase as well (ie bottom or charm production).
 c One may use equivalently a condition on maxjetflavor
               tmpmass=pmass(i)
             endif

=== modified file 'Template/NLO/SubProcesses/ajob_template'
--- Template/NLO/SubProcesses/ajob_template	2013-11-21 19:07:16 +0000
+++ Template/NLO/SubProcesses/ajob_template	2014-03-30 11:08:12 +0000
@@ -86,6 +86,7 @@
     link1up FKS_params.dat
     link1up OLE_order.olc
     link1up param_card.dat
+    link1up libLHAPDF.so
     
 #   check where is the basic file for the creation of input_app.txt
 #       

=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2014-01-09 16:50:33 +0000
+++ Template/NLO/SubProcesses/cuts.f	2014-03-30 11:08:12 +0000
@@ -1,195 +1,107 @@
-      logical function pass_point(p)
+c
+c This file contains the default cuts (as defined in the run_card.dat)
+c and can easily be extended by the user to include other.  This
+c function should return true if event passes cuts
+c (passcuts_user=.true.) and false otherwise (passcuts_user=.false.).
+c
+c NOTE THAT ONLY IRC-SAFE CUTS CAN BE APPLIED OTHERWISE THE INTEGRATION
+c MIGHT NOT CONVERGE
+c
+      logical function passcuts_user(p,istatus,ipdg)
       implicit none
-      double precision p
-      logical passcuts
-      external passcuts
-      pass_point = .true.
-c      pass_point = passcuts(p)
-      end
-
-
-      LOGICAL FUNCTION PASSCUTS(P,rwgt)
-C**************************************************************************
-C     INPUT:
-C            P(0:3,1)           MOMENTUM OF INCOMING PARTON
-C            P(0:3,2)           MOMENTUM OF INCOMING PARTON
-C            P(0:3,3)           MOMENTUM OF d
-C            P(0:3,4)           MOMENTUM OF b
-C            P(0:3,5)           MOMENTUM OF bbar
-C            P(0:3,6)           MOMENTUM OF e+
-C            P(0:3,7)           MOMENTUM OF ve
-C            COMMON/JETCUTS/   CUTS ON JETS
-C     OUTPUT:
-C            TRUE IF EVENTS PASSES ALL CUTS LISTED
-C**************************************************************************
-C
-C *WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING*
-C
-C In MadFKS, the momenta given in input to this function are in the
-C reduced parton c.m. frame. If need be, boost them to the lab frame.
-C The rapidity of this boost is
-C
-C       YBST_TIL_TOLAB
-C
-C given in the common block /PARTON_CMS_STUFF/
-C
-C This is the rapidity that enters in the arguments of the sinh() and
-C cosh() of the boost, in such a way that
-C       ylab = ycm - ybst_til_tolab
-C where ylab is the rapidity in the lab frame and ycm the rapidity
-C in the center-of-momentum frame.
-C
-C *WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING*
-c
-      IMPLICIT NONE
-c
-c     Constants
-c
-      include 'genps.inc'
-      include "nexternal.inc"
-C
-C     ARGUMENTS
-C
-      REAL*8 P(0:3,nexternal),rwgt
-
-C
-C     LOCAL
-C
-      LOGICAL FIRSTTIME
-      DATA FIRSTTIME/.TRUE./
+c This includes the 'nexternal' parameter that labels the number of
+c particles in the (n+1)-body process
+      include 'nexternal.inc'
+c This include file contains common blocks filled with the cuts defined
+c in the run_card.dat
+      include 'cuts.inc'
+c
+c This is an array which is '-1' for initial state and '1' for final
+c state particles
+      integer istatus(nexternal)
+c This is an array with (simplified) PDG codes for the particles. Note
+c that channels that are combined (i.e. they have the same matrix
+c elements) are given only 1 set of PDG codes. This means, e.g., that
+c when using a 5-flavour scheme calculation (massless b quark), no
+c b-tagging can be applied.
+      integer iPDG(nexternal)
+c The array of the momenta and masses of the initial and final state
+c particles in the lab frame. The format is "E, px, py, pz, mass", while
+c the second dimension loops over the particles in the process. Note
+c that these are the (n+1)-body particles; for the n-body there is one
+c momenta equal to all zero's (this is not necessarily the last particle
+c in the list). If one uses IR-safe obserables only, there should be no
+c difficulty in using this.
+      double precision p(0:4,nexternal)
+c
+C     external functions that can be used. Some are defined in this
+C     file, others are in ./Source/kin_functions.f
+      REAL*8 R2_04,invm2_04,pt_04,eta_04,pt,eta
+      external R2_04,invm2_04,pt_04,eta_04,pt,eta
+c local integers
       integer i,j
-C
-C     EXTERNAL
-C
-      REAL*8 R2,DOT,ET,RAP,DJ,SumDot,pt,rewgt,eta
-      logical cut_bw
-      external cut_bw,rewgt,eta,r2,dot,et,rap,dj,sumdot,pt
-C
-C     GLOBAL
-C
-      include 'run.inc'
-      include 'cuts.inc'
-c For boosts
-      double precision ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
-      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,
-     #                        sqrtshat,shat
-      double precision pjetlab(0:3,nexternal)
-      double precision chybst,shybst,chybstmo
-      double precision xd(1:3)
-      data (xd(i),i=1,3)/0,0,1/
-c Jets and charged leptons
-      LOGICAL  IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
-      LOGICAL  IS_A_PH(NEXTERNAL)
-      COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
-      include 'coupl.inc'
 c jet cluster algorithm
       integer nQCD,NJET,JET(nexternal)
-      double precision plab(0:3, nexternal)
       double precision pQCD(0:3,nexternal),PJET(0:3,nexternal)
       double precision rfj,sycut,palg,amcatnlo_fastjetdmerge
       integer njet_eta
+      integer mm
 c Photon isolation
       integer nph,nem,k,nin
       double precision ptg,chi_gamma_iso,iso_getdrv40
       double precision Etsum(0:nexternal)
       real drlist(nexternal)
       double precision pgamma(0:3,nexternal),pem(0:3,nexternal)
-     $     ,pgammalab(0:3)
       logical alliso
 c Sort array of results: ismode>0 for real, isway=0 for ascending order
       integer ismode,isway,izero,isorted(nexternal)
       parameter (ismode=1)
       parameter (isway=0)
       parameter (izero=0)
-
-
-
-      integer mm
-
-C-----
-C  BEGIN CODE
-C-----
-      PASSCUTS=.TRUE.             !EVENT IS OK UNLESS OTHERWISE CHANGED
-      IF (FIRSTTIME) THEN
-         FIRSTTIME=.FALSE.
-         write (*,*) '================================================='
-         write (*,*) 'From cuts.f'
-         if (jetalgo.eq.1) then
-            write (*,*) 'Jets are defined with the kT algorithm'
-         elseif (jetalgo.eq.0) then
-            write (*,*) 'Jets are defined with the C/A algorithm'
-         elseif (jetalgo.eq.-1) then
-            write (*,*) 'Jets are defined with the anti-kT algorithm'
-         else
-            write (*,*) 'Jet algorithm not defined in the run_card.dat,'
-     &           //'or not correctly processed by the code.',jetalgo
-         endif
-         write (*,*) 'with a mimumal pT of ',ptj,'GeV'
-         if (etaj.gt.0) then
-            write (*,*) 'and maximal pseudo-rapidity of ',etaj,'.'
-         else
-            write (*,*) 'and no maximal pseudo-rapidity.'
-         endif
-         write (*,*) 'Charged leptons are required to have at least',ptl
-     &        ,'GeV of transverse momentum and'
-         if (etal.gt.0) then
-            write (*,*) 'pseudo rapidity of maximum',etal,'.'
-         else
-            write (*,*) 'no maximum for the pseudo rapidity.'
-         endif
-         write (*,*) 'Opposite charged lepton pairs need to be'//
-     &        ' separated by at least ',drll
-         write (*,*) 'and have an invariant mass of',mll,' GeV'
-         write (*,*) '================================================='
-      ENDIF
-c
-c     Make sure have reasonable 4-momenta
-c
-      if (p(0,1) .le. 0d0) then
-         passcuts=.false.
-         return
-      endif
-
-c     Also make sure there's no INF or NAN
-      do i=1,nexternal
-         do j=0,3
-            if(p(j,i).gt.1d32.or.p(j,i).ne.p(j,i))then
-               passcuts=.false.
-               return
-            endif
-         enddo
-      enddo
-
-      rwgt=1d0
-
-c Uncomment for bypassing charged lepton cuts
-c$$$      goto 124
-
-c Boost the momenta p(0:3,nexternal) to the lab frame plab(0:3,nexternal)
-      chybst=cosh(ybst_til_tolab)
-      shybst=sinh(ybst_til_tolab)
-      chybstmo=chybst-1.d0
-      do i=1,nexternal
-         call boostwdir2(chybst,shybst,chybstmo,xd,
-     &        p(0,i),plab(0,i))
-      enddo
-
+c logicals that define if particles are leptons, jets or photons. These
+c are filled from the PDG codes (iPDG array) in this function.
+      logical is_a_lp(nexternal),is_a_lm(nexternal),is_a_j(nexternal)
+     $     ,is_a_ph(nexternal)
+
+      passcuts_user=.true. ! event is okay; otherwise it is changed
+
+C***************************************************************
+C***************************************************************
+C Cuts from the run_card.dat
+C***************************************************************
+C***************************************************************
 c
 c CHARGED LEPTON CUTS
 c
+c find the charged leptons (also used in the photon isolation cuts below)
+      do i=1,nexternal
+         if(istatus(i).eq.1 .and.
+     &    (ipdg(i).eq.11 .or. ipdg(i).eq.13 .or. ipdg(i).eq.15)) then
+            is_a_lm(i)=.true.
+         else
+            is_a_lm(i)=.false.
+         endif
+         if(istatus(i).eq.1 .and.
+     &    (ipdg(i).eq.-11 .or. ipdg(i).eq.-13 .or. ipdg(i).eq.-15)) then
+            is_a_lp(i)=.true.
+         else
+            is_a_lp(i)=.false.
+         endif
+      enddo
+c apply the charged lepton cuts
       do i=nincoming+1,nexternal
          if (is_a_lp(i).or.is_a_lm(i)) then
 c transverse momentum
             if (ptl.gt.0d0) then
-               if (pt(p(0,i)).lt.ptl) then
-                  passcuts=.false.
+               if (pt_04(p(0,i)).lt.ptl) then
+                  passcuts_user=.false.
                   return
                endif
             endif
 c pseudo-rapidity
             if (etal.gt.0d0) then
-               if (abs(eta(plab(0,i))).gt.etal) then
-                  passcuts=.false.
+               if (abs(eta_04(p(0,i))).gt.etal) then
+                  passcuts_user=.false.
                   return
                endif
             endif
@@ -198,14 +110,14 @@
                do j=nincoming+1,nexternal
                   if (is_a_lm(j)) then
                      if (drll.gt.0d0) then
-                        if (R2(plab(0,i),plab(0,j)).lt.drll**2) then
-                           passcuts=.false.
+                        if (R2_04(p(0,i),p(0,j)).lt.drll**2) then
+                           passcuts_user=.false.
                            return
                         endif
                      endif
                      if (mll.gt.0d0) then
-                        if (sumdot(p(0,i),p(0,j),1d0).lt.mll**2) then
-                           passcuts=.false.
+                        if (invm2_04(p(0,i),p(0,j),1d0).lt.mll**2) then
+                           passcuts_user=.false.
                            return
                         endif
                      endif
@@ -214,19 +126,22 @@
             endif
          endif
       enddo
-
- 124  continue
-
 c
 c JET CUTS
 c
-c Uncomment for bypassing jet algo and cuts, and photon isolation
-c$$$      goto 123
+c find the jets
+      do i=1,nexternal
+         if (istatus(i).eq.1 .and.
+     &        (abs(ipdg(i)).le.maxjetflavor .or. ipdg(i).eq.21)) then
+            is_a_j(i)=.true.
+         else
+            is_a_j(i)=.false.
+         endif
+      enddo
 
 c If we do not require a mimimum jet energy, there's no need to apply
 c jet clustering and all that.
       if (ptj.ne.0d0.or.ptgmin.ne.0d0) then
-
 c Put all (light) QCD partons in momentum array for jet clustering.
 c From the run_card.dat, maxjetflavor defines if b quark should be
 c considered here (via the logical variable 'is_a_jet').  nQCD becomes
@@ -237,16 +152,12 @@
             if (is_a_j(j)) then
                nQCD=nQCD+1
                do i=0,3
-                  pQCD(i,nQCD)=p(i,j) ! Use C.o.M. frame momenta
+                  pQCD(i,nQCD)=p(i,j)
                enddo
             endif
          enddo
-
       endif
 
-
-c Uncomment for bypassing jet algo and cuts
-c$$$      goto 122
       if (ptj.gt.0d0.and.nQCD.gt.1) then
 
 c Cut some peculiar momentum configurations, i.e. two partons very soft.
@@ -258,7 +169,7 @@
             if(abs(pQCD(0,j)/p(0,1)).lt.1.d-8) mm=mm+1
          enddo
          if(mm.gt.1)then
-            passcuts=.false.
+            passcuts_user=.false.
             return
          endif
 
@@ -279,31 +190,24 @@
 c     jet algorithm:               palg, 1.0=kt, 0.0=C/A, -1.0 = anti-kt
 c
 c     OUTPUT:
-c     jet momenta:                             pjet(0:3,nexternal), E is 0th cmpnt
-c     the number of jets (with pt > SYCUT):    njet
-c     the jet for a given particle 'i':        jet(i),   note that this is
-c     the particle in pQCD, which doesn't necessarily correspond to the particle
-c     label in the process
+c     jet momenta:                           pjet(0:3,nexternal), E is 0th cmpnt
+c     the number of jets (with pt > SYCUT):  njet
+c     the jet for a given particle 'i':      jet(i),   note that this is the
+c                                            particle in pQCD, which doesn't
+c                                            necessarily correspond to the particle
+c                                            label in the process
 c
          call amcatnlo_fastjetppgenkt_timed(pQCD,nQCD,rfj,sycut,palg,
-     &pjet,njet,jet)
+     $        pjet,njet,jet)
 c
 c******************************************************************************
 
 c Apply the maximal pseudo-rapidity cuts on the jets:      
          if (etaj.gt.0d0) then 
-c Boost the jets to the lab frame for the pseudo-rapidity cut
-            chybst=cosh(ybst_til_tolab)
-            shybst=sinh(ybst_til_tolab)
-            chybstmo=chybst-1.d0
-            do i=1,njet
-               call boostwdir2(chybst,shybst,chybstmo,xd,
-     &              pjet(0,i),pjetlab(0,i))
-            enddo
 c Count the number of jets that pass the pseud-rapidity cut
             njet_eta=0
             do i=1,njet
-               if (abs(eta(pjetlab(0,i))).lt.ETAJ) then
+               if (abs(eta(pjet(0,i))).lt.ETAJ) then
                   njet_eta=njet_eta+1
                endif
             enddo
@@ -312,132 +216,223 @@
 
 c Apply the jet cuts
          if (njet .ne. nQCD .and. njet .ne. nQCD-1) then
-            passcuts=.false.
+            passcuts_user=.false.
             return
          endif
       endif
-
- 122  continue
-
-c Begin photon isolation
-c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE 
-c   Use is made of parton cm frame momenta. If this must be
-c   changed, pQCD used below must be redefined
-c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE 
+c
+c PHOTON (ISOLATION) CUTS
+c
+c find the photons
+      do i=1,nexternal
+         if (istatus(i).eq.1 .and. ipdg(i).eq.22) then
+            is_a_ph(i)=.true.
+         else
+            is_a_ph(i)=.false.
+         endif
+      enddo
       if (ptgmin.ne.0d0) then
-        nph=0
-        do j=nincoming+1,nexternal
-          if (is_a_ph(j)) then
-            nph=nph+1
-            do i=0,3
-              pgamma(i,nph)=p(i,j) ! Use C.o.M. frame momenta
-            enddo
-          endif
-        enddo
-        if(nph.eq.0)goto 444
-
-        if(isoEM)then
-          nem=nph
-          do k=1,nem
-            do i=0,3
-              pem(i,k)=pgamma(i,k)
-            enddo
-          enddo
-          do j=nincoming+1,nexternal
-            if (is_a_lp(j).or.is_a_lm(j)) then
-              nem=nem+1
-              do i=0,3
-                pem(i,nem)=p(i,j) ! Use C.o.M. frame momenta
-              enddo
+         nph=0
+         do j=nincoming+1,nexternal
+            if (is_a_ph(j)) then
+               nph=nph+1
+               do i=0,3
+                  pgamma(i,nph)=p(i,j)
+               enddo
             endif
-          enddo
-        endif
-
-        alliso=.true.
-
-        j=0
-        dowhile(j.lt.nph.and.alliso)
+         enddo
+         if(nph.eq.0)goto 444
+         
+         if(isoEM)then
+            nem=nph
+            do k=1,nem
+               do i=0,3
+                  pem(i,k)=pgamma(i,k)
+               enddo
+            enddo
+            do j=nincoming+1,nexternal
+               if (is_a_lp(j).or.is_a_lm(j)) then
+                  nem=nem+1
+                  do i=0,3
+                     pem(i,nem)=p(i,j)
+                  enddo
+               endif
+            enddo
+         endif
+         
+         alliso=.true.
+
+         j=0
+         do while(j.lt.nph.and.alliso)
 c Loop over all photons
-          j=j+1
-
-          ptg=pt(pgamma(0,j))
-          if(ptg.lt.ptgmin)then
-            passcuts=.false.
-            return
-          endif
-
-          if (etagamma.gt.0d0) then
-c     for rapidity cut, boost this one gamma to the lab frame
-             call boostwdir2(chybst,shybst,chybstmo,xd,
-     &            pgamma(0,j),pgammalab)
-             if (abs(eta(pgammalab)).gt.etagamma) then
-                passcuts=.false.
-                return
-             endif
-          endif
-
+            j=j+1
+            
+            ptg=pt(pgamma(0,j))
+            if(ptg.lt.ptgmin)then
+               passcuts_user=.false.
+               return
+            endif
+            if (etagamma.gt.0d0) then
+               if (abs(eta(pgamma(0,j))).gt.etagamma) then
+                  passcuts_user=.false.
+                  return
+               endif
+            endif
+         
 c Isolate from hadronic energy
-          do i=1,nQCD
-            drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i)))
-          enddo
-          call sortzv(drlist,isorted,nQCD,ismode,isway,izero)
-          Etsum(0)=0.d0
-          nin=0
-          do i=1,nQCD
-            if(dble(drlist(isorted(i))).le.R0gamma)then
-              nin=nin+1
-              Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i)))
-            endif
-          enddo
-          do i=1,nin
-            alliso=alliso .and.
-     #        Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
-     #                                  R0gamma,xn,epsgamma,ptg)
-          enddo
-
-c Isolate from EM energy
-          if(isoEM.and.nem.gt.1)then
-            do i=1,nem
-              drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
+            do i=1,nQCD
+               drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i)))
             enddo
-            call sortzv(drlist,isorted,nem,ismode,isway,izero)
-c First of list must be the photon: check this, and drop it
-            if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
-              write(*,*)'Error #1 in photon isolation'
-              write(*,*)j,isorted(1),drlist(isorted(1))
-              stop
-            endif
+            call sortzv(drlist,isorted,nQCD,ismode,isway,izero)
             Etsum(0)=0.d0
             nin=0
-            do i=2,nem
-              if(dble(drlist(isorted(i))).le.R0gamma)then
-                nin=nin+1
-                Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
-              endif
+            do i=1,nQCD
+               if(dble(drlist(isorted(i))).le.R0gamma)then
+                  nin=nin+1
+                  Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i)))
+               endif
             enddo
             do i=1,nin
-              alliso=alliso .and.
-     #          Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
-     #                                    R0gamma,xn,epsgamma,ptg)
+               alliso=alliso .and.
+     $              Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
+     $              R0gamma,xn,epsgamma,ptg)
             enddo
-
-          endif
-
+            
+c Isolate from EM energy
+            if(isoEM.and.nem.gt.1)then
+               do i=1,nem
+                  drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
+               enddo
+               call sortzv(drlist,isorted,nem,ismode,isway,izero)
+c First of list must be the photon: check this, and drop it
+               if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
+                  write(*,*)'Error #1 in photon isolation'
+                  write(*,*)j,isorted(1),drlist(isorted(1))
+                  stop
+               endif
+               Etsum(0)=0.d0
+               nin=0
+               do i=2,nem
+                  if(dble(drlist(isorted(i))).le.R0gamma)then
+                     nin=nin+1
+                     Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
+                  endif
+               enddo
+               do i=1,nin
+                  alliso=alliso .and.
+     $               Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
+     $               R0gamma,xn,epsgamma,ptg)
+               enddo
+            endif
 c End of loop over photons
-        enddo
-
-        if(.not.alliso)then
-          passcuts=.false.
-          return
-        endif
-
- 444    continue
+         enddo
+         if(.not.alliso)then
+            passcuts_user=.false.
+            return
+         endif
+ 444     continue
 c End photon isolation
       endif
 
-
- 123  continue
-
+C***************************************************************
+C***************************************************************
+C PUT HERE YOUR USER-DEFINED CUTS
+C***************************************************************
+C***************************************************************
+C
+c$$$C EXAMPLE: cut on top quark pT
+c$$$      do i=1,nexternal   ! loop over all external particles
+c$$$         if (istatus(i).eq.1    ! final state particle
+c$$$     &        .and. abs(ipdg(i)).eq.6) then    ! top quark
+c$$$C apply the pT cut (pT should be large than 200 GeV for the event to
+c$$$C pass cuts)
+c$$$            if ( p(1,i)**2+p(2,i)**2 .lt. 200d0**2 ) then
+c$$$C momenta do not pass cuts. Set passcuts_user to false and return
+c$$$               passcuts_user=.false.
+c$$$               return
+c$$$            endif
+c$$$         endif
+c$$$      enddo
+c
+      return
+      end
+
+
+
+
+
+
+
+C***************************************************************
+C***************************************************************
+C NO NEED TO CHANGE ANY OF THE FUNCTIONS BELOW
+C***************************************************************
+C***************************************************************
+      logical function passcuts(p,rwgt)
+      implicit none
+      include "nexternal.inc"
+      include 'run.inc'
+      include 'genps.inc'
+      REAL*8 P(0:3,nexternal),rwgt
+      integer i,j,istatus(nexternal),iPDG(nexternal)
+c For boosts
+      double precision ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,
+     #                        sqrtshat,shat
+      double precision chybst,shybst,chybstmo
+      double precision xd(1:3)
+      data (xd(i),i=1,3)/0,0,1/
+c Momenta of the particles
+      double precision plab(0:3, nexternal),pp(0:4, nexternal)
+c Masses of external particles
+      double precision pmass(nexternal)
+      common/to_mass/pmass
+c PDG codes of particles
+      integer maxflow
+      parameter (maxflow=999)
+      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
+     &     icolup(2,nexternal,maxflow)
+      common /c_leshouche_inc/idup,mothup,icolup
+      logical passcuts_user
+      external passcuts_user
+c Make sure have reasonable 4-momenta
+      if (p(0,1) .le. 0d0) then
+         passcuts=.false.
+         return
+      endif
+c Also make sure there's no INF or NAN
+      do i=1,nexternal
+         do j=0,3
+            if(p(j,i).gt.1d32.or.p(j,i).ne.p(j,i))then
+               passcuts=.false.
+               return
+            endif
+         enddo
+      enddo
+      rwgt=1d0
+c Boost the momenta p(0:3,nexternal) to the lab frame plab(0:3,nexternal)
+      chybst=cosh(ybst_til_tolab)
+      shybst=sinh(ybst_til_tolab)
+      chybstmo=chybst-1.d0
+      do i=1,nexternal
+         call boostwdir2(chybst,shybst,chybstmo,xd,
+     &        p(0,i),plab(0,i))
+      enddo
+c Fill the arrays (momenta, status and PDG):
+      do i=1,nexternal
+         if (i.le.nincoming) then
+            istatus(i)=-1
+         else
+            istatus(i)=1
+         endif
+         do j=0,3
+            pp(j,i)=plab(j,i)
+         enddo
+         pp(4,i)=pmass(i)
+         ipdg(i)=idup(i,1)
+      enddo
+c Call the actual cuts function  
+      passcuts = passcuts_user(pp,istatus,ipdg)
       RETURN
       END
 
@@ -733,3 +728,116 @@
       iso_getdelphi=tmp
       return
       end
+
+
+
+      DOUBLE PRECISION FUNCTION R2_04(P1,P2)
+c************************************************************************
+c     Distance in eta,phi between two particles.
+c************************************************************************
+      IMPLICIT NONE
+c
+c     Arguments
+c
+      double precision p1(0:4),p2(0:4),p1a(0:3),p2a(0:3)
+      integer i
+c
+c     External
+c
+      double precision eta,DELTA_PHI
+      external eta,delta_phi
+c-----
+c  Begin Code
+c-----
+      do i=0,3
+         p1a(i)=p1(i)
+         p2a(i)=p2(i)
+      enddo
+      R2_04 = (DELTA_PHI(P1a,P2a))**2+(eta(p1a)-eta(p2a))**2
+      RETURN
+      END
+
+      double precision function pt_04(p)
+c************************************************************************
+c     Returns transverse momentum of particle
+c************************************************************************
+      IMPLICIT NONE
+c
+c     Arguments
+c
+      double precision p(0:4)
+c-----
+c  Begin Code
+c-----
+
+      pt_04 = dsqrt(p(1)**2+p(2)**2)
+
+      return
+      end
+
+
+      double precision function eta_04(p)
+c************************************************************************
+c     Returns pseudo rapidity of particle
+c************************************************************************
+      IMPLICIT NONE
+c
+c     Arguments
+c
+      double precision p(0:4),pa(0:3)
+      integer i
+c
+c     external
+c
+      double precision theta,tp,pi
+      parameter (pi=3.14159265358979323846264338327950d0)
+      external theta
+c-----
+c  Begin Code
+c-----
+      do i=0,3
+         pa(i)=p(i)
+      enddo
+      tp=theta(pa)
+      if (abs(tp).lt.1d-5) then
+         eta_04=25d0
+      elseif (abs(tp-pi).lt.1d-5) then
+         eta_04=-25d0
+      else
+         eta_04=-dlog(dtan(theta(pa)/2d0))
+      endif
+
+      return
+      end
+
+
+
+      DOUBLE PRECISION FUNCTION invm2_04(P1,P2,dsign)
+c************************************************************************
+c     Invarient mass of 2 particles
+c************************************************************************
+      IMPLICIT NONE
+c
+c     Arguments
+c
+      double precision p1(0:4),p2(0:4),dsign
+c
+c     Local
+c      
+      integer i
+      double precision ptot(0:3)
+c
+c     External
+c
+      double precision dot
+      external dot
+c-----
+c  Begin Code
+c-----
+
+      do i=0,3
+         ptot(i)=p1(i)+dsign*p2(i)
+      enddo
+      invm2_04 = dot(ptot,ptot)
+      RETURN
+      END

=== modified file 'Template/NLO/SubProcesses/cuts.inc'
--- Template/NLO/SubProcesses/cuts.inc	2013-11-21 14:36:04 +0000
+++ Template/NLO/SubProcesses/cuts.inc	2014-03-30 11:08:12 +0000
@@ -1,5 +1,5 @@
 C
-C     NEW CUTS
+C     DEFAULT CUTS FROM RUN_CARD>DAt
 C
       REAL*8 JETALGO,JETRADIUS,PTJ,ETAJ  ! jet cuts
       REAL*8 PTL,ETAL,DRLL,MLL           ! lepton cuts
@@ -10,74 +10,10 @@
       LOGICAL ISOEM                           ! photon cuts
       COMMON /TO_ISOGAMMA_CUTS/PTGMIN,R0GAMMA,XN,EPSGAMMA,ETAGAMMA
       COMMON /TO_ISOGAMMA_EM/ISOEM
-
-C
-C     AUTOMATIC CUTS
-C
-c      REAL*8 PTJ,PTB,PTA,PTL,PTONIUM  ! ptmin
-      REAL*8 PTB,PTA,PTONIUM  ! ptmin
-      REAL*8 EJ,EB,EA,EL,ETAONIUM      ! min energy
-c      REAL*8 ETAJ,ETAB,ETAA,ETAL ! min rapidity
-      REAL*8 ETAB,ETAA ! min rapidity
-c      REAL*8 DRJJ,DRBB,DRLL,DRAA,DRBJ,DRAJ,DRJL,DRAB,DRBL,DRAL ! min distance
-      REAL*8 DRJJ,DRBB,DRAA,DRBJ,DRAJ,DRJL,DRAB,DRBL,DRAL ! min distance
-      REAL*8 MMJJ,MMLL,MMAA,MMBB !min inv mass
-      COMMON /TO_AUTO_CUTS/
-c     &     PTJ ,PTB ,PTA ,PTL ,PTONIUM, ! pt
-     &     PTB ,PTA ,PTONIUM, ! pt
-     &     EJ  ,EB  ,EA  ,EL  , ! energy
-c     &     ETAJ,ETAB,ETAA,ETAL,ETAONIUM, ! rapidity 
-     &     ETAB,ETAA,ETAONIUM, ! rapidity 
-     &     DRJJ,DRAJ,DRBJ,DRJL, ! delta R
-     &     DRBB,DRAB,DRBL,
-     &     DRAA,DRAL,
-c     &     DRLL,
-     &     MMJJ,MMLL,MMAA,MMBB ! invariant masses of opposite sign
-
-C     
-C     MIN/MAX cuts
-C
-      REAL*8 PTJmax ,PTBmax ,PTAmax ,PTLmax    ! max pt
-      REAL*8 EJmax  ,EBmax  ,EAmax  ,ELmax     ! max rap
-      REAL*8 ETAJmin,ETABmin,ETAAmin,ETALmin   ! min rap
-      REAL*8 DRJJmax,DRBBmax,DRLLmax,DRAAmax
-      REAL*8 DRBJmax,DRAJmax,DRJLmax,DRABmax,DRBLmax,DRALmax
-      REAL*8 MMJJmax,MMLLmax,MMAAmax,MMBBmax !max inv mass 
-      REAL*8 MMNL,MMNLMAX ! invariant mass of all leptons
-      REAL*8 cutuse
-      REAL*8 ptj1min,ptj2min,ptj3min,ptj4min
-      REAL*8 ptj1max,ptj2max,ptj3max,ptj4max
-      REAL*8 ht2min,ht3min,ht4min
-      REAL*8 ht2max,ht3max,ht4max
-      REAL*8 htjmin,htjmax
-      REAL*8 misset,missetmax,ptheavy
       integer maxjetflavor
-	
-      COMMON /to_min_max_cuts/
-     &     PTJmax,PTBmax,PTAmax,PTLmax,
-     &     EJmax,EBmax,EAmax,ELmax,
-     &     ETAJmin,ETABmin,ETAAmin,ETALmin,
-     &     DRJJmax,DRBBmax,DRLLmax,DRAAmax,
-     &     DRBJmax,DRAJmax,DRJLmax,DRABmax,DRBLmax,DRALmax,
-     &     MMJJmax,MMLLmax,MMAAmax,MMBBmax,
-     &     MMNL,MMNLmax,
-     &     ptj1min,ptj1max,
-     &     ptj2min,ptj2max,
-     &     ptj3min,ptj3max,
-     &     ptj4min,ptj4max,
-     &     cutuse,
-     &     ht2min,ht3min,ht4min,
-     &     ht2max,ht3max,ht4max,
-     &     htjmin,htjmax,
-     &     misset,missetmax,ptheavy,
-     &     maxjetflavor
-
-C
-C     SPECIAL CUTS
-C
-C-- this are somewhat old kind of cuts that where used for specific
-c   purposes
-c
-      real*8 xptj,xptb,xpta,xptl,xmtc
-      real*8 xetamin,xqcut,deltaeta
-      common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
+      COMMON /to_min_max_cuts/maxjetflavor
+C
+C     SPECIAL CUTS (FOR MERGING)
+C
+      real*8 xmtc,xqcut
+      common /to_specxpt/xmtc,xqcut

=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2013-12-22 08:19:15 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2014-03-30 11:08:12 +0000
@@ -10,6 +10,7 @@
       parameter       (ZERO = 0d0)
       include 'nexternal.inc'
       include 'genps.inc'
+      include 'reweight.inc'
       INTEGER    ITMAX,   NCALL
 
       common/citmax/itmax,ncall
@@ -203,12 +204,15 @@
       call addfil(dum)
       if (imode.eq.-1.or.imode.eq.0) then
          if(imode.eq.0)then
+c Don't safe the reweight information when just setting up the grids.
+            doreweight=.false.
             do j=0,nintervals
                do i=1,ndimmax
                   xgrid(j,i)=0.d0
                enddo
             enddo
          else
+            doreweight=do_rwgt_scale.or.do_rwgt_pdf
 c to restore grids:
             open (unit=12, file='mint_grids',status='old')
             do j=0,nintervals
@@ -427,9 +431,6 @@
          endif
       enddo
 
-      sigint=0d0
-
-c Find the nFKSprocess for which we compute the Born-like contributions
       if (firsttime) then
          firsttime=.false.
          foundB(1)=.false.
@@ -450,7 +451,8 @@
          write (*,*) 'For the Born we use nFKSprocesses  #',
      &        nFKSprocessBorn
       endif
-         
+
+      sigint=0d0
 c
 c Compute the Born-like contributions with nbody=.true.
 c THIS CAN BE OPTIMIZED

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2013-12-22 08:19:15 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2014-03-30 11:08:12 +0000
@@ -10,6 +10,7 @@
       parameter       (ZERO = 0d0)
       include 'nexternal.inc'
       include 'genps.inc'
+      include 'reweight.inc'
       INTEGER    ITMAX,   NCALL
 
       common/citmax/itmax,ncall
@@ -102,8 +103,9 @@
       integer ifold_energy,ifold_phi,ifold_yij
       common /cifoldnumbers/ifold_energy,ifold_phi,ifold_yij
       logical putonshell
+      logical only_virt
       integer imode
-      common /c_imode/imode
+      common /c_imode/imode,only_virt
       logical unwgt
       double precision evtsgn
       common /c_unwgt/evtsgn,unwgt
@@ -181,6 +183,18 @@
       write(*,*) "getting user params"
       call get_user_params(ncall,itmax,iconfig,imode,
      &     ixi_i,iphi_i,iy_ij,SHsep)
+c Only do the reweighting when actually generating the events
+      if (imode.eq.2) then
+         doreweight=do_rwgt_scale.or.do_rwgt_pdf
+      else
+         doreweight=.false.
+      endif
+      if (abrv(1:4).eq.'virt') then
+         only_virt=.true.
+      else
+         only_virt=.false.
+      endif
+
       if(imode.eq.0)then
         flat_grid=.true.
       else
@@ -398,7 +412,8 @@
 
          weight=(ans(1)+ans(5))/ncall
 
-         if (abrv(1:3).ne.'all' .and. abrv(1:4).ne.'born') then
+         if (abrv(1:3).ne.'all' .and. abrv(1:4).ne.'born' .and.
+     $        abrv(1:4).ne.'virt') then
             write (*,*) 'CANNOT GENERATE EVENTS FOR ABRV',abrv
             stop 1
          endif
@@ -411,7 +426,7 @@
                vn=3
                call gen(sigintF,ndim,xgrid,ymax,ymax_virt,1,x,vn)
             else
-               if (ran2().lt.ans(5)/(ans(1)+ans(5))) then
+               if (ran2().lt.ans(5)/(ans(1)+ans(5)) .or. only_virt) then
                   abrv='virt'
                   vn=1
                   call gen(sigintF,ndim,xgrid,ymax,ymax_virt,1,x,vn)
@@ -602,9 +617,6 @@
       write(*,*)'  Enter alpha>0 to set G_azi=0 (no azi corr)'
       read(*,*)alazi,beazi
       write (*,*) 'for G_azi: alpha=',alazi,', beta=',beazi
-
-c$$$      write (*,*) "H-events (0), or S-events (1)"
-c$$$      read (*,*) i
       i=2
       if (i.eq.0) then
          Hevents=.true.
@@ -774,10 +786,8 @@
       integer sum
       parameter (sum=3)
       data firsttime /.true./
-      logical foundB(2),j_fks_initial(fks_configs),found_ini1,found_ini2
+      logical j_fks_initial(fks_configs),found_ini1,found_ini2
      $     ,found_fnl,j_fks_initial_found,j_fks_final_found
-      integer nFKSprocessBorn(2)
-      save nFKSprocessBorn,foundB
       double precision vol1,sigintR,res,f_tot,rfract
       integer proc_map(0:fks_configs,0:fks_configs)
      $     ,i_fks_proc(fks_configs),j_fks_proc(fks_configs)
@@ -802,20 +812,9 @@
 c Find the nFKSprocess for which we compute the Born-like contributions
       if (firsttime) then
          firsttime=.false.
-         foundB(1)=.false.
-         foundB(2)=.false.
          maxproc_save=0
          do nFKSprocess=1,fks_configs
             call fks_inc_chooser()
-            if (PDG_type(i_fks).eq.21) then
-               if (j_fks.le.nincoming) then
-                  foundB(1)=.true.
-                  nFKSprocessBorn(1)=nFKSprocess
-               else
-                  foundB(2)=.true.
-                  nFKSprocessBorn(2)=nFKSprocess
-               endif
-            endif
 c Set Bjorken x's to some random value before calling the dlum() function
             xbk(1)=0.5d0
             xbk(2)=0.5d0
@@ -828,50 +827,8 @@
             endif
          enddo
          write (*,*) 'Total number of FKS directories is', fks_configs
-         write (*,*) 'For the Born we use nFKSprocesses  #',
-     &        nFKSprocessBorn
 c For sum over identical FKS pairs, need to find the identical structures
-         if (sum.eq.0) then
-c MC over FKS directories (1 FKS directory per nbody PS point)
-            proc_map(0,0)=fks_configs
-            do i=1,fks_configs
-               proc_map(i,0)=1
-               proc_map(i,1)=i
-            enddo
-         elseif (sum.eq.1) then
-c Sum over FKS directories (all FKS directories per nbody PS point)
-            proc_map(0,0)=1
-            proc_map(1,0)=fks_configs
-            do i=1,fks_configs
-               proc_map(1,i)=i
-            enddo
-         elseif (sum.eq.2) then
-c Sum over all FKS pairs that have the same i_fks and j_fks
-            proc_map(0,0)=0
-            do i=1,fks_configs
-               proc_map(i,0)=0
-               i_fks_proc(i)=0
-               j_fks_proc(i)=0
-            enddo
-            do nFKSprocess=1,fks_configs
-               call fks_inc_chooser()
-               i=1
-               do while ( i.le.proc_map(0,0) )
-                  if (i_fks.eq.i_fks_proc(i)
-     &              .and. j_fks.eq.j_fks_proc(i) ) then
-                     exit
-                  endif
-                  i=i+1
-               enddo
-               proc_map(i,0)=proc_map(i,0)+1
-               proc_map(i,proc_map(i,0))=nFKSprocess
-               if (i.gt.proc_map(0,0)) then
-                  proc_map(0,0)=proc_map(0,0)+1
-                  i_fks_proc(proc_map(0,0))=i_fks
-                  j_fks_proc(proc_map(0,0))=j_fks
-               endif
-            enddo
-         elseif (sum.eq.3) then
+         if (sum.eq.3) then
 c MC over FKS pairs that have soft singularity
             proc_map(0,0)=0
             do i=1,fks_configs
@@ -987,26 +944,6 @@
                   endif
                endif
             enddo
-         elseif(sum.eq.4) then
-c Sum over all j_fks initial (final) state
-            proc_map(0,0)=0
-            do i=1,2
-               proc_map(i,0)=0
-            enddo
-            j_fks_initial_found=.false.
-            j_fks_final_found=.false.
-            do nFKSprocess=1,fks_configs
-               call fks_inc_chooser()
-               if (j_fks.gt.nincoming) then
-                  proc_map(1,0)=proc_map(1,0)+1
-                  proc_map(1,proc_map(1,0))=nFKSprocess
-               else
-                  proc_map(2,0)=proc_map(2,0)+1
-                  proc_map(2,proc_map(2,0))=nFKSprocess
-               endif
-            enddo
-            if (proc_map(1,0).ne.0) proc_map(0,0)=proc_map(0,0)+1
-            if (proc_map(2,0).ne.0) proc_map(0,0)=proc_map(0,0)+1
          else
             write (*,*) 'sum not know in driver_mintMC.f',sum
          endif
@@ -1049,7 +986,7 @@
             double_check(i)=0d0
          enddo
          do i=1,99
-            if (abrv.eq.'grid'.or.abrv.eq.'born'.or.abrv(1:2).eq.'vi')
+            if (abrv.eq.'grid'.or.abrv.eq.'born')
      &           then
                if(i.le.ndim-3)then
                   x(i)=xx(i)
@@ -1070,29 +1007,10 @@
 c Compute the Born-like contributions with nbody=.true.
 c     
          call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
-         nFKSprocess=proc_map(proc_map(0,1),1) ! just pick the first
-                                               ! because it only matters
-                                               ! which parton is j_fks
+c Pick the first one because that's the one with the soft singularity.
+         nFKSprocess=proc_map(proc_map(0,1),1)
          nFKSprocess_all=nFKSprocess
          call fks_inc_chooser()
-         if (j_fks.le.nincoming) then
-            if (.not.foundB(1)) then
-               write(*,*) 'Trying to generate Born momenta with '/
-     &              /'initial state j_fks, but there is no '/
-     &              /'configuration with i_fks a gluon and j_fks '/
-     &              /'initial state'
-               stop
-            endif
-            nFKSprocess=nFKSprocessBorn(1)
-         else
-            if (.not.foundB(2)) then
-               write(*,*) 'Trying to generate Born momenta with '/
-     &              /'final state j_fks, but there is no configuration'/
-     &              /' with i_fks a gluon and j_fks final state'
-               stop
-            endif
-            nFKSprocess=nFKSprocessBorn(2)
-         endif
          nbody=.true.
          fillh=.false. ! this is set to true in BinothLHA if doing MC over helicities
          nFKSprocess_used=nFKSprocess
@@ -1168,7 +1086,7 @@
          f_check=f(2)
          if (f_check.ne.0d0.or.sigintF.ne.0d0) then
             if (abs(sigintF-f_check)/max(abs(f_check),abs(sigintF))
-     $           .gt.1d-1) then
+     $           .gt.1d0) then
                write (*,*) 'Error inaccuracy in unweight table 1'
      $              ,sigintF,f_check
                stop
@@ -1265,8 +1183,9 @@
       common/cproc_combination/iproc_save,eto,etoi,maxproc_found
       double precision virtual_over_born
       common/c_vob/virtual_over_born
+      logical only_virt
       integer imode
-      common /c_imode/imode
+      common /c_imode/imode,only_virt
 c Trivial check on the Born contribution
       do i=1,iproc_save(nFKSprocess_used_born)
          if (unwgt_table(0,2,i).ne.0d0) then
@@ -1372,8 +1291,9 @@
                      do j=1,iproc_save(nFKSprocess_used_born)
                         if (eto(j,nFKSprocess_used_born).eq.i) then
 c If when computing upper bounding envelope (imode.eq.1) do not include
-c the virtual corrections, because a separate bound is computed for them
-                           if (imode.eq.1) then
+c the virtual corrections, because a separate bound is computed for them.
+c Exception: when computing only the virtual, do include it here!
+                           if (imode.eq.1 .and. .not. only_virt) then
                               f_unwgt(nFKSprocess_soft,i) =
      $                             f_unwgt(nFKSprocess_soft,i) +
      $                             unwgt_table(0,1,i)+unwgt_table(0,2,i)
@@ -1593,7 +1513,8 @@
                if (eto(j,nFKSprocess_used_born).eq.i) then
 c If when computing upper bounding envelope (imode.eq.1) do not include
 c the virtual corrections, because a separate bound is computed for them
-                  if (imode.eq.1) then
+c Exception: when computing only the virtual, do include it here!
+                  if (imode.eq.1 .and. .not. only_virt) then
                      f_unwgt(nFKSprocess_used_born,i)=
      &                    f_unwgt(nFKSprocess_used_born,i)+
      &                    unwgt_table(0,1,i)

=== modified file 'Template/NLO/SubProcesses/fill_MC_mshell.f'
--- Template/NLO/SubProcesses/fill_MC_mshell.f	2012-09-03 13:51:10 +0000
+++ Template/NLO/SubProcesses/fill_MC_mshell.f	2014-03-30 11:08:12 +0000
@@ -3,13 +3,13 @@
       integer i
 c Monte Carlo masses: use PDG conventions.
 c May be given in input eventually
-      double precision mcmass(-5:21)
+      double precision mcmass(-16:21)
       common/cmcmass/mcmass
 c
       character*10 MonteCarlo
       common/cMonteCarloType/MonteCarlo
 c
-      do i=-5,21
+      do i=-16,21
         mcmass(i)=-1.d10
       enddo
       if (MonteCarlo.eq.'HERWIG6') then
@@ -29,5 +29,8 @@
       do i=-5,-1
          mcmass(i)=mcmass(-i)
       enddo
+      do i=-16,-11
+         mcmass(i)=mcmass(-i)
+      enddo
       return      
       end

=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2014-02-18 07:43:12 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2014-03-30 11:08:12 +0000
@@ -257,7 +257,7 @@
      #                 prefact_deg_sxi,prefact_deg_slxi,deg_wgt,deg_swgt,
      #                 deg_xi_c,deg_lxi_c,deg_xi_sc,deg_lxi_sc,
      #                 cnt_swgt,cnt_wgt,xlum_ev,xlum_c,xlum_s,xlum_sc,xsec,
-     #                 bpower
+     #                 bpower,cpower
       integer i,j,iplot
 
       integer izero,ione,itwo,mohdr,iplot_ev,iplot_cnt,iplot_born
@@ -406,6 +406,7 @@
       external setclscales,rewgt
 
       double precision pmass(nexternal)
+      double precision rwgt_muR_dep_fac
       include "pmass.inc"
 
       vegas_weight=vegaswgt
@@ -437,6 +438,23 @@
 c Put here call to compute bpower
          call compute_bpower(p_born,bpower)
          wgtbpower=bpower
+
+c Compute cpower done for bottom Yukawa, routine needs to be adopted
+c for other muR-dependendent factors
+         call compute_cpower(p_born,cpower)
+         if(dabs(cpower+1d0).lt.tiny) then
+            wgtcpower=0d0
+         else
+            wgtcpower=cpower
+         endif
+c Check that things are done consistently
+         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'
+           stop
+         endif
+
          firsttime=.false.
       endif
 
@@ -541,7 +559,8 @@
            call sreal(pp,xi_i_fks_ev,y_ij_fks_ev,fx_ev)
            xlum_ev = dlum()
            xsec = fx_ev*s_ev*ffact*wgt*prefact*rwgt
-           ev_wgt = xlum_ev*xsec
+           ev_wgt = xlum_ev*xsec * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
            if(doreweight)then
              wgtmuR2(1)=muR2_current/muR_over_ref**2
              wgtmuF12(1)=muF12_current/muF1_over_ref**2
@@ -609,12 +628,14 @@
             call sreal(p1_cnt(0,1,1),xi_i_fks_cnt(ione),one,fx_c)
             xlum_c = dlum()
             xsec = fx_c*s_c*jac_cnt(1)*(prefact_c+prefact_coll)*rwgt
-            cnt_wgt_c=cnt_wgt_c-xlum_c*xsec
+            cnt_wgt_c=cnt_wgt_c-xlum_c*xsec * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             call sreal_deg(p1_cnt(0,1,1),xi_i_fks_cnt(ione),one,
      #                     deg_xi_c,deg_lxi_c)
             deg_wgt=deg_wgt+( deg_xi_c+deg_lxi_c*log(xi_i_fks_cnt(ione)) )*
      #                      jac_cnt(1)*prefact_deg*rwgt/(shat/(32*pi**2))*
-     #                      xlum_c
+     #                      xlum_c * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             iplot=1
             if(doreweight)then
               call reweight_fillkin(pp,ithree)
@@ -653,7 +674,11 @@
               xsec=fx_s*s_s*jac_cnt(0)
               cnt_s=xlum_s*xsec
               cnt_wgt_s=cnt_wgt_s-cnt_s*prefact*rwgt
+     f             * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               cnt_swgt_s=cnt_swgt_s-cnt_s*prefact_cnt_ssc*rwgt
+     f             * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               if(doreweight)
      #          wgtwreal(2)=-xsec*(prefact+prefact_cnt_ssc)*rwgt/
      #                      g**(nint(2*wgtbpower+2.d0))
@@ -688,10 +713,15 @@
                 wgtwnsmur(2)=wgtwnstmpmur*xsec/g**(nint(2*wgtbpower
      &               +2.d0))
               endif
-              bsv_wgt=bsv_wgt*xnormsv
-              virt_wgt=virt_wgt*xnormsv
-              born_wgt=born_wgt*xnormsv
+              bsv_wgt=bsv_wgt*xnormsv * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
+              virt_wgt=virt_wgt*xnormsv * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
+              born_wgt=born_wgt*xnormsv * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               born_wgt_ao2pi=born_wgt_ao2pi*xnormsv
+     f             * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             endif
  548        continue
             iplot=0
@@ -719,17 +749,23 @@
             xsec=fx_sc*s_sc*jac_cnt(2)
             cnt_sc=xlum_sc*xsec
             cnt_wgt_sc=cnt_wgt_sc+cnt_sc*(prefact_c+prefact_coll)*rwgt
+     f           * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             cnt_swgt_sc=cnt_swgt_sc+
      &           cnt_sc*(prefact_cnt_ssc_c+prefact_coll_c)*rwgt
+     f           * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             call sreal_deg(p1_cnt(0,1,2),zero,one,
      #                     deg_xi_sc,deg_lxi_sc)
             deg_wgt=deg_wgt-( deg_xi_sc+deg_lxi_sc*log(xi_i_fks_cnt(ione)) )*
      #                     jac_cnt(2)*prefact_deg*rwgt/(shat/(32*pi**2))*
-     #                     xlum_sc
+     #                     xlum_sc * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             deg_swgt=deg_swgt-( deg_xi_sc*prefact_deg_sxi +
      #                     deg_lxi_sc*prefact_deg_slxi )*
      #                     jac_cnt(2)*rwgt/(shat/(32*pi**2))*
-     #                     xlum_sc
+     #                     xlum_sc * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             if(iplot.ne.0)iplot=2
             if(doreweight)then
               call reweight_fillkin(pp,ifour)
@@ -1009,7 +1045,7 @@
      &     ,gfactcl,xmcMC ,xmcME,SxmcMC,SxmcME,HxmcMC,HxmcME, xlum_c
      &     ,xlum_s,xlum_sc ,xlum_mc,xlum_mc_save, dummy,Sev_wgt,Hev_wgt
      &     ,fx_ev,probne ,sevmc ,xlum_ev,get_ptrel, xlum_mc_fact,xnormsv
-     &     ,xsec,bpower ,dsigS ,dsigH,totH_wgt,virt_wgt
+     &     ,xsec,bpower,cpower ,dsigS ,dsigH,totH_wgt,virt_wgt
       integer i,j
 
       integer izero,ione,itwo,mohdr,iplot_ev,iplot_cnt,iplot_born
@@ -1201,6 +1237,8 @@
       double precision rewgt
       external setclscales,rewgt
 
+      double precision rwgt_muR_dep_fac
+
       double precision pmass(nexternal)
       include "pmass.inc"
 
@@ -1305,6 +1343,22 @@
 c Put here call to compute bpower
          call compute_bpower(p_born,bpower)
          wgtbpower=bpower
+
+c Compute cpower done for bottom Yukawa, routine needs to be adopted
+c for other muR-dependendent factors
+         call compute_cpower(p_born,cpower)
+         if(dabs(cpower+1d0).lt.tiny) then
+            wgtcpower=0d0
+         else
+            wgtcpower=cpower
+         endif
+c Check that things are done consistently
+        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'
+           stop
+         endif
       endif
 c
       prefact=xinorm_ev/xi_i_fks_ev*
@@ -1454,10 +1508,12 @@
             do j=1,IPROC
                unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess,1
      &              ,j)+xmcxsec(i)*PD(j)*xlum_mc_fact*sevmc*wgt*prefact
-     &              *rwgt*CONV
+     &              *rwgt*CONV * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
                unwgt_table(nFKSprocess,2,j)=unwgt_table(nFKSprocess,2
      &              ,j)-xmcxsec(i)*PD(j)*xlum_mc_fact*sevmc*wgt*prefact
-     &              *rwgt*CONV
+     &              *rwgt*CONV * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             enddo
             if(doreweight)then
                wgtwmcxsec_all(i,nFKSprocess*2)=-xsec*xlum_mc_fact
@@ -1471,8 +1527,10 @@
             endif
           endif
         enddo
-        SxmcMC=xmcMC*sevmc*wgt*prefact*rwgt
-        HxmcMC=-xmcMC*sevmc*wgt*prefact*rwgt
+        SxmcMC=xmcMC*sevmc*wgt*prefact*rwgt * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
+        HxmcMC=-xmcMC*sevmc*wgt*prefact*rwgt * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
       endif
 c
       if( (.not.flagmc).and.gfactsf.eq.1.d0 .and.
@@ -1560,6 +1618,8 @@
             do j=1,IPROC
                unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess,1
      &              ,j)+xsec*PD(j)*(1-gfactsf)*probne*CONV
+     f              * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here          
             enddo
             if ((gfactsf.lt.1.d0.and.gfactcl.lt.1.d0 .and.
      &           probne.gt.0.d0) .and. pmass(j_fks).eq.0.d0) then
@@ -1567,6 +1627,8 @@
                do j=1,IPROC
                   unwgt_table(nFKSprocess,2,j)=unwgt_table(nFKSprocess
      &                 ,2,j)-xsec*PD(j)*(1-gfactsf)*probne*CONV
+     f                 * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here             
                enddo
             endif
             if(doreweight) then
@@ -1578,21 +1640,26 @@
             endif
             if( y_ij_fks_ev.gt.1d0-deltaS )then
                xsec = fx_c*s_c*jac_cnt(1)*(prefact_c+prefact_coll)*rwgt
-               cnt_wgt_c=cnt_wgt_c-xlum_c*xsec
+               cnt_wgt_c=cnt_wgt_c-xlum_c*xsec * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
                do j=1,IPROC
                   unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess
-     &                 ,1,j)-xsec*PD(j)*CONV
+     &                 ,1,j)-xsec*PD(j)*CONV * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
                enddo
                call sreal_deg(p1_cnt(0,1,1),xi_i_fks_cnt(ione),one,
      #                        deg_xi_c,deg_lxi_c)
                deg_wgt=deg_wgt+( deg_xi_c+deg_lxi_c*log(xi_i_fks_cnt(ione)) )*
      #                         jac_cnt(1)*prefact_deg*rwgt/(shat/(32*pi**2))*
-     #                         xlum_c
+     #                         xlum_c * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
                do j=1,IPROC
                   unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess
      &                 ,1,j)+PD(j)*( deg_xi_c+deg_lxi_c
      &                 *log(xi_i_fks_cnt(ione)) )* jac_cnt(1)
      &                 *prefact_deg*rwgt/(shat/(32*pi**2))*CONV
+     f                 * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here             
                enddo
                if(doreweight)then
                   wgtwreal_all(3,nFKSprocess*2-1)=wgtwreal_all(3
@@ -1655,6 +1722,8 @@
             do j=1,IPROC
                unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess,1
      &              ,j)+xsec*PD(j)*(1-gfactsf)*probne*CONV
+     f              * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here          
             enddo
             if(doreweight)then
                wgtwreal_all(2,nFKSprocess*2-1)=xsec*(1-gfactsf)*probne/
@@ -1665,6 +1734,8 @@
                do j=1,IPROC
                   unwgt_table(nFKSprocess,2,j)=unwgt_table(nFKSprocess
      &                 ,2,j)-xsec*PD(j)*(1-gfactsf)*probne*CONV
+     f                 * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
                enddo
                if(doreweight)wgtwreal_all(2,nFKSprocess*2)=
      &              xsec/g**(nint(2*wgtbpower+2.d0))
@@ -1673,11 +1744,16 @@
               xsec=fx_s*s_s*jac_cnt(0)
               cnt_s=xlum_s*xsec
               cnt_wgt_s=cnt_wgt_s-cnt_s*prefact*rwgt
+     f             * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               cnt_swgt_s=cnt_swgt_s-cnt_s*prefact_cnt_ssc*rwgt
+     f             * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               do j=1,IPROC
                  unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess,1
      $                ,j)-PD(j)*xsec*(prefact+prefact_cnt_ssc)*rwgt
-     $                *CONV
+     $                *CONV * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               enddo
               if(doreweight)wgtwreal_all(2,nFKSprocess*2-1)
      &             =wgtwreal_all(2,nFKSprocess*2-1)-xsec*(prefact
@@ -1714,15 +1790,22 @@
               endif
               do j=1,IPROC
                  unwgt_table(0,1,j)=unwgt_table(0,1,j)+PD(j)*bsv_wgt
-     &                *xsec*CONV
+     &                *xsec*CONV * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
                  unwgt_table(0,3,j)=unwgt_table(0,3,j)+PD(j)*virt_wgt
-     &                *xsec*CONV
+     &                *xsec*CONV * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
                  unwgt_table(1,3,j)=unwgt_table(1,3,j)+PD(j)*born_wgt
      &                *xsec*CONV*g**2/(8d0*PI**2)
+     f                * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               enddo
-              bsv_wgt=bsv_wgt*xnormsv
-              virt_wgt=virt_wgt*xnormsv
-              born_wgt=born_wgt*xnormsv
+              bsv_wgt=bsv_wgt*xnormsv * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
+              virt_wgt=virt_wgt*xnormsv * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
+              born_wgt=born_wgt*xnormsv * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             endif
  548        continue
          endif
@@ -1762,6 +1845,8 @@
             do j=1,IPROC
                unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess,1
      &              ,j)-PD(j)*xsec*(1-gfactsf)*probne*CONV
+     f              * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             enddo
             if(doreweight)then
                wgtwreal_all(4,nFKSprocess*2-1)=-xsec*(1-gfactsf)*probne/
@@ -1773,6 +1858,8 @@
                do j=1,IPROC
                   unwgt_table(nFKSprocess,2,j)=unwgt_table(nFKSprocess
      &                 ,2,j)+PD(j)*xsec*(1-gfactsf)*probne*CONV
+     f                 * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here             
                enddo
                if(doreweight)wgtwreal_all(4,nFKSprocess*2)=
      &              -xsec/g**(nint(2*wgtbpower+2.d0))
@@ -1783,18 +1870,24 @@
               xsec=fx_sc*s_sc*jac_cnt(2)
               cnt_sc=xlum_sc*xsec
               cnt_wgt_sc=cnt_wgt_sc+cnt_sc*(prefact_c+prefact_coll)*rwgt
+     f             * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               cnt_swgt_sc=cnt_swgt_sc+
      &             cnt_sc*(prefact_cnt_ssc_c+prefact_coll_c)*rwgt
+     f             * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               call sreal_deg(p1_cnt(0,1,2),zero,one,
      #                       deg_xi_sc,deg_lxi_sc)
               deg_wgt=deg_wgt-
      #                    ( deg_xi_sc+deg_lxi_sc*log(xi_i_fks_cnt(ione)) )*
      #                    jac_cnt(2)*prefact_deg*rwgt/(shat/(32*pi**2))*
-     #                    xlum_sc
+     #                    xlum_sc * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               deg_swgt=deg_swgt-( deg_xi_sc*prefact_deg_sxi +
      #                       deg_lxi_sc*prefact_deg_slxi )*
      #                       jac_cnt(2)*rwgt/(shat/(32*pi**2))*
-     #                       xlum_sc
+     #                       xlum_sc * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               do j=1,IPROC
                  unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess,1
      &                ,j)+PD(j)*(xsec*(prefact_c+prefact_coll
@@ -1803,7 +1896,8 @@
      &                jac_cnt(2)*prefact_deg*rwgt/(shat/(32*pi**2))
      &                -(deg_xi_sc*prefact_deg_sxi + deg_lxi_sc
      &                *prefact_deg_slxi )* jac_cnt(2)*rwgt/(shat/(32*pi
-     &                **2)))*CONV
+     &                **2)))*CONV * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               enddo
               if(doreweight)then
                  wgtwreal_all(4,nFKSprocess*2-1)=wgtwreal_all(4
@@ -1823,8 +1917,10 @@
            endif
         endif
       endif
-      SxmcME=SxmcME*(1-gfactsf)*probne
-      HxmcME=-HxmcME*(1-gfactsf)*probne
+      SxmcME=SxmcME*(1-gfactsf)*probne * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
+      HxmcME=-HxmcME*(1-gfactsf)*probne * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
       if(doreweight)then
          xsec=(1-gfactsf)*probne
          wgtwreal_all(2,nFKSprocess*2)=
@@ -1879,13 +1975,17 @@
           call sreal(pp,xi_i_fks_ev,y_ij_fks_ev,fx_ev)
           xlum_ev = dlum()
           xsec = fx_ev*s_ev*ffact*wgt*prefact*rwgt
-          Sev_wgt = xlum_ev*xsec*(1-probne)
-          Hev_wgt = xlum_ev*xsec*probne
+          Sev_wgt = xlum_ev*xsec*(1-probne) * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
+          Hev_wgt = xlum_ev*xsec*probne * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
           do j=1,IPROC
              unwgt_table(nFKSprocess,1,j)=unwgt_table(nFKSprocess,1,j)
-     &            +PD(j)*xsec*(1-probne)*CONV
+     &            +PD(j)*xsec*(1-probne)*CONV * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
              unwgt_table(nFKSprocess,2,j)=unwgt_table(nFKSprocess,2,j)
-     &            +PD(j)*xsec*probne*CONV
+     &            +PD(j)*xsec*probne*CONV * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
           enddo
           if(doreweight)then
              if(ifill1H.eq.0)then
@@ -4242,8 +4342,8 @@
       parameter (ComputePoles=.false.)
       parameter (fksprefact=.true.)
 
-      double precision beta0
-      common/cbeta0/beta0
+      double precision beta0,ren_group_coeff
+      common/cbeta0/beta0,ren_group_coeff
 
       logical calculatedBorn
       common/ccalculatedBorn/calculatedBorn
@@ -4470,20 +4570,20 @@
 
 c eq.(MadFKS.C.13)
          if(abrv.eq.'viSA'.or.abrv.eq.'viSB')then
-           bsv_wgt=bsv_wgt + 2*pi*beta0*wgtbpower*log(shat/QES2)*
-     #                       ao2pi*dble(wgt1(1))
+           bsv_wgt=bsv_wgt + 2*pi*(beta0*wgtbpower
+     #      +ren_group_coeff*wgtcpower)*log(shat/QES2)*ao2pi*dble(wgt1(1))
          elseif(abrv.eq.'novA'.or.abrv.eq.'novB')then
-           bsv_wgt=bsv_wgt + 2*pi*beta0*wgtbpower*log(q2fact(1)/shat)*
-     #                       ao2pi*dble(wgt1(1))
+           bsv_wgt=bsv_wgt + 2*pi*(beta0*wgtbpower
+     #      +ren_group_coeff*wgtcpower)*log(q2fact(1)/shat)*ao2pi*dble(wgt1(1))
          elseif(abrv.ne.'virt' .and. abrv.ne.'viSC' .and.
      #          abrv.ne.'viLC')then
-           bsv_wgt=bsv_wgt + 2*pi*beta0*wgtbpower*log(q2fact(1)/QES2)*
-     #                       ao2pi*dble(wgt1(1))
+           bsv_wgt=bsv_wgt + 2*pi*(beta0*wgtbpower
+     #      +ren_group_coeff*wgtcpower)*log(q2fact(1)/QES2)*ao2pi*dble(wgt1(1))
          endif
 c eq.(MadFKS.C.14)
          if(abrv(1:2).ne.'vi')then
-           bsv_wgt=bsv_wgt - 2*pi*beta0*wgtbpower*
-     #                       log(q2fact(1)/scale**2)*ao2pi*dble(wgt1(1))
+           bsv_wgt=bsv_wgt - 2*pi*(beta0*wgtbpower
+     #      +ren_group_coeff*wgtcpower)*log(q2fact(1)/scale**2)*ao2pi*dble(wgt1(1))
          endif
 
 
@@ -4507,7 +4607,8 @@
                  endif
                enddo
                wgtwnstmpmuf=ao2pi*wgtwnstmpmuf*dble(wgt1(1))
-               wgtwnstmpmur=2*pi*beta0*wgtbpower*ao2pi*dble(wgt1(1))
+               wgtwnstmpmur=2*pi*(beta0*wgtbpower
+     #      +ren_group_coeff*wgtcpower)*ao2pi*dble(wgt1(1))
              endif
 c bsv_wgt here always contains the Born; must subtract it, since 
 c we need the pure NLO terms only
@@ -4652,6 +4753,102 @@
       return
       end
 
+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 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)
+      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
+
+      logical calculatedBorn
+      common/ccalculatedBorn/calculatedBorn
+
+      double precision tiny
+      parameter (tiny=1d-6)
+c comment these lines to calculate cpower
+      cpower = -1d0
+      return
+c comment these lines to calculate cpower
+
+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:
+c$$$      born_wgt=dble(wgt1(1))
+c$$$      
+c$$$c Multiply the Yukawa by 10 (If you use this,
+c$$$c double check that GC_33 is the yukawa! (also below))
+c$$$      if (GC_33.ne.0d0) then
+c$$$         GC_33 = GC_33 * 10d0
+c$$$      else
+c$$$         write(*,*)'Warning In Bornsoftvirtual'
+c$$$         Write(*,*)'Yukawa Is Zero - Cpower Set To Zero'
+c$$$         Cpower = 0d0
+c$$$         Return
+c$$$      Endif
+c$$$
+c$$$c recompute the Born with the new Yukawa
+c$$$      calculatedBorn=.false.
+c$$$      call sborn(p_born,wgt1)
+c$$$
+c$$$c Compute cpower
+c$$$      cpower=Log10(dble(wgt1(1))/born_wgt)
+c$$$      if(abs(cpower-dble(nint(cpower))) .gt. tiny) then
+c$$$         write(*,*)'Error in computation of cpower:'
+c$$$         write(*,*)' not an integer',cpower
+c$$$         stop
+c$$$      elseif (cpower.lt.-tiny) then
+c$$$         write(*,*)'Error in computation of cpower:'
+c$$$         write(*,*)' negative value',cpower
+c$$$         stop
+c$$$      else
+c$$$c set it to the integer exactly
+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$$$     &        abs(cpower-wgtcpower).gt.tiny )then
+c$$$c$$$            write(*,*)'Error in compute_cpower'
+c$$$c$$$            write(*,*)'cpower(s) are:',cpower,wgtcpower
+c$$$c$$$            stop
+c$$$c$$$         endif
+c$$$      endif
+c$$$
+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$$$
+c$$$      return
+      end
+
 
       subroutine eikonal_Ireg(p,m,n,xicut_used,eikIreg)
       implicit none
@@ -5130,6 +5327,11 @@
 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
 
@@ -5144,8 +5346,8 @@
       double precision c(0:1),gamma(0:1),gammap(0:1)
       common/fks_colors/c,gamma,gammap
 
-      double precision beta0
-      common/cbeta0/beta0
+      double precision beta0,ren_group_coeff
+      common/cbeta0/beta0,ren_group_coeff
 
       logical softtest,colltest
       common/sctests/softtest,colltest
@@ -5253,6 +5455,8 @@
             
 c Beta_0 defined according to (MadFKS.C.5)
       beta0=gamma(0)/(2*pi)
+c ren_group_coeff defined accordingly
+      ren_group_coeff=ren_group_coeff_in/(2*pi)
 
       if (firsttime_nFKSprocess(nFKSprocess)) then
          firsttime_nFKSprocess(nFKSprocess)=.false.

=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
--- Template/NLO/SubProcesses/genps_fks.f	2013-12-22 08:19:15 +0000
+++ Template/NLO/SubProcesses/genps_fks.f	2014-03-30 11:08:12 +0000
@@ -1689,7 +1689,7 @@
         if(rat.gt.-tiny)then
           tmp=0.d0
         else
-          write(6,*)'Error #2 in function Lambda:',s,ma2,mb2
+          write(6,*)'Error #2 in function Lambda:',s,ma2,mb2,rat
         endif
       endif
       LAMBDA=tmp
@@ -2061,7 +2061,11 @@
      &              (qwidth(i)+cBW_width(i,-1)) ! b(-1) is negative here
                b(-1)=qmass(i)+b(-1)*qwidth(i)
                b(-1)=b(-1)**2
-               if (x(-i).lt.0.5d0) then
+
+               if (b(-1).gt.smax) then
+                   s(i)=(smax-smin)*x(-i)+smin
+                   xjac0=xjac0*(smax-smin)
+               elseif (x(-i).lt.0.5d0) then
                   x0=2d0*x(-i)
                   s(i)=(b(-1)-smin)*x0+smin
                   xjac0=2d0*xjac0*(b(-1)-smin)

=== modified file 'Template/NLO/SubProcesses/handling_lhe_events.f'
--- Template/NLO/SubProcesses/handling_lhe_events.f	2013-12-12 16:05:18 +0000
+++ Template/NLO/SubProcesses/handling_lhe_events.f	2014-03-30 11:08:12 +0000
@@ -45,7 +45,7 @@
       subroutine write_lhef_header_banner(ifile,nevents,MonteCarlo,path)
       implicit none 
       integer ifile,nevents,iseed,i,pdf_set_min,pdf_set_max,idwgt
-      double precision mcmass(-5:21),rw_Rscale_down,rw_Rscale_up
+      double precision mcmass(-16:21),rw_Rscale_down,rw_Rscale_up
      $     ,rw_Fscale_down,rw_Fscale_up
       character*10 MonteCarlo
       character*100 path
@@ -159,6 +159,9 @@
       do i=1,5
          write (ifile,'(2x,i6,3x,e12.6)')i,mcmass(i)
       enddo
+      write (ifile,'(2x,i6,3x,e12.6)')11,mcmass(11)
+      write (ifile,'(2x,i6,3x,e12.6)')13,mcmass(13)
+      write (ifile,'(2x,i6,3x,e12.6)')15,mcmass(15)
       write (ifile,'(2x,i6,3x,e12.6)')21,mcmass(21)
       write(ifile,'(a)') '  </MonteCarloMasses>'
 c Write here the reweight information if need be
@@ -280,7 +283,7 @@
       character*80 string,string0
       character*3 event_norm
       common/cevtnorm/event_norm
-      double precision temp,remcmass(-5:21)
+      double precision temp,remcmass(-16:21)
       common/cremcmass/remcmass
       ipart=-1000000
       nevents = -1
@@ -309,7 +312,7 @@
           dowhile( index(string,'</montecarlomasses>').eq.0 .and.
      #             index(string,'</MonteCarloMasses>').eq.0 )
             read(string,*)ipart,temp
-            if(ipart.lt.-5.or.ipart.gt.21)then
+            if(ipart.lt.-16.or.ipart.gt.21)then
               write(*,*)'Error in read_lhef_header:'
               write(*,*)' incomprehensible list of parton masses',ipart
               stop
@@ -495,14 +498,15 @@
             write(ifile,442)wgtwmcxsecE(i),
      #                      wgtmcxbjE(1,i),wgtmcxbjE(2,i)
           enddo
-          if(jwgtinfo.eq.4) write(ifile,'(1x,e14.8,1x,i4,1x,i4)')
-     &         wgtbpower,nFKSprocess_used,nFKSprocess_used_born
+          if(jwgtinfo.eq.4) write(ifile,
+     f         '(1x,e14.8,1x,e14.8,1x,i4,1x,i4)')
+     &       wgtbpower,wgtcpower,nFKSprocess_used,nFKSprocess_used_born
           write(ifile,'(a)') '  </rwgt>'
          elseif(jwgtinfo.eq.5) then
            write(ifile,'(a)')'  <rwgt>'
            if (iSorH_lhe.eq.1) then ! S-event
-              write(ifile,'(1x,e14.8,i4,i4)') wgtbpower,nScontributions
-     $             ,i_process
+              write(ifile,'(1x,e14.8,1x,e14.8,i4,i4)') 
+     f             wgtbpower,wgtcpower,nScontributions,i_process
               write(ifile,'(1x,i4,1x,e14.8)') nFKSprocess_used_born
      &             ,wgtref_nbody_all(i_process)
               do i=1,mexternal
@@ -545,7 +549,8 @@
      $                ,iFKS),wgtmuF22_all(2,iFKS)
               enddo
            elseif (iSorH_lhe.eq.2) then ! H-event
-              write(ifile,'(1x,e14.8,i4)') wgtbpower,i_process
+              write(ifile,'(1x,e14.8,1x,e14.8,i4)')
+     f             wgtbpower,wgtcpower,i_process
               iFKS=nFKSprocess_used*2
               write(ifile,'(1x,i4)') nFKSprocess_used
               write(ifile,'(1x,e14.8,1x,i4)') wgtref_all(iFKS,i_process)
@@ -695,13 +700,15 @@
             read(ifile,442)wgtwmcxsecE(i),
      #                     wgtmcxbjE(1,i),wgtmcxbjE(2,i)
           enddo
-          if(jwgtinfo.eq.4) read(ifile,'(1x,e14.8,1x,i4,1x,i4)')
-     &         wgtbpower,nFKSprocess_used,nFKSprocess_used_born
+          if(jwgtinfo.eq.4) read(ifile,
+     f     '(1x,e14.8,1x,e14.8,1x,i4,1x,i4)')
+     &     wgtbpower,wgtcpower,nFKSprocess_used,nFKSprocess_used_born
           read(ifile,'(a)')string
         elseif(jwgtinfo.eq.5) then
            read(ifile,'(a)')string
            if (iSorH_lhe.eq.1) then ! S-event
-              read(ifile,'(1x,e14.8,i4,i4)') wgtbpower,nScontributions
+              read(ifile,'(1x,e14.8,1x,e14.8,i4,i4)')
+     f             wgtbpower,wgtcpower,nScontributions
      $             ,i_process
               read(ifile,'(1x,i4,1x,e14.8)') nFKSprocess_used_born
      &             ,wgtref_nbody_all(i_process)
@@ -747,7 +754,8 @@
      $                ,iFKS),wgtmuF22_all(2,iFKS)
               enddo
            elseif (iSorH_lhe.eq.2) then ! H-event
-              read(ifile,'(1x,e14.8,i4)') wgtbpower,i_process
+              read(ifile,'(1x,e14.8,1x,e14.8,i4)')
+     f             wgtbpower,wgtcpower,i_process
               read(ifile,'(1x,i4)') nFKSprocess_used
               iFKS=nFKSprocess_used*2
               read(ifile,'(1x,e14.8,1x,i4)') wgtref_all(iFKS,i_process)
@@ -908,13 +916,15 @@
             read(ifile,442)wgtwmcxsecE(i),
      #                     wgtmcxbjE(1,i),wgtmcxbjE(2,i)
           enddo
-          if(jwgtinfo.eq.4) read(ifile,'(1x,e14.8,1x,i4,1x,i4)')
-     &         wgtbpower,nFKSprocess_used,nFKSprocess_used_born
+          if(jwgtinfo.eq.4) read(ifile,
+     f      '(1x,e14.8,1x,e14.8,1x,i4,1x,i4)')
+     &      wgtbpower,wgtcpower,nFKSprocess_used,nFKSprocess_used_born
           read(ifile,'(a)')string
         elseif(jwgtinfo.eq.5) then
            read(ifile,'(a)')string
            if (iSorH_lhe.eq.1) then ! S-event
-              read(ifile,'(1x,e14.8,i4,i4)') wgtbpower,nScontributions
+              read(ifile,'(1x,e14.8,1x,e14.8,i4,i4)') 
+     f             wgtbpower,wgtcpower,nScontributions
      $             ,i_process
               read(ifile,'(1x,i4,1x,e14.8)') nFKSprocess_used_born
      &             ,wgtref_nbody_all(i_process)
@@ -960,7 +970,8 @@
      $                ,iFKS),wgtmuF22_all(2,iFKS)
               enddo
            elseif (iSorH_lhe.eq.2) then ! H-event
-              read(ifile,'(1x,e14.8,i4)') wgtbpower,i_process
+              read(ifile,'(1x,e14.8,1x,e14.8,i4)') 
+     f             wgtbpower,wgtcpower,i_process
               read(ifile,'(1x,i4)') nFKSprocess_used
               iFKS=nFKSprocess_used*2
               read(ifile,'(1x,e14.8,1x,i4)') wgtref_all(iFKS,i_process)
@@ -1063,13 +1074,13 @@
 
 
       subroutine fill_MC_mshell_wrap(MC,masses)
-      double precision mcmass(-5:21),masses(-5:21)
+      double precision mcmass(-16:21),masses(-16:21)
       common/cmcmass/mcmass
       character*10 MonteCarlo,MC
       common/cMonteCarloType/MonteCarlo
       MonteCarlo=MC
       call fill_MC_mshell()
-      do i=-5,21
+      do i=-16,21
          masses(i)=mcmass(i)
       enddo
       return

=== modified file 'Template/NLO/SubProcesses/madfks_plot.f'
--- Template/NLO/SubProcesses/madfks_plot.f	2013-11-28 14:01:26 +0000
+++ Template/NLO/SubProcesses/madfks_plot.f	2014-03-30 11:08:12 +0000
@@ -6,7 +6,10 @@
       integer nwgt,max_weight
       parameter (max_weight=maxscales*maxscales+maxpdfs+1)
       character*15 weights_info(max_weight)
-      integer i,npdfs
+      integer i,npdfs,ii,jj,n
+      double precision xsecScale_acc(maxscales,maxscales)
+     $     ,xsecPDFr_acc(0:maxPDFs)
+      common /scale_pdf_print/xsecScale_acc,xsecPDFr_acc
       nwgt=1
       weights_info(nwgt)="central value  "
       if (do_rwgt_scale) then
@@ -48,12 +51,24 @@
          nwgt=nwgt+npdfs
       endif
       call analysis_begin(nwgt,weights_info)
+c To keep track of the accumulated results:
+      do ii=1,numscales
+         do jj=1,numscales
+            xsecScale_acc(jj,ii)=0d0
+         enddo
+      enddo
+      do n=0,npdfs
+         xsecPDFr_acc(n)=0d0
+      enddo
       return
       end
 
 
       subroutine topout
       implicit none
+      include 'reweight0.inc'
+      include 'reweightNLO.inc'
+      integer ii,jj,n
       logical usexinteg,mint
       common/cusexinteg/usexinteg,mint
       integer itmax,ncall
@@ -61,6 +76,9 @@
       logical useitmax
       common/cuseitmax/useitmax
       real*8 xnorm
+      double precision xsecScale_acc(maxscales,maxscales)
+     $     ,xsecPDFr_acc(0:maxPDFs)
+      common /scale_pdf_print/xsecScale_acc,xsecPDFr_acc
 c
       if(usexinteg.and..not.mint) then
          xnorm=1.d0/float(itmax)
@@ -71,6 +89,24 @@
       endif
       if(useitmax)xnorm=xnorm/float(itmax)
       call analysis_end(xnorm)
+c Write the accumulated results to a file
+      open (unit=34,file='scale_pdf_dependence.dat',status='unknown')
+      if (.not.useitmax) xnorm=xnorm/float(itmax)
+      write (34,*) numscales**2
+      if (numscales.gt.0) then
+         write (34,*) ((xsecScale_acc(ii,jj)*xnorm,ii=1
+     $        ,numscales),jj=1,numscales)
+      else
+         write (34,*) ''
+      endif
+      if (numPDFs.gt.0) then
+         write (34,*) numPDFs
+         write (34,*) (xsecPDFr_acc(n)*xnorm,n=0,numPDFs-1)
+      else
+         write(34,*) numPDFs
+         write (34,*) ''
+      endif
+      close(34)
       return                
       end
 
@@ -117,7 +153,10 @@
       common /c_leshouche_inc/idup,mothup,icolup
       integer nwgt,max_weight
       parameter (max_weight=maxscales*maxscales+maxpdfs+1)
-      double precision wgts(max_weight),wgtden
+      double precision wgts(max_weight),wgtden,ratio
+      double precision xsecScale_acc(maxscales,maxscales)
+     $     ,xsecPDFr_acc(0:maxPDFs)
+      common /scale_pdf_print/xsecScale_acc,xsecPDFr_acc
 c Born, n-body or (n+1)-body contribution:
       if(itype.eq.11) then
          ibody=1 ! (n+1)-body
@@ -156,20 +195,49 @@
 c The weights comming from reweighting:
       nwgt=1
       wgts(1)=www
+      if (wgtden.eq.0d0) then
+         if (www.eq.0d0) then
+            ratio=0d0
+         else
+            if (doreweight) then
+               write (*,*) 'ERROR in madfks_plot.f', wgtden,www
+            else
+               ratio=1d0
+            endif
+         endif
+      else
+c this ratio should essentially be the weight from vegas
+         ratio=www/wgtden
+      endif
       if (do_rwgt_scale) then
          do i=1,numscales
             do j=1,numscales
                nwgt=nwgt+1
-               wgts(nwgt)=wgtNLOxsecmu(ibody,i,j)*www/wgtden
+               wgts(nwgt)=wgtNLOxsecmu(ibody,i,j)*ratio
             enddo
          enddo
       endif
       if (do_rwgt_pdf) then
-         do i=1,numPDFs
+         do i=1,numPDFs-1 ! exclude the central set, so 'numPDFs-1'
             nwgt=nwgt+1
-            wgts(nwgt)=wgtNLOxsecPDF(ibody,i)*www/wgtden
+            wgts(nwgt)=wgtNLOxsecPDF(ibody,i)*ratio
          enddo
       endif
       call analysis_fill(p,istatus,ipdg,wgts,ibody)
+c Fill the accumulated results
+      if (do_rwgt_scale) then
+         do i=1,numscales
+            do j=1,numscales
+               xsecScale_acc(i,j)=xsecScale_acc(i,j)+
+     &              wgtNLOxsecmu(ibody,i,j)*ratio
+            enddo
+         enddo
+      endif
+      if (do_rwgt_pdf) then
+         xsecPDFr_acc(0)=xsecPDFr_acc(0)+www
+         do i=1,numPDFs-1
+            xsecPDFr_acc(i)=xsecPDFr_acc(i)+wgtNLOxsecPDF(ibody,i)*ratio
+         enddo
+      endif
  999  return      
       end

=== modified file 'Template/NLO/SubProcesses/reweight.f'
--- Template/NLO/SubProcesses/reweight.f	2013-12-11 10:33:34 +0000
+++ Template/NLO/SubProcesses/reweight.f	2014-03-30 11:08:12 +0000
@@ -267,6 +267,9 @@
       integer ipdg(n_max_cl),ipart(2,n_max_cl)
       logical isjet
       external isjet
+      integer iddgluon, iddother, idgluon, idother
+      logical isqcd
+      external isqcd
 
       idmo=ipdg(imo)
       idda1=ipdg(ida1)
@@ -349,19 +352,43 @@
             ipart(1,imo)=ipart(1,ida2)
             ipart(2,imo)=ipart(1,ida1)
          endif
-      else if(idmo.eq.21)then
-c     gluon -> higgs gluon
-         if ( idda1.eq.21 ) then
-            ipart(1,imo)=ipart(1,ida1)
-            ipart(2,imo)=ipart(2,ida1)
-         elseif ( idda2.eq.21 ) then
-            ipart(1,imo)=ipart(1,ida2)
-            ipart(2,imo)=ipart(2,ida2)
-         else
-            write (*,*) "ERROR in ipartupdate: unknown coupling",idmo
-     $           ,idda1,idda2
-            stop
-         endif
+      else if(idmo.eq.21.and.(idda1.eq.21.or.idda2.eq.21))then
+         if(idda1.eq.21) then
+            iddgluon = idda1
+            idgluon = ida1
+            iddother = idda2
+            idother = ida2
+         else
+            iddgluon = idda2
+            iddother = idda1
+            idgluon = ida2
+            idother = ida1
+         endif
+         if (isqcd(idother))then
+c        gluon -> gluon + scalar octet Choose hardest one
+            if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
+     $         p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
+               ipart(1,imo)=ipart(1,ida1)
+               ipart(2,imo)=ipart(2,ida1)
+            else
+               ipart(1,imo)=ipart(1,ida2)
+               ipart(2,imo)=ipart(2,ida2)
+            endif
+         else
+c        gluon -> gluon + Higgs use the gluon one
+               ipart(1,imo)=ipart(1,idgluon)
+               ipart(2,imo)=ipart(2,idgluon)
+         endif
+      else if(idmo.eq.21) then
+c     gluon > octet octet Choose hardest one
+            if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
+     $         p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
+               ipart(1,imo)=ipart(1,ida1)
+               ipart(2,imo)=ipart(2,ida1)
+            else
+               ipart(1,imo)=ipart(1,ida2)
+               ipart(2,imo)=ipart(2,ida2)
+            endif
       else if(idmo.eq.idda1.or.idmo.eq.idda1+sign(1,idda2))then
 c     quark -> quark-gluon or quark-Z or quark-h or quark-W
          ipart(1,imo)=ipart(1,ida1)
@@ -422,8 +449,8 @@
         return
       endif        
 c     Final State clustering
-      if(isjet(idda1).and.(isjet(idmo).or.idmo.eq.idda2).or.
-     $   isjet(idda2).and.(isjet(idmo).or.idmo.eq.idda1)) then
+      if((isjet(idda1).and.(isjet(idmo).or.idmo.eq.idda2)).or.(
+     $   isjet(idda2).and.(isjet(idmo).or.idmo.eq.idda1))) then
          isjetvx=.true.
       else
          isjetvx=.false.
@@ -814,7 +841,7 @@
 c     Flag mother as good jet if PDG is jet and both daughters are jets
             goodjet(imocl(n))=
      $           (isjet(ipdgcl(imocl(n),igraphs(1),nFKSprocess)).and.
-     $           goodjet(idacl(n,1)).and.goodjet(idacl(n,1)))
+     $           goodjet(idacl(n,1)).and.goodjet(idacl(n,2)))
          endif
       enddo
 

=== modified file 'Template/NLO/SubProcesses/reweight.inc'
--- Template/NLO/SubProcesses/reweight.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/SubProcesses/reweight.inc	2014-03-30 11:08:12 +0000
@@ -4,11 +4,7 @@
       include 'reweight1.inc'
 
       logical doreweight
-      parameter (doreweight=.true.)
-
-c Save reweight information for NLO computations
-      logical doNLOreweight
-      parameter (doNLOreweight=.true.)
+      common/c_doreweight/doreweight
 
 c Check that save weights are consistent with the reference weight wgtref
       logical check_reweight

=== modified file 'Template/NLO/SubProcesses/reweight0.inc'
--- Template/NLO/SubProcesses/reweight0.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/SubProcesses/reweight0.inc	2014-03-30 11:08:12 +0000
@@ -94,3 +94,33 @@
 c nFKSprocess used for the nbody and the .not.nbody
       integer nFKSprocess_used,nFKSprocess_used_Born
       common /c_nFKSproc_file/ nFKSprocess_used,nFKSprocess_used_Born
+
+c cpower: power of muR-dependent factor in the born
+      double precision wgtcpower
+      common/cwgtcpower/wgtcpower
+
+c input of cpower (checked against calculated value)
+      double precision cpowerinput
+      parameter (cpowerinput=0.d0)
+
+c switch for running muR-dependent factor runfac=1(running)/0(fixed)
+      integer runfac
+      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
+      double precision ren_group_coeff_in
+      parameter (ren_group_coeff_in=0d0)
+
+

=== modified file 'Template/NLO/SubProcesses/reweight_xsec.f'
--- Template/NLO/SubProcesses/reweight_xsec.f	2013-11-28 14:01:26 +0000
+++ Template/NLO/SubProcesses/reweight_xsec.f	2014-03-30 11:08:12 +0000
@@ -623,7 +623,6 @@
       include 'run.inc'
       include "reweight.inc"
 
-      logical passcuts
       double precision compute_rwgt_wgt_NLO
       double precision xmuR_over_ref,xmuF1_over_ref,
      #                 xmuF2_over_ref,xQES_over_ref
@@ -649,6 +648,7 @@
       logical setclscales
       double precision rewgt
       external setclscales,rewgt
+      double precision rwgt_muR_dep_fac
 c
       save_murrat=muR_over_ref
       save_muf1rat=muF1_over_ref
@@ -667,10 +667,13 @@
       xsec11=0.d0
       xsec12=0.d0
       xsec20=0.d0
+
+      if (wgtwreal(1).eq.0d0) goto 541
+
       call set_cms_stuff(mohdr)
       if( (kwgtinfo.eq.1.and.wgtmuR2(1).ne.0.d0) .or.
-     #    ((kwgtinfo.ge.3.or.kwgtinfo.le.5).and.
-     #     passcuts(wgtkin(0,1,1),rwgt)) )then
+     $     ((kwgtinfo.ge.3.or.kwgtinfo.le.5).and.wgtkin(0,1,1).gt.0.d0)
+     $     )then
         if(kwgtinfo.eq.1)then
           scale=muR_over_ref*sqrt(wgtmuR2(1))
           g=sqrt(4d0*pi*alphas(scale))
@@ -707,13 +710,25 @@
           nFKSprocess=nFKSprocess_used
           xlum = dlum()
           xsec11=xsec11+xlum*wgtwreal(1)*g**(2*wgtbpower+2.d0)
+     f         * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
         endif
       endif
 c
+ 541  continue
+      if ( wgtwreal(2).eq.0d0 .and.
+     $     wgtwreal(3).eq.0d0 .and. wgtwreal(4).eq.0d0 .and.
+     $     wgtwdeg(2).eq.0d0 .and.
+     $     wgtwdeg(3).eq.0d0 .and. wgtwdeg(4).eq.0d0 .and.
+     $     wgtwdegmuf(2).eq.0d0 .and.
+     $     wgtwdegmuf(3).eq.0d0 .and. wgtwdegmuf(4).eq.0d0 .and.
+     $     wgtwborn(2).eq.0d0 .and. wgtwns(2).eq.0d0 .and.
+     $     wgtwnsmuf(2).eq.0d0 .and. wgtwnsmur(2).eq.0d0) goto 542
+
       call set_cms_stuff(izero)
       if( (kwgtinfo.eq.1.and.wgtmuR2(2).ne.0.d0) .or.
-     #    ((kwgtinfo.ge.3.or.kwgtinfo.le.5).and.
-     #     passcuts(wgtkin(0,1,2),rwgt)) )then
+     $     ((kwgtinfo.ge.3.or.kwgtinfo.le.5).and.wgtkin(0,1,2).gt.0.d0)
+     $     )then
         if(kwgtinfo.eq.1)then
           scale=muR_over_ref*sqrt(wgtmuR2(2))
           g=sqrt(4d0*pi*alphas(scale))
@@ -737,12 +752,6 @@
           stop
         endif
         QES2_local=wgtqes2(2)
-        if(abs(QES2/QES2_local-1.d0).gt.tiny.and.
-     &       (kwgtinfo.ge.3.or.kwgtinfo.le.5))then
-          write(*,*)'Error in compute_rwgt_wgt_NLO'
-          write(*,*)' Mismatch in ES scale',QES2,QES2_local
-          stop
-        endif
         if(QES2_local.eq.0.d0)then
           if(wgtwdegmuf(3).ne.0.d0.or.
      #       wgtwdegmuf(4).ne.0.d0.or.
@@ -757,6 +766,12 @@
           xlgmuf=0.d0
           xlgmur=0.d0
         else
+           if(abs(QES2/QES2_local-1.d0).gt.tiny.and.
+     &          (kwgtinfo.ge.3.or.kwgtinfo.le.5))then
+              write(*,*)'Error in compute_rwgt_wgt_NLO'
+              write(*,*)' Mismatch in ES scale',QES2,QES2_local
+              stop
+           endif
           xlgmuf=log(q2fact(1)/QES2_local)
           xlgmur=log(scale**2/QES2_local)
         endif
@@ -784,24 +799,32 @@
             xlum = dlum()
             xsec12=xsec12+xlum*( wgtwreal(k)+wgtwdeg(k)+
      #                           wgtwdegmuf(k)*xlgmuf )*
-     #                g**(2*wgtbpower+2.d0)
+     #                g**(2*wgtbpower+2.d0) * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             if(k.eq.2)then
                nFKSprocess=nFKSprocess_used_born
                xlum = dlum()
               if(wgtbpower.gt.0)then
                 xsec20=xsec20+xlum*wgtwborn(k)*g**(2*wgtbpower)
+     f                * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               else
-                xsec20=xsec20+xlum*wgtwborn(k)
+                xsec20=xsec20+xlum*wgtwborn(k) * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
               endif
               xsec12=xsec12+xsec20
               xsec12=xsec12+xlum*( wgtwns(k)+
      #                         wgtwnsmuf(k)*xlgmuf+
      #                         wgtwnsmur(k)*xlgmur )*
-     #                  g**(2*wgtbpower+2.d0)
+     #                  g**(2*wgtbpower+2.d0) * rwgt_muR_dep_fac(scale)
+com-- muR-dependent fac is reweighted here
             endif
           endif
         enddo
       endif
+
+ 542  continue
+
 c
       muR_over_ref=save_murrat
       muF1_over_ref=save_muf1rat
@@ -833,12 +856,11 @@
       include "reweight.inc"
       include 'nFKSconfigs.inc'
 
-      logical passcuts
       double precision compute_rwgt_wgt_Hev
       double precision xmuR_over_ref,xmuF1_over_ref,
      #                 xmuF2_over_ref,xQES_over_ref