Merge lp:~maddevelopers/mg5amcnlo/2.3.3_PY8_install_cmd into lp:~maddevelopers/mg5amcnlo/2.3.4

[Valentin Hirschi]

This branch features two things:

1) A more advanced implementation of the 'install' command for the more complicated common dependencies of MG5_aMC: lhapdf5, lhapdf6, Pythia8, hepmc as well as the 'mg5amc_pythia8_interface'

2) A streamlined interface at LO to the showering of the produced .lhe events with Pythia8.

I shortly describe each aspect:

1)
This advanced installation is done via a separate set of installer scripts, hosted on (please quickly review these as well):

lp:~maddevelopers/mg5amcnlo/HEPToolsInstallers

These installers are automatically updated on the MG servers, and when the user wants to install a dependency, he firsts downloads these installers which will do the job for him.
The advantage is that if we want to modify the installation routines (to make it work on some specific setup we had overlook or simply to update the version used for some dependency), we can directly update HEPToolsInstallers without having to release a new version of MG5.
The script I wrote in HEPToolsInstallers is flexible so that it will be easy in the future to automatically install more complicated dependencies. In essence, it is a mini-package manager.
As a by-product of this implementation, it now also possible to install from madgraph the following dependencies (with the same 'install' command):

boost, zlib and hepmc

Special care was given to nicely handle the mutual dependency of Pythia8 and MG5 to LHAPDF. In other terms, when installing pythia8 automatically, MG5 makes several check and efforts in order to insure that Pythia8 is linked to the same LHAPDF version as the one specified in MG5 configuration.

2) When launching a process, one no longer turn 'pythia=ON/OFF' but instead modify the parameter 'SHOWER=' which can take the value 'PYTHIA6' or 'PYTHIA8' depending on what is installed (thanks Olivier for that). Same for detector simulation.

Both MLM and CKKW-L merging are supported as well as the automatic variation of the merging scale to evaluate systematics. This used to be done by SysCalc in Pythia8 and is now entirely performed by the Pythia8 shower driver. However, the choice of the variational merging scales to consider still happens via the run_card parameter 'sys_matchscale'.

Speaking of which, this driver is hosted on the launchpad branch:
   lp:~maddevelopers/mg5amcnlo/MG5aMC_PY8_interface
for which synchronized copy is kept on MadGraph servers and which is used for the function
'install mg5amc_pythia8_interface'
Int his way, both us and the pythia8 developers can update it without having to rerelease either Pythia8 or MadGraph5. Version compatibility by a dedicated check upon both the installation and usage of the this interface.

This implementation should also serve as a basis for future improvements for handling Pythia8.
In particular, the implementation of the pythia8 card in the class PYCard of 'madgraph/various/banner.py' should allow to easily tweak PY8 parameters (and their visibility to the user) in future developments. For example, the 'subruns' structure of the pythia8 card is fully supported (i.e. listing several .lhe to be processed within one single card, and each with specific parameters)

In order to test this branch you can start by simply doing:

MG5_aMC>install hepmc
MG5_aMC>install pythia8
MG5_aMC>install mg5amc_py8_interface

(normally the installer would automatically install hepmc and pythia8 for you when installing 'mg5amc_py8_interface' but it is better to force the installation by hand here otherwise my installer will find your system-wide installation of these tools and use it, even though they are not the development versions compatible with this LO interface)

Then you can do

MG5_aMC>generate p p > z > mu+ mu-
MG5_aMC>add process p p > z > mu+ mu- j
MG5_aMC>add process p p > z > mu+ mu- j j
MG5_aMC>output MyMLMTest
MG5_aMC>launch
>shower=PYTHIA8 (or roll the options [assuming you have installed pythia6 too] by pressing '1' several times)
And possibly change some options (commands below are not mandatory)
>set xqcut 12
>set sys_matchscale 20 30 50 60
>set JetMatching:nJetMax 2
>set JetMatching:qCut 40
(you can try to tab-complete 'set pythia8_card ' to see all PY8 parameter visible to the user and you can edit the card directly if you want to add one that MG5aMC doesn't know about)

Ok now you can launch, and if everything went well you should be able to follow the shower run live via its log file and eventually get the merging plots, directly available from within the HTML page.

[lp:~maddevelopers/mg5amcnlo/2.3.3_PY8_install_cmd]

355. By Valentin Hirschi on 2016-02-20

1. Temporarily commented out qCutMin since the PY8 driver is not yet ready for it.

356. By Valentin Hirschi on 2016-02-20

1. Added an IOTest for the histograms.py functionality of reading DJR plots and writing a gnuplot output
   for them.

[Paolo Torrielli (paolo-torrielli)]

Hi Valentin,
thanks for this interesting branch.

I tried to install lhapdf5 on my laptop (mac os 10.9), and it doesn't with the following error

---------------------------------------------------------------
   You are installing 'lhapdf5', please cite ref(s).
         [arXiv:0605240]
   when using results produced with this tool.
---------------------------------------------------------------
Downloading the HEPToolInstaller at:
   http://madgraph.hep.uiuc.edu/Downloads/HEPToolsInstaller/HEPToolsInstaller_V26.tar.gz
Now installing lhapdf5. Be patient...
Installer HEPToolInstaller.py is now processing the following command:
   /Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/HEPToolsInstallers/HEPToolInstaller.py lhapdf5 --prefix=/Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools --mg5_path=/Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd --mg5amc_py8_interface_tarball=http://madgraph.phys.ucl.ac.be/Downloads/MG5aMC_PY8_interface/MG5aMC_PY8_interface_V1.0.tar.gz
Downloading 'lhapdf5' sources...
Fetching data with command:
  curl -OL http://www.hepforge.org/archive/lhapdf/lhapdf-5.9.0.tar.gz
  % Total % Received % Xferd Average Speed Time Time Time Current
                                 Dload Upload Total Spent Left Speed
100 349 100 349 0 0 1354 0 --:--:-- --:--:-- --:--:-- 2624
100 1207k 100 1207k 0 0 537k 0 0:00:02 0:00:02 --:--:-- 746k
Installing tool 'lhapdf5'...
(you can follow the installation progress by running the command below in a separate terminal)
  tail -f /Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/lhapdf5/lhapdf5_install.log
Traceback (most recent call last):
  File "/Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/HEPToolsInstallers/HEPToolInstaller.py", line 715, in <module>
    install_with_dependencies(target_tool,is_main_target=True)
  File "/Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/HEPToolsInstallers/HEPToolInstaller.py", line 693, in install_with_dependencies
    exec('install_%s(tmp_path)'%target)
  File "<string>", line 1, in <module>
  File "/Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/HEPToolsInstallers/HEPToolInstaller.py", line 364, in install_lhapdf5
    stderr=lhapdf5_log)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/subprocess.py", line 524, in call
    return Popen(*popenargs, **kwargs).wait()
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/subprocess.py", line 711, in __init__
    errread, errwrite)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/subprocess.py", line 1308, in _execute_child
    raise child_exception
OSError: [Errno 2] No such file or directory
Command "install lhapdf5" interrupted with error:
InvalidCmd : Installation of lhapdf5 failed.
Please report this bug on https://bugs.launchpad.net/madgraph5
More information is found in 'MG5_debug'.
Please attach this file to your report.

Installing pythia8 does not work with the following error

MG5...

review: Needs Fixing

[Valentin Hirschi (valentin-hirschi)]

> Hi Valentin,
> thanks for this interesting branch.

Hi Paolo,

Thanks for having quickly started to look at this.

>
> I tried to install lhapdf5 on my laptop (mac os 10.9), and it doesn't with the
> following error
>
> ---------------------------------------------------------------
> You are installing 'lhapdf5', please cite ref(s).
> [arXiv:0605240]
> when using results produced with this tool.
> ---------------------------------------------------------------
> Downloading the HEPToolInstaller at:
> http://madgraph.hep.uiuc.edu/Downloads/HEPToolsInstaller/HEPToolsInstaller_
> V26.tar.gz
> Now installing lhapdf5. Be patient...
> Installer HEPToolInstaller.py is now processing the following command:
> /Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/
> HEPToolsInstallers/HEPToolInstaller.py lhapdf5 --prefix=/Users/paolotorrielli/
> WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools
> --mg5_path=/Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd -
> -mg5amc_py8_interface_tarball=http://madgraph.phys.ucl.ac.be/Downloads/MG5aMC_
> PY8_interface/MG5aMC_PY8_interface_V1.0.tar.gz
> Downloading 'lhapdf5' sources...
> Fetching data with command:
> curl -OL http://www.hepforge.org/archive/lhapdf/lhapdf-5.9.0.tar.gz
> % Total % Received % Xferd Average Speed Time Time Time
> Current
> Dload Upload Total Spent Left Speed
> 100 349 100 349 0 0 1354 0 --:--:-- --:--:-- --:--:-- 2624
> 100 1207k 100 1207k 0 0 537k 0 0:00:02 0:00:02 --:--:-- 746k
> Installing tool 'lhapdf5'...
> (you can follow the installation progress by running the command below in a
> separate terminal)
> tail -f /Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HE
> PTools/lhapdf5/lhapdf5_install.log
> Traceback (most recent call last):
> File "/Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPT
> ools/HEPToolsInstallers/HEPToolInstaller.py", line 715, in <module>
> install_with_dependencies(target_tool,is_main_target=True)
> File "/Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPT
> ools/HEPToolsInstallers/HEPToolInstaller.py", line 693, in
> install_with_dependencies
> exec('install_%s(tmp_path)'%target)
> File "<string>", line 1, in <module>
> File "/Users/paolotorrielli/WORK/Research/MG5_aMC/2.3.3_PY8_install_cmd/HEPT
> ools/HEPToolsInstallers/HEPToolInstaller.py", line 364, in install_lhapdf5
> stderr=lhapdf5_log)
> File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7
> /subprocess.py", line 524, in call
> return Popen(*popenargs, **kwargs).wait()
> File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7
> /subprocess.py", line 711, in __init__
> errread, errwrite)
> File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7
> /subprocess.py", line 1308, in _execute_child
> raise child_exception
> OSError: [Errno 2] No such file or directory
> Command "install lhapdf5" interrupted with error:
> InvalidCmd : Installation of lhapdf5 failed.
> Please report this bug ...

[Valentin Hirschi (valentin-hirschi)]

Dear reviewers,

We looked at the issue with Stefan and the error you were getting when installing Pythia8 was because of a part of both the main89-madgraph.cc example code and the MG5_aMC_PY8 interface drivers which was wrong (and not very well tested because it only gets compiled when using the non-hacked version of HepMC for writing named weights).

This means that it would have been fine if you had installed HEPMC with MG5aMC as well prior to your Pythia8 installation.
Anyway, now this issue has been fixed (both in the pythia8 tarball [which is lighter now that it has been stripped from the .lhe examples]) and in the online interface driver. So it should work with any version of HepMC.

Please give it another try and let me know.
Thanks again for your testing,

Cheers

[Rikkert Frederix (frederix)]

Hi Valentin,

This didn't fix the problem on my LINUX desktop machine. It still doesn't compile (problem with rsync); see e-mail I sent a couple of days ago.

On my mac laptop it now compiles correctly---apart from the fact that it cannot find any dependencies automatically. If I run the installer by hand, it does seem to find the dependencies:

[ Rikkert@MacBook-Air ~/Documents/Physics/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/HEPToolsInstallers ]$ ./HEPToolInstaller.py pythia8
Installer HEPToolInstaller.py is now processing the following command:
   ./HEPToolInstaller.py pythia8
'pythia8' dependency 'hepmc' found at:
  /usr/local
'pythia8' dependency 'zlib' found at:
  /usr
'pythia8' dependency 'lhapdf6' found at:
  /usr/local
Downloading 'pythia8' sources...

Cheers,
Rik

review: Needs Fixing

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

Very nice branch. It's really cool all the interesting way the code is improving those days.
I particularly like the change in the card edition part (maybe I did those but I do not remember)

Here is my review:
First some functionalities which does not work on this branch

1) The output pythia8 does not work in this branch. This need to be fixed for sure.
   - You change the comment for the PY8 path in the configuration file, is that related?

2) The standard plot routines (MA4) does not work anymore:
   - I do not have any plot either at parton level
   - I do not have any plot for after the shower (even if the lhe file in order to do such plot is created)

3) if I type the function shower, I have as question:
Which programs do you want to run?
    0 / auto : running existing cards
    1 / pythia : Pythia
    2 / pgs : Pythia + PGS
    3 / delphes : Pythia + Delphes.
 [0, 1, 2, 3, auto, pythia, pgs, delphes][60s to answer]
This need to be updated to the new format.
At the end of the run, the code ask me a second question:
Which programs do you want to run?
    0 / auto : running existing cards
    1 / pythia8 : Pythia8
    3 / delphes : Pythia8 + Delphes.
 [0, 1, 3, auto, pythia8, delphes][60s to answer]
 (This probably happens only if you have both PY6 and PY8 installed)

One general question (having those answer will allow me to review that part of the modification)

1) I have a bunch of question related to the CKKW-L support that you have added in MG:
   a) Where can I find some documentation for how to run CKKW-L?
      - In particular what is the meaning of is_pdg_for_merging_cut
      - Can I do p p > w+ [0-3]j, w+ > j j
        Is the jet comming from the W are correctly not handle by the Durham kt/ PT lund cut?
        (If so we have to update the cut_decays comment)
   b) Is the modification of durham kt cut retro compatible?
      - if ATLAS use a old run_card, will they have the same results?
   c) What is the official support for CKKWL?
      - in term of user visibility
      - in term of user support
      - in term of code validation
      - in term of long term plan

Following some more technical/code related remark

1) In the class common_run_interface you define the function:
def help_pythia8(self):
def help_shower(self):
I would naively expect to have those in the madevent interface (like check_pythia8)

2) I would suggest to refactorise the code between the following comment:
+ ########################################################################
+ ### Special setup of the Hidden parameters in the card for each MG type
+ # of run
in the do_pythia8 command
and to put those command in the PY8 card directly. This would make more sense and allow to
have a cleaner do_pythia8

3) Why do you recommend to install:
install mg5amc_py8_interface
if someone install Delphes2/pythia-pgs
should not it be "install pythia8"?
What is the logic behind this?

4) One problem with the automatic download of PY8/... couple to the auto-update of MG5
is that we can have recent version of MG5 with old version of PY8/...
I would propose to add a way to track the latest version online at the same l...

review: Needs Fixing

[Valentin Hirschi (valentin-hirschi)]

Hi Rik,

> Hi Valentin,
>
> This didn't fix the problem on my LINUX desktop machine. It still doesn't
> compile (problem with rsync); see e-mail I sent a couple of days ago.

This does not seem related to my branch, although you claim that you could install Pythia8.2 on this machine yourself, so it somehow must be related to what I do. Your problem seems to be:

rsync -a bin/* /tuph/t30/space/frederix/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/pythia8//bin --exclude .svn
rsync: chgrp "/tuph/t30/space/frederix/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/pythia8/bin/.pythia8-config.QThxt3" failed: Operation not permitted (1)
rsync error: some files/attrs were not transferred (see previous errors) (code 23) at main.c(1183) [sender=3.1.0]
make: *** [install] Error 23

So could I ask you to download the pythia8 tarball

http://slac.stanford.edu/~prestel/pythia82151.tar.gz

untar it somewhere and then install it by runnning:

./configure --prefix=/tuph/t30/space/frederix/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/pythia8/ --with-hepmc2=/space/rikkert/rivet/local --with-hepmc2-include=/space/rikkert/rivet/local/include --with-gzip=/tuph/t30/space/frederix/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/zlib --with-lhapdf6=/home/t30/ben/ga63zot/Documents/Physics/LHAPDF-6.1.6/build --with-lhapdf6-plugin=LHAPDF6.h --with-boost=default
make
make install

If it goes through this means that the problem comes from the fact that I do this installation from an automatically generated directory typically placed in /tmp/ or /var/ (it is done with 'tempfile.mkdtemp') and that you have only limited rights from there.
So I guess that if this is the problem I should simply use some directory directly in HEPTools then.

> On my mac laptop it now compiles correctly---apart from the fact that it
> cannot find any dependencies automatically. If I run the installer by hand, it
> does seem to find the dependencies:
>
> [ Rikkert@MacBook-Air
> ~/Documents/Physics/MG5_aMC/2.3.3_PY8_install_cmd/HEPTools/HEPToolsInstallers
> ]$ ./HEPToolInstaller.py pythia8
> Installer HEPToolInstaller.py is now processing the following command:
> ./HEPToolInstaller.py pythia8
> 'pythia8' dependency 'hepmc' found at:
> /usr/local
> 'pythia8' dependency 'zlib' found at:
> /usr
> 'pythia8' dependency 'lhapdf6' found at:
> /usr/local
> Downloading 'pythia8' sources...

Ok I added /usr/local/lib and /usr/lib as default lookup path, this should resolve the problem.
I don't want to update the servers before the issue above regarding rsync is solved, so don't test the above yet but tell me if you managed the installation by hand of the development version of Pythia 8 I specified.

>
> Cheers,
> Rik

[Valentin Hirschi (valentin-hirschi)]

Dea rreviewers,

We have still small fixes pending for MLM and for runs without SysCalc, so please as a first testing, I suggest you run the following commands on top of the installation of the mg5amc_pythia8_interface:

MG5_aMC>install SysCalc

Please report any failure of the automatic installation of SysCalc, because this is something that must work in this branch now (notice that this one does not proceed through the HEPToolsInstaller manager).

Once both SysCalc and the mg5amc_pythia8_interface are successfully installed, I suggest you try the following simulation (hence avoiding the pending problems in MLM and non-SysCalc runs).
Having your feedback on this type of runs would already help a lot towards making progress in this review:

MG5_aMC>generate p p > z > mu+ mu-
MG5_aMC>add process p p > z > mu+ mu- j
MG5_aMC>add process p p > z > mu+ mu- j j
MG5_aMC>output MyCKKWLTest
MG5_aMC>launch

Now turn on Pythia8 shower with the command below, or by rolling the options by pressing '1' several times:
(Several options are available if you have installed pythia6 too).

>shower=PYTHIA8
>[enter] or '0'

Now setup lhapdf (necessary for SysCalc)

>set pdlabel lhapdf
>set lhaid 11000
>set sys_pdf CT10nlo.LHgrid

And change the following run_card parameters so as to turn on CKKWL merging

>set xqcut 0
>set ickkw 0
>set ktdurham 20
>set use_syst True
>set sys_matchscale 20 30 50 60

And then the following Pythia8 shower card options related to the details of the CKKWL merging

>set Merging:Process 'pp>LEPTONS,NEUTRINOS'
>set Merging:TMS 40

You can now start the run as usual by typing enter or '0'.

Most of the time should then be spent in the shower which you should be able to follow from the log file, as suggested. Beware that it generates the heavy hepmc file as it goes.

Finally, you can assess the quality of the CKKW-L merging with the plots automatically generated and available on the HTML summary webpage, at:

<MG_ROOT_DIR>/MyCKKWLTest/index.html

where the plot pages is accessible by clicking on:

> Results and Event Database

leading you to a summary page of the various runs with links to pages showing specific aspects of the run. One of special interest is the 'DJR(plot)' link sending you to a page where all generated merging plots (with all uncertainties) are available in a table.

Let us know your comments on the above and how the runs went.

[lp:~maddevelopers/mg5amcnlo/2.3.3_PY8_install_cmd]

357. By Valentin Hirschi on 2016-02-25

1. Fixed a small issue in the plot creation of common_run_interface.py

358. By Olivier Mattelaer on 2016-02-26

fix the problem that I point for pythia/delphes command to select the type of run

359. By Olivier Mattelaer on 2016-02-26

adding help for the card edition

360. By Valentin Hirschi on 2016-02-27

1. Merged histograms.py with the large changes Rik did in this file for
   his branch improved pdf/scale reweighting.

361. By Valentin Hirschi on 2016-02-27

1. Updated default values for some options in the plotter.

362. By Valentin Hirschi on 2016-02-27

1. Updated the histograms.py IOTest

363. By Valentin Hirschi on 2016-03-03

1. Updated histograms.py and added JetMatching:ShowerKt parameter

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

One additional problem that I found:
I can not run any PY8 analysis with the PY8 version installed:
(This is annoying if I want to check the problem of the PY8 analysis which did not handle the
reweighting correctly)

Sorry for the annoying report

Olivier

Here is the detail of the crash:
INFO: Compiling MCatNLO for PYTHIA8...
WARNING: 'dl' was added to extralibs from the shower_card.dat.
It is needed for the correct running of PY8.2xx.
If this library cannot be found on your system, a crash will occur.
Undefined symbols for architecture x86_64:
  "_gzclose", referenced from:
      Pythia8::gzstreambuf::close() in libpythia8.a(Streams.o)
      Pythia8::gzstreambase::~gzstreambase() in libpythia8.a(Streams.o)
      Pythia8::gzstreambase::close() in libpythia8.a(Streams.o)
      Pythia8::gzstreambuf::~gzstreambuf() in libpythia8.a(Streams.o)
  "_gzopen", referenced from:
      Pythia8::gzstreambuf::open(char const*, int) in libpythia8.a(Streams.o)
      Pythia8::gzstreambase::open(char const*, int) in libPythia8 compilation did not succeed, exiting
from:
      Pythia8::gzstreambuf::underflow() in libpythia8.a(Streams.o)
  "_gzwrite", referenced from:
      Pythia8::gzstreambuf::flush_buffer() in libpythia8.a(Streams.o)
      Pythia8::gzstreambuf::overflow(int) in libpythia8.a(Streams.o)
      Pythia8::gzstreambuf::sync() in libpythia8.a(Streams.o)
      Pythia8::gzstreambuf::~gzstreambuf() in libpythia8.a(Streams.o)
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make: *** [Pythia82] Error 1

[Valentin Hirschi (valentin-hirschi)]

> Hi Valentin,
>
> One additional problem that I found:
> I can not run any PY8 analysis with the PY8 version installed:
> (This is annoying if I want to check the problem of the PY8 analysis which did
> not handle the
> reweighting correctly)
>

This is because the Pythia8 version we have is compiled so as to support zipped event files so that when compiling the driver, you need to add the 'zlib' library.
This is not a problem for our LO interface because the driver is compiled directly with the Pythia8 makefile which knowns what dependencies where turned on when it was installed.

Unfortunately, as you know the aMCatNLO interface to shower compiles its driver completely independently of Pythia8, so that you have to specify yourself what libraries to use when compiling the driver. You can therefore cure your problem by adding the 'libz' library to the list of external libraries to link in the shower_card, i.e.

EXTRALIBS = stdhep Fmcfio
->
EXTRALIBS = stdhep Fmcfio z

And with the change above it works (at least it does for me).

I don't know what would be the best way to have this be done automatically (and also for possible other dependencies of Pythia8). I don't think that we can do it in an elegant way within the aMC@NLO way of interfacing to shower, so for now it is best to leave it to the user to do it by hand.

> Sorry for the annoying report
>
> Olivier
>
> Here is the detail of the crash:
> INFO: Compiling MCatNLO for PYTHIA8...
> WARNING: 'dl' was added to extralibs from the shower_card.dat.
> It is needed for the correct running of PY8.2xx.
> If this library cannot be found on your system, a crash will occur.
> Undefined symbols for architecture x86_64:
> "_gzclose", referenced from:
> Pythia8::gzstreambuf::close() in libpythia8.a(Streams.o)
> Pythia8::gzstreambase::~gzstreambase() in libpythia8.a(Streams.o)
> Pythia8::gzstreambase::close() in libpythia8.a(Streams.o)
> Pythia8::gzstreambuf::~gzstreambuf() in libpythia8.a(Streams.o)
> "_gzopen", referenced from:
> Pythia8::gzstreambuf::open(char const*, int) in libpythia8.a(Streams.o)
> Pythia8::gzstreambase::open(char const*, int) in libPythia8 compilation
> did not succeed, exiting
> from:
> Pythia8::gzstreambuf::underflow() in libpythia8.a(Streams.o)
> "_gzwrite", referenced from:
> Pythia8::gzstreambuf::flush_buffer() in libpythia8.a(Streams.o)
> Pythia8::gzstreambuf::overflow(int) in libpythia8.a(Streams.o)
> Pythia8::gzstreambuf::sync() in libpythia8.a(Streams.o)
> Pythia8::gzstreambuf::~gzstreambuf() in libpythia8.a(Streams.o)
> ld: symbol(s) not found for architecture x86_64
> clang: error: linker command failed with exit code 1 (use -v to see
> invocation)
> make: *** [Pythia82] Error 1

[lp:~maddevelopers/mg5amcnlo/2.3.3_PY8_install_cmd]

364. By Valentin Hirschi on 2016-03-04

1. Changed the suggestion to install pythia8

[Valentin Hirschi (valentin-hirschi)]

> Hi Valentin,
>
> Very nice branch. It's really cool all the interesting way the code is
> improving those days.
> I particularly like the change in the card edition part (maybe I did those but
> I do not remember)

Well it is OUR common branch, which underlies our work for the PY8-MG5aMC enhanced automated interface+merging at LO.

> Here is my review:
> First some functionalities which does not work on this branch
>
> 1) The output pythia8 does not work in this branch. This need to be fixed for
> sure.
> - You change the comment for the PY8 path in the configuration file, is
> that related?

I don't think so, because the information on screen shows that it correctly substitutes the '.' with the path of the MG5 distribution. So this will need some investigation, but it seems that it is just that it tries to compile the 'main' file from within the root path of Pythia8 and not its correct subdirectory.
I tried to compile by hand the main file generated by MGaMC but it crashes for various reasons (mainly some include files missing).
Anyway, this command doesn't seem to work in 2.3.4 either and I have never looked at its implementation before, so I cannot (don't want to) fix it now.
But I agree that it would make sense to have it fixed for this version.

>
> 2) The standard plot routines (MA4) does not work anymore:
> - I do not have any plot either at parton level
> - I do not have any plot for after the shower (even if the lhe file in
> order to do such plot is created)

I installed MA4 and pythia-pgs and I could successfully run the full chain but indeed I don't get any MA4 plots.
This must be related to some checks or return statement in the function 'create_plot' of the common_run_interface. Nothing big a priori.
Notice that when running the shower 'by hand' with the command 'shower run_01' in the madevent interactive interface, the plots are created.

> 3) if I type the function shower, I have as question:
> Which programs do you want to run?
> 0 / auto : running existing cards
> 1 / pythia : Pythia
> 2 / pgs : Pythia + PGS
> 3 / delphes : Pythia + Delphes.
> [0, 1, 2, 3, auto, pythia, pgs, delphes][60s to answer]
> This need to be updated to the new format.
> At the end of the run, the code ask me a second question:
> Which programs do you want to run?
> 0 / auto : running existing cards
> 1 / pythia8 : Pythia8
> 3 / delphes : Pythia8 + Delphes.
> [0, 1, 3, auto, pythia8, delphes][60s to answer]
> (This probably happens only if you have both PY6 and PY8 installed)

Well the idea there is that if the user doesn't specify which shower it wants, i.e. with

shower pythia8 run_01

or

shower pythia run_01

then it does all of the available ones. I guess it should rather ask which one the user wants, but this is most efficiently done with the tab completetion.
So if you don't like the current behavior then I think that forcing the use to specify the shower is the simplest.

>
> One general question (having those answer will allow me to review that part of
> the modification)
>
> 1) I have a bunch of question related to the CKKW-L support that you have
> added in MG:
> a) Where can...

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

Here is some reply for the point who need them:
Thanks,

Olivier

> Well the idea there is that if the user doesn't specify which shower it wants, i.e. with
>
> shower pythia8 run_01
>
> or
>
> shower pythia run_01
>
> then it does all of the available ones. I guess it should rather ask which one the user wants, but this is most efficiently done with the tab completetion.
> So if you don't like the current behavior then I think that forcing the use to specify the shower is the simplest.

Clearly the current behavior is not satisfactory. Either having a question asking which one to run or force to specify is fine.
(or always choose PY8 is also fine) but like this no.

>> 1) I have a bunch of question related to the CKKW-L support that you have
>> added in MG:
>> a) Where can I find some documentation for how to run CKKW-L?
>
> On the pythia8 online manual:
> http://home.thep.lu.se/~torbjorn/pythia81html/Welcome.html
> (look for CKKW-L merging on the right hand side of the page)

We should create a launchpad FAQ and/or wiki page with minimal instruction (and refer to PY8)

> Again, please Stefan comment on that because I don't know much. The subtleties of the merging of decay products in CKKW-L are controlled by a combination of the two shower parameters ‘Merging:Process' and 'Merging:may_remove_decay_products'.

You need to have this clear since we need at least one of us knowing all the details of the MG side of the implementation.

> I agree that if the more general list support should be used for this as well, but it can be a bit tricky because you must make sure that corresponding fortran definition in run_card.inc is correctly written.

Indeed, but this is the same in Rik branch. It is important to use the same format. If you do not like the way Rik is writting this information it is then the perfect time to speak.

> On Mar 4, 2016, at 02:02, Valentin Hirschi <email address hidden> wrote:
>
>> Hi Valentin,
>>
>> Very nice branch. It's really cool all the interesting way the code is
>> improving those days.
>> I particularly like the change in the card edition part (maybe I did those but
>> I do not remember)
>
> Well it is OUR common branch, which underlies our work for the PY8-MG5aMC enhanced automated interface+merging at LO.
>
>> Here is my review:
>> First some functionalities which does not work on this branch
>>
>> 1) The output pythia8 does not work in this branch. This need to be fixed for
>> sure.
>> - You change the comment for the PY8 path in the configuration file, is
>> that related?
>
> I don't think so, because the information on screen shows that it correctly substitutes the '.' with the path of the MG5 distribution. So this will need some investigation, but it seems that it is just that it tries to compile the 'main' file from within the root path of Pythia8 and not its correct subdirectory.
> I tried to compile by hand the main file generated by MGaMC but it crashes for various reasons (mainly some include files missing).
> Anyway, this command doesn't seem to work in 2.3.4 either and I have never looked at its implementation bef...

[Rikkert Frederix (frederix)]

Hi Valentin,

I tried:

./bin/mg5_aMC
generate p p > e+ ve ; add process e+ ve j ; output ; launch
shower=PYTHIA8

and it crashes when starting showering with "KeyError : 'showerkt'" . Indeed, the showerkt parameter is not in the run_card.dat, while line 3333 of the madevent_interface.py requires it.

Cheers,
Rik

review: Needs Fixing

[lp:~maddevelopers/mg5amcnlo/2.3.3_PY8_install_cmd]

365. By Valentin Hirschi on 2016-03-31

1. Added the possibility of passing any option to the HEPToolsInstallers via the command install in MG5.
2. Fixed two bugs in the PY8-MG5aMC interface. One related to the new label structure in Histograms.py and one because showerkt is not in the run_card. This will need to be fetched elsewhere.

366. By Olivier Mattelaer on 2016-04-09

improve API for matplotlib output

367. By Valentin Hirschi on 2016-04-12

1. Fixed normalization issue in our .lhe by writing a function enforce_lha_strategy in lhe_parser.py.
2. Fixed a couple issues in the merging and properly defined 'qWeed'.

368. By Valentin Hirschi on 2016-04-12

1. Use a different lhe file for the enforcement of the lha_strategy -4.

369. By Valentin Hirschi on 2016-04-13

1. Fixed an issue with ptlund merging.

370. By Olivier Mattelaer on 2016-04-14

allow to generate event with event normalise to the average

371. By Olivier Mattelaer on 2016-04-14

update html for PY8 without matching

372. By Olivier Mattelaer on 2016-04-15

just make a warning for Syscalc but let him run (so far)

373. By Valentin Hirschi on 2016-04-18

1. Improved the MGPY8 interface by allowing to specify HEPMC output and
   also checking if PY8 crash.
2. Adding the pt observable to compare merging schemes whose DJR definition
   don't match.

374. By Valentin Hirschi on 2016-04-18

1. Revert the advanced_install hack to install local copies.
2. Added a Pythia8 option:
     HEPMCoutput:file
   Which can be set to /dev/null to avoid produce an lhe.

375. By Valentin Hirschi on 2016-04-19

1. Added the completion for the option specifying the path of the pythia8 or mg5amc tarball.

376. By Stefan Prestel on 2016-06-21

Fixed merging cut so that states with one massive particle are correctly handled. Also, incorrect do_cuts replaced by more suitable from_decay

377. By Stefan Prestel on 2016-06-21

Merged with latest revision, after updating merging cut.

378. By Valentin Hirschi on 2016-06-21

1. Added the missing option 'madanalysis5_path' in common_run_interface.py
2. Reinstated the warning about SysCalc not supporting the 'average' lhe convention.

379. By Stefan Prestel on 2016-06-21

Changed ktdurham to ensure that no ktdurham cut is applied to states with a no massless outgoing partons at all

380. By Valentin Hirschi on 2016-07-21

1. commented hack for local HEPToolsInstaller

381. By Valentin Hirschi on 2016-07-29

1. Added the completion for the new option '--no_root_in_MA5' for the advanced install function of MA5 (which vetoes the use of ROOT in MA5)

382. By Valentin Hirschi on 2016-08-11

1. Temporary disabling the generation of the MA5 hadron-level card since MA5 is broken as of now :(.
2. Added the central value to be part of the merging scale variations when those are automatically set by MG5aMC.

383. By Valentin Hirschi on 2016-08-11

1. Fixed an issue when setting the pythia8_path option in set_configuration.
2. Change all path settings in advanced_install to be absolute.
3. Remove the hack to pick a specific MG5 server in do_install().

384. By Valentin Hirschi on 2016-08-11

1. Insured that the central merging scale is added to the list of variations, unless the user directly specified this list in the PY8 card with the option 'SysCalc:qCutList' for MLM or 'SysCalc:tmsList' for CKKWL.

385. By Olivier Mattelaer on 2016-08-18

merge with 2.5.0

386. By Olivier Mattelaer on 2016-08-18

fix the output pythia8 command to make it compatible with the version of py8 installed via install PY8

387. By Olivier Mattelaer on 2016-08-18

fix the display of the options to put hidden the retro-compatible option like pythia

388. By Valentin Hirschi on 2016-08-18

1. First round of my fixes.

389. By Olivier Mattelaer on 2016-08-18

put back more decent printout for the refs

390. By Valentin Hirschi on 2016-08-19

1. Restructured the card edition

391. By Valentin Hirschi on 2016-08-19

1. Fixed typo

392. By Valentin Hirschi on 2016-08-19

1. Cleaned up the way Pythia8 cross-sections are stored.

393. By Valentin Hirschi on 2016-08-19

1. Added '-lstdc++' or '-lc++' to the MCatNLO compilation for Pythia8 if (and only if)
   these were specified when compiling Pythia8.
2. Fixed the specification of the _cpp_dependencies to HEPToolsInstaller when the user uses
   clang and didn't specify his compiler in the MG5 options.

[Olivier Mattelaer (olivier-mattelaer)]

This branch is now out of date. So I close the merge request.
As far as I know Valentin (or I) implemented all the comment made in this review.
The code obviously re-need a full review.

review: Disapprove

Unmerged revisions

393. By Valentin Hirschi on 2016-08-19

1. Added '-lstdc++' or '-lc++' to the MCatNLO compilation for Pythia8 if (and only if)
   these were specified when compiling Pythia8.
2. Fixed the specification of the _cpp_dependencies to HEPToolsInstaller when the user uses
   clang and didn't specify his compiler in the MG5 options.

392. By Valentin Hirschi on 2016-08-19

1. Cleaned up the way Pythia8 cross-sections are stored.

391. By Valentin Hirschi on 2016-08-19

1. Fixed typo

390. By Valentin Hirschi on 2016-08-19

1. Restructured the card edition

389. By Olivier Mattelaer on 2016-08-18

put back more decent printout for the refs

388. By Valentin Hirschi on 2016-08-18

1. First round of my fixes.

387. By Olivier Mattelaer on 2016-08-18

fix the display of the options to put hidden the retro-compatible option like pythia

386. By Olivier Mattelaer on 2016-08-18

fix the output pythia8 command to make it compatible with the version of py8 installed via install PY8

385. By Olivier Mattelaer on 2016-08-18

merge with 2.5.0

384. By Valentin Hirschi on 2016-08-11

1. Insured that the central merging scale is added to the list of variations, unless the user directly specified this list in the PY8 card with the option 'SysCalc:qCutList' for MLM or 'SysCalc:tmsList' for CKKWL.

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2016-03-18 13:28:29 +0000
+++ MadSpin/decay.py	2016-08-19 23:07:19 +0000
@@ -1224,6 +1224,7 @@
                                   'P%s' % me.get('processes')[0].shell_string())
             topo['decays'] = []
             topo['total_br'] = 0 
+            topo['Pid'] = proc.get('id')
         
             if skip:
                 self.add_me_symmetric(skip, topo)
@@ -1636,7 +1637,7 @@
                 (' '.join([str(i) for i in opts['particles']]),
                  ' '.join('--%s=%s' % (key,value) for (key,value) in opts.items()
                         if key not in ['model', 'force', 'particles'] and value))
-        cmd.exec_cmd('import model %s' % model.get('name'))
+        cmd.exec_cmd('import model %s' % model.get('modelpath+restriction'))
 #        cmd._curr_model = model
 #        cmd._curr_fortran_model = helas_call_writers.FortranUFOHelasCallWriter(model)
         cmd.exec_cmd(line)
@@ -2115,7 +2116,8 @@
         # launch the decay and reweighting
         self.mscmd.update_status('MadSpin: Decaying Events')
         efficiency = self.decaying_events(self.inverted_decay_mapping)
-        if  efficiency != 1:
+        self.efficiency = efficiency
+        if  efficiency != 1 and any(v==-1 for v in self.br_per_id.values()):
             # need to change the banner information [nb_event/cross section]
             files.cp(self.outputfile.name, '%s_tmp' % self.outputfile.name)
             self.outputfile = open(self.outputfile.name, 'w')
@@ -2612,7 +2614,7 @@
         return decay_mapping
     
 
-    #@misc.mute_logger()
+    @misc.mute_logger()
     @test_aloha.set_global()
     def generate_all_matrix_element(self):
         """generate the full series of matrix element needed by Madspin.
@@ -2651,13 +2653,12 @@
                       if re.search('^\s*add\s+process', line)]
         
         mgcmd = self.mgcmd
-        modelpath = self.model.get('modelpath')
-        if os.path.basename(modelpath) != mgcmd._curr_model['name']:
-            name, restrict = mgcmd._curr_model['name'].rsplit('-',1)
-            if os.path.exists(pjoin(os.path.dirname(modelpath),name, 'restrict_%s.dat' % restrict)):
-                modelpath = pjoin(os.path.dirname(modelpath), mgcmd._curr_model['name'])
+        modelpath = self.model.get('modelpath+restriction')
+
+        commandline="import model %s" % modelpath
+        if not self.model.mg5_name:
+            commandline += ' --modelname'
             
-        commandline="import model %s " % modelpath
         mgcmd.exec_cmd(commandline)
         # Handle the multiparticle of the banner        
         #for name, definition in self.mscmd.multiparticles:
@@ -2685,13 +2686,13 @@
                 proc_nb = '@ %i' % proc_nb 
             else:
                 proc_nb = ''
-            
+                        
             if ',' in proc:
                 raise MadSpinError, 'MadSpin can not decay event which comes from a decay chain.'+\
                         '\n  The full decay chain should either be handle by MadGraph or by Masdspin.'
             
             if '[' not in proc:
-                commandline+="add process %s  --no_warning=duplicate;" % proc
+                commandline+="add process %s %s  --no_warning=duplicate;" % (proc, proc_nb)
             else:
                 process, order, final = re.split('\[\s*(.*)\s*\]', proc)
                 commandline+="add process %s %s --no_warning=duplicate;" % (process, proc_nb)
@@ -2723,6 +2724,7 @@
                         commandline +='add process %s pert_%s  %s --no_warning=duplicate;' % (process,order, proc_nb)                                       
         commandline = commandline.replace('add process', 'generate',1)
         logger.info(commandline)
+        
         mgcmd.exec_cmd(commandline, precmd=True)
         commandline = 'output standalone_msP %s %s' % \
         (pjoin(path_me,'production_me'), ' '.join(self.list_branches.keys()))        
@@ -3110,7 +3112,7 @@
 
         
         probe_weight = []
-        
+       
 
         starttime = time.time()
         ev = -1
@@ -3148,6 +3150,7 @@
             atleastonedecay=False
             for decay in self.all_ME[production_tag]['decays']:
                 tag = decay['decay_tag']
+
                 if decay_mapping and not tag in decay_mapping:
                     continue
                 if not tag:
@@ -3161,10 +3164,11 @@
                 else:
                     max_decay[tag] = weight
                     #print weight, max_decay[name]
-                    #raise Exception   
+                    #raise Exception 
+                      
             if not atleastonedecay:
                 # NO decay [one possibility is all decay are identical to their particle]
-                logger.info('No independent decay for one type of final states -> skip those events')
+                logger.info('No independent decay for one type of final states -> skip those events for the maximum weight computation')
                 nb_decay[decaying] = numberev
                 ev += numberev -1
                 continue
@@ -3295,7 +3299,7 @@
 
         return max_weight
         
-    def loadfortran(self, mode, path, stdin_text):
+    def loadfortran(self, mode, path, stdin_text, first=True):
         """ call the fortran executable """
 
         tmpdir = ''
@@ -3312,8 +3316,18 @@
             self.calculator[('full',path,)] = external 
             self.calculator_nbcall[('full',path)] = 1 
 
-        external.stdin.write(stdin_text)
-        
+        try:
+            external.stdin.write(stdin_text)
+        except IOError:
+            if not first:
+                raise
+            try:
+                external.terminate()
+            except:
+                pass
+            del self.calculator[('full',path,)]
+            return self.loadfortran(mode, path, stdin_text, first=False)
+
         if mode == 'maxweight':
             maxweight=float(external.stdout.readline())
             output = maxweight
@@ -3864,16 +3878,24 @@
         ms_banner = ""
         cross_section = True # tell if possible to write the cross-section in advance
         total_br = []
+        self.br_per_id = {}
         for production in self.all_ME.values():
             one_br = 0
+            partial_br = 0
             for decay in production['decays']:
                 if not decay['decay_tag']:
                     cross_section = False
                     one_br += decay['br']
                     continue
+                partial_br += decay['br']
                 ms_banner += "# %s\n" % ','.join(decay['decay_tag']).replace('\n',' ')
                 ms_banner += "# BR: %s\n# max_weight: %s\n" % (decay['br'], decay['max_weight'])
                 one_br += decay['br']
+            
+            if production['Pid'] not in self.br_per_id:
+                self.br_per_id[production['Pid']] = partial_br
+            elif self.br_per_id[production['Pid']] != partial_br:
+                self.br_per_id[production['Pid']] = -1
             total_br.append(one_br)
         
         if __debug__:
@@ -3881,7 +3903,6 @@
                 assert production['total_br'] - min(total_br) < 1e-4
         
         self.branching_ratio = max(total_br) * eff
-
         #self.banner['madspin'] += ms_banner
         # Update cross-section in the banner
         if 'mggenerationinfo' in self.banner:
@@ -3912,18 +3933,36 @@
                     mg_info[i] = '%s : %s' % (info, value * self.branching_ratio)
                 self.banner['mggenerationinfo'] = '\n'.join(mg_info)
                 
-            
-        
-        if 'init' in self.banner:
+        self.cross = 0
+        self.error = 0
+        if 'init' in self.banner and (eff!=1 or not any(v==-1 for v in self.br_per_id.values())):
             new_init =''
+            curr_proc = 0
+            has_missing=False
             for line in self.banner['init'].split('\n'):
                 if len(line.split()) != 4:
                     new_init += '%s\n' % line
                 else:
+                    curr_proc += 1 
                     data = [float(nb) for nb in line.split()]
-                    data[:3] = [ data[i] * self.branching_ratio for i  in range(3)]
+                    id = int(data[-1])
+                    if id in self.br_per_id and not any(v==-1 for v in self.br_per_id.values()):
+                        data[:3] = [data[i] * self.br_per_id[id] for i  in range(3)]
+                    else:
+                        data[:3] = [ data[i] * self.branching_ratio for i  in range(3)]
+                        has_missing=True
                     new_init += ' %.12E %.12E %.12E %i\n' % tuple(data)
+                    cross, error = [float(d) for d in data[:2]]
+                    self.cross += cross
+                    self.error += error**2
+                    
+                    
             self.banner['init'] = new_init
+            self.error = math.sqrt(self.error)
+            if has_missing and curr_proc not in [0,1]:
+                logger.warning('''The partial cross section for each subprocess can not be determine. due
+    Reason: multiple final state in the same subprocess (and the presence of multiple BR)
+    Consequence: the <init> information of the lhe will therefore be incorrect. Please correct it if needed.''')
         self.banner.write(self.outputfile, close_tag=False)        
         
     def terminate_fortran_executables(self, path_to_decay=0 ):

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2016-03-18 13:28:29 +0000
+++ MadSpin/interface_madspin.py	2016-08-19 23:07:19 +0000
@@ -209,10 +209,10 @@
             mg_names = True
         model_name = self.banner.get('proc_card', 'model')
         if model_name:
+            model_name = os.path.expanduser(model_name)
             self.load_model(model_name, mg_names, complex_mass)
         else:
             raise self.InvalidCmd('Only UFO model can be loaded in MadSpin.')
-        
         # check particle which can be decayed:
         self.final_state = set()
         final_model = False
@@ -243,7 +243,7 @@
     def import_model(self, args):
         """syntax: import model NAME CARD_PATH
             args didn't include import model"""
-            
+        
         bypass_check = False
         if '--bypass_check' in args:
             args.remove('--bypass_check')
@@ -392,6 +392,9 @@
         args = self.split_arg(line)
         self.check_set(args)
         
+        if args[0] not in ['ms_dir', 'run_card']:
+            args[1] = args[1].lower()
+        
         if args[0] in  ['max_weight', 'BW_effect','ms_dir', 'spinmode']:
             self.options[args[0]] = args[1]
             if args[0] == 'ms_dir':
@@ -460,11 +463,11 @@
         """ """
 
         try:
-            self.mg5cmd.do_define(line)
+            self.mg5cmd.exec_cmd('define %s' % line)
         except:
             #cleaning if the error is recover later
             key = line.split()[0]
-            if hasattr(self, 'multiparticles_ms' and key in self.multiparticles_ms):
+            if hasattr(self, 'multiparticles_ms') and key in self.multiparticles_ms:
                 del self.multiparticles_ms[key]
             raise
            
@@ -566,6 +569,9 @@
         generate_all.run()
                         
         self.branching_ratio = generate_all.branching_ratio
+        self.cross = generate_all.cross
+        self.error = generate_all.error
+        self.efficiency = generate_all.efficiency
         try:
             self.err_branching_ratio = generate_all.err_branching_ratio
         except Exception:
@@ -662,11 +668,15 @@
             #seed is specified need to use that one:
             open(pjoin(self.options['ms_dir'],'seeds.dat'),'w').write('%s\n'%self.seed)
             #remove all ranmar_state
-            for name in glob.glob(pjoin(self.options['ms_dir'], '*', 'SubProcesses','*','ranmar_state.dat')):
+            for name in misc.glob(pjoin('*', 'SubProcesses','*','ranmar_state.dat'), 
+                                                        self.options['ms_dir']):
                 os.remove(name)    
         
         generate_all.ending_run()
         self.branching_ratio = generate_all.branching_ratio
+        self.cross = generate_all.cross
+        self.error = generate_all.error
+        self.efficiency = generate_all.efficiency
         try:
             self.err_branching_ratio = generate_all.err_branching_ratio
         except Exception:
@@ -691,8 +701,7 @@
         #   of event to generate for each type of particle.
         # 2. Generate the events requested
         # 3. perform the merge of the events.
-        #    if not enough events. re-generate the missing one.
-                
+        #    if not enough events. re-generate the missing one.  
         # First define an utility function for generating events when needed
         def generate_events(pdg, nb_event, mg5, restrict_file=None, cumul=False):
             """generate new events for this particle
@@ -786,7 +795,7 @@
                 os.mkdir(self.path_me) 
         else:
             # cleaning
-            for name in glob.glob(pjoin(self.path_me, "decay_*_*")):
+            for name in misc.glob("decay_*_*", self.path_me):
                 shutil.rmtree(name)
 
         self.events_file.close()
@@ -821,8 +830,9 @@
         # 2. Generate the events requested
         with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
             mg5 = self.mg5cmd
-            modelpath = self.model.get('modelpath+restriction')
-            mg5.exec_cmd("import model %s" % modelpath)      
+            if not self.model:
+                modelpath = self.model.get('modelpath+restriction')
+                mg5.exec_cmd("import model %s" % modelpath)      
             to_event = {}
             for pdg, nb_needed in to_decay.items():
                 #check if a splitting is needed
@@ -880,6 +890,12 @@
             else:
                 raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
         self.branching_ratio = br
+        self.efficiency = 1
+        self.cross, self.error = self.banner.get_cross(witherror=True)
+        self.cross *= br
+        self.error *= br
+        
+        
         # modify the cross-section in the init block of the banner
         self.banner.scale_init_cross(self.branching_ratio)
                     
@@ -982,7 +998,6 @@
     def load_model(self, name, use_mg_default, complex_mass=False):
         """load the model"""
         
-        
         loop = False
         #if (name.startswith('loop_')):
         #    logger.info("The model in the banner is %s" % name)

=== added file 'Template/LO/Cards/pythia8_card_default.dat'
--- Template/LO/Cards/pythia8_card_default.dat	1970-01-01 00:00:00 +0000
+++ Template/LO/Cards/pythia8_card_default.dat	2016-08-19 23:07:19 +0000
@@ -0,0 +1,62 @@
+!
+! Pythia8 cmd card automatically generated by MadGraph5_aMC@NLO
+!
+! ==================
+! General parameters 
+! ==================
+!
+Main:numberOfEvents      = -1
+!
+! -------------------------------------------------------------------
+! Specify the HEPMC output of the Pythia8 shower. You can set it to:
+!   auto  : to let MG5aMC automatically assign it
+!   fifo  : to have MG5aMC setup the piping of the PY8 output to analysis
+!           tools such as MadAnalysis5. You can also use the syntax 
+!           fifo@<fifo_path> to select where the fifo file should be
+!           created (useful to select a mounted drive that supports fifo)
+!   /dev/null : to turn off the HEPMC output.
+!   <path>    : to select where the HEPMC file must written. It will 
+!               therefore *not* be available in the process output.
+! -------------------------------------------------------------------
+!
+HEPMCoutput:file         = auto
+!
+! --------------------------------------------------------------------
+! Parameters relevant only when performing MLM merging
+! For details, see section 'Jet Matching' on the left-hand menu of 
+!    http://home.thep.lu.se/~torbjorn/pythia81html/Welcome.html
+! --------------------------------------------------------------------
+! If equal to -1.0, MadGraph5_aMC@NLO will set it automatically based 
+! on the parameter 'xqcut' of the run_card.dat 
+JetMatching:qCut         = -1.0
+! Use default kt-MLM to match parton level jets to those produced by the
+! shower. But the other Shower-kt scheme is available too with this option.
+JetMatching:doShowerKt   = off
+! A value of -1 means that it is automatically guessed by MadGraph.
+! It is however always safer to explicitly set it.
+JetMatching:nJetMax      = -1
+!
+! --------------------------------------------------------------------
+! Parameters relevant only when performing CKKWL merging
+! For details, see section 'CKKW-L Merging' on the left-hand menu of 
+!    http://home.thep.lu.se/~torbjorn/pythia81html/Welcome.html
+! --------------------------------------------------------------------
+! Central merging scale values you want to be used.
+! If equal to -1.0, then MadGraph5_aMC@NLO will set this automatically
+! based on the parameter 'ktdurham' of the run_card.dat
+Merging:TMS              = -1.0
+! This must be set manually, according to Pythia8 directives.
+! An example of possible value is 'pp>LEPTONS,NEUTRINOS'
+Merging:Process          = <set_by_user>
+! A value of -1 means that it is automatically guessed by MadGraph.
+! It is however always safer to explicitly set it.
+Merging:nJetMax  		 = -1
+!
+! For all merging schemes, decide wehter you want the merging scale
+! variation computed for only the central weights or all other 
+! PDF and scale variation weights as well
+SysCalc:fullCutVariation = off
+!
+! ==========================
+! User customized parameters 
+! ==========================

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2015-08-11 00:35:53 +0000
+++ Template/LO/Cards/run_card.dat	2016-08-19 23:07:19 +0000
@@ -21,10 +21,6 @@
 #*********************************************************************
   %(run_tag)s     = run_tag ! name of the run 
 #*********************************************************************
-# Run to generate the grid pack                                      *
-#*********************************************************************
-  %(gridpack)s     = gridpack  !True = setting up the grid pack
-#*********************************************************************
 # Number of events and rnd seed                                      *
 # Warning: Do not generate more than 1M events in a single run       *
 # If you want to run Pythia, avoid more than 50k events in a run.    *
@@ -61,51 +57,52 @@
  %(dynamical_scale_choice)s = dynamical_scale_choice ! Choose one of the preselected dynamical choices
  %(scalefact)s  = scalefact        ! scale factor for event-by-event scales
 #*********************************************************************
-# Time of flight information. (-1 means not run)
-#*********************************************************************
- %(time_of_flight)s = time_of_flight ! threshold below which info is not written
-#*********************************************************************
-# Matching - Warning! ickkw > 1 is still beta
-#*********************************************************************
- %(ickkw)s = ickkw            ! 0 no matching, 1 MLM, 2 CKKW matching
- %(highestmult)s = highestmult      ! for ickkw=2, highest mult group
- %(ktscheme)s = ktscheme         ! for ickkw=1, 1 Durham kT, 2 Pythia pTE
+# Type and output format
+#*********************************************************************
+  %(gridpack)s     = gridpack  !True = setting up the grid pack
+  %(time_of_flight)s = time_of_flight ! threshold below which info is not written (-1 means not written)
+  %(lhe_version)s = lhe_version       ! Change the way clustering information pass to shower.        
+  %(clusinfo)s = clusinfo         ! include clustering tag in output
+  %(event_norm)s =  event_norm       ! average/sum. Normalization of the weight in the LHEF
+
+#*********************************************************************
+# Matching parameter (MLM only)
+#*********************************************************************
+ %(ickkw)s = ickkw            ! 0 no matching, 1 MLM
  %(alpsfact)s = alpsfact         ! scale factor for QCD emission vx
  %(chcluster)s = chcluster        ! cluster only according to channel diag
- %(pdfwgt)s = pdfwgt           ! for ickkw=1, perform pdf reweighting
  %(asrwgtflavor)s = asrwgtflavor     ! highest quark flavor for a_s reweight
- %(clusinfo)s = clusinfo         ! include clustering tag in output
- %(lhe_version)s = lhe_version       ! Change the way clustering information pass to shower.        
+ %(auto_ptj_mjj)s  = auto_ptj_mjj  ! Automatic setting of ptj and mjj if xqcut >0
+                                   ! (turn off for VBF and single top processes) 
+ %(xqcut)s   = xqcut   ! minimum kt jet measure between partons
 #*********************************************************************
-#**********************************************************
 #
-#**********************************************************
-# Automatic ptj and mjj cuts if xqcut > 0
-# (turn off for VBF and single top processes)
-#**********************************************************
-   %(auto_ptj_mjj)s  = auto_ptj_mjj  ! Automatic setting of ptj and mjj
-#**********************************************************
-#                                                                    
-#**********************************
-# BW cutoff (M+/-bwcutoff*Gamma)
-#**********************************
-  %(bwcutoff)s  = bwcutoff      ! (M+/-bwcutoff*Gamma)
-#**********************************************************
-# Apply pt/E/eta/dr/mij/kt_durham cuts on decay products or not
-# (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
-#*************************************************************
-   %(cut_decays)s  = cut_decays    ! Cut decay products 
-#*************************************************************
+#*********************************************************************
 # Number of helicities to sum per event (0 = all helicities)
 # 0 gives more stable result, but longer run time (needed for
 # long decay chains e.g.).
 # Use >=2 if most helicities contribute, e.g. pure QCD.
-#*************************************************************
+#*********************************************************************
    %(nhel)s  = nhel          ! Number of helicities used per event
-#*******************                                                 
-# Standard Cuts
-#*******************                                                 
+#*********************************************************************
+#
+#
+#*******************************                                                 
+# Parton level cuts definition *
+#*******************************                                     
 #                                                                    
+#
+#*********************************************************************
+# BW cutoff (M+/-bwcutoff*Gamma) ! Define on/off-shell for "$" and decay  
+#*********************************************************************
+  %(bwcutoff)s  = bwcutoff      ! (M+/-bwcutoff*Gamma)
+#*********************************************************************
+# Apply pt/E/eta/dr/mij/kt_durham cuts on decay products or not
+# (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
+#*********************************************************************
+   %(cut_decays)s  = cut_decays    ! Cut decay products 
+#*********************************************************************
+# Standard Cuts                                                      *
 #*********************************************************************
 # Minimum and maximum pt's (for max, -1 means no cut)                *
 #*********************************************************************
@@ -246,21 +243,20 @@
 #*********************************************************************
  %(xetamin)s   = xetamin ! minimum rapidity for two jets in the WBF case  
  %(deltaeta)s   = deltaeta ! minimum rapidity for two jets in the WBF case 
-#*********************************************************************
-# KT DURHAM CUT                                                      *
-#*********************************************************************
- %(ktdurham)s    =  ktdurham        
- %(dparameter)s  =  dparameter 
+#***********************************************************************
+# Turn on either the ktdurham or ptlund cut to activate                *
+# CKKW(L) merging with Pythia8 [arXiv:1410.3012, arXiv:1109.4829]      *
+#***********************************************************************
+ %(ktdurham)s  =  ktdurham        
+ %(dparameter)s   =  dparameter
+ %(ptlund)s  =  ptlund
+ %(pdgs_for_merging_cut)s  =  pdgs_for_merging_cut ! PDGs for two cuts above   
 #*********************************************************************
 # maximal pdg code for quark to be considered as a light jet         *
 # (otherwise b cuts are applied)                                     *
 #*********************************************************************
  %(maxjetflavor)s = maxjetflavor    ! Maximum jet pdg code
 #*********************************************************************
-# Jet measure cuts                                                   *
-#*********************************************************************
- %(xqcut)s   = xqcut   ! minimum kt jet measure between partons
-#*********************************************************************
 #
 #*********************************************************************
 # Store info for systematics studies                                 *
@@ -281,3 +277,4 @@
 # PDF sets and number of members (0 or none for all members).
 %(sys_pdf)s = sys_pdf # matching scales
 # MSTW2008nlo68cl.LHgrid 1  = sys_pdf
+#

=== modified file 'Template/LO/Source/PDF/makefile'
--- Template/LO/Source/PDF/makefile	2015-11-03 13:40:56 +0000
+++ Template/LO/Source/PDF/makefile	2016-08-19 23:07:19 +0000
@@ -12,9 +12,9 @@
 
 ifdef lhapdf
   ifeq ($(lhapdfversion),5)
-    PDF         = pdfwrap_lhapdf.o pdf_lhapdf.o pdg2pdf_lhapdf.o opendata.o
+    PDF         = pdfwrap_lhapdf.o pdf_lhapdf.o pdg2pdf_lhapdf.o opendata.o  PhotonFlux.o
   else
-    PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o
+    PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o
   endif
 else
   PDF         = Ctq4Fn.o Ctq5Par.o Ctq5Pdf.o Partonx5.o Ctq6Pdf.o cteq3.o \

=== modified file 'Template/LO/Source/PDF/pdf_lhapdf.f'
--- Template/LO/Source/PDF/pdf_lhapdf.f	2014-06-19 11:20:33 +0000
+++ Template/LO/Source/PDF/pdf_lhapdf.f	2016-08-19 23:07:19 +0000
@@ -11,6 +11,7 @@
       INTEGER IH,I
       double precision  photon
       LOGICAL has_photon
+      double precision epa_electron,epa_proton
 C
 C     Include
 C
@@ -28,11 +29,11 @@
          do i=-6,6
             pdf(i)=f(i)/x
          enddo
+      elseif(ih .eq. 2) then  !from a proton without breaking
+          pdf(7) = epa_proton(x, q * q)
       else
          write (*,*) 'beam type not supported in lhadpf'
-         do i=-6,6
-            pdf(i)=0d0
-         enddo
+         stop 1
       endif
 
       return	

=== modified file 'Template/LO/Source/PDF/pdg2pdf_lhapdf6.f'
--- Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2016-02-02 17:38:59 +0000
+++ Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2016-08-19 23:07:19 +0000
@@ -16,6 +16,7 @@
       integer i,j,ihlast(20),ipart,iporg,ireuse,imemlast(20),iset,imem
      &     ,i_replace,ii,ipartlast(20)
       double precision xlast(20),xmulast(20),pdflast(20)
+      double precision epa_proton
       save ihlast,xlast,xmulast,pdflast,imemlast,ipartlast
       data ihlast/20*-99/
       data ipartlast/20*-99/
@@ -99,8 +100,15 @@
 
 c     Call lhapdf and give the current values to the arrays that should
 c     be saved
-      call evolvepart(ipart,x,xmu,pdg2pdf)
-      pdg2pdf=pdg2pdf/x
+      if(ih.eq.1) then
+         call evolvepart(ipart,x,xmu,pdg2pdf)
+         pdg2pdf=pdg2pdf/x
+      else if(ih.eq.2) then ! photon from a proton without breaking
+          pdg2pdf = epa_proton(x,xmu*xmu)
+      else
+         write (*,*) 'beam type not supported in lhadpf'
+         stop 1
+      endif
       pdflast(i_replace)=pdg2pdf
       xlast(i_replace)=x
       xmulast(i_replace)=xmu

=== modified file 'Template/LO/Source/cuts.inc'
--- Template/LO/Source/cuts.inc	2015-04-09 01:31:09 +0000
+++ Template/LO/Source/cuts.inc	2016-08-19 23:07:19 +0000
@@ -88,6 +88,10 @@
       REAL*8 KT_DURHAM, D_PARAMETER
       LOGICAL DO_KT_DURHAM
       COMMON /JET_MEASURE_CUTS/ KT_DURHAM, D_PARAMETER
-      
-      
-      
\ No newline at end of file
+
+C COMMON BLOCK FOR PYTHIA JET MEASURE CUT
+      REAL*8 PT_LUND
+      COMMON /PYTHIA_MEASURE_CUTS/ PT_LUND
+      
+      
+      

=== modified file 'Template/LO/Source/make_opts'
--- Template/LO/Source/make_opts	2016-05-06 09:17:15 +0000
+++ Template/LO/Source/make_opts	2016-08-19 23:07:19 +0000
@@ -1,9 +1,9 @@
 DEFAULT_F2PY_COMPILER=f2py
 DEFAULT_F_COMPILER=gfortran
-MACFLAG=-mmacosx-version-min=10.7
-DEFAULT_CPP_COMPILER=clang
-STDLIB=-lc++
-STDLIB_FLAG=-stdlib=libc++
+MACFLAG=
+DEFAULT_CPP_COMPILER=g++
+STDLIB=-lstdc++
+STDLIB_FLAG=
 #end_of_make_opts_variables
 # Rest of the makefile
 
@@ -37,10 +37,10 @@
 
 # Set FC unless it's defined by an environment variable
 ifeq ($(origin FC),default)
-FC=gfortran
+FC=$(DEFAULT_F_COMPILER)
 endif
 ifeq ($(origin F2PY), undefined)
-F2PY=f2py
+F2PY=$(DEFAULT_F2PY_COMPILER)
 endif
 
 UNAME := $(shell uname -s)
@@ -83,7 +83,7 @@
 
 # Option lhapdf
 
-ifdef lhapdf
+ifneq ($(lhapdf),)
 CXXFLAGS += $(shell $(lhapdf) --cppflags)
 alfas_functions=alfas_functions_lhapdf
 llhapdf+= -lLHAPDF

=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2015-04-02 22:56:24 +0000
+++ Template/LO/Source/run.inc	2016-08-19 23:07:19 +0000
@@ -63,3 +63,11 @@
 C
       logical MC_grouped_subproc
       common/to_MC_grouped_subproc/MC_grouped_subproc
+
+C
+C Controls what are the PDGs included in the CKKWl merging procedure, i.e. what
+C are the PDGs subject to the ktdurham cut
+C
+      integer pdgs_for_merging_cut(0:1000)
+      common/TO_MERGE/pdgs_for_merging_cut
+

=== modified file 'Template/LO/SubProcesses/cuts.f'
--- Template/LO/SubProcesses/cuts.f	2015-08-11 00:35:53 +0000
+++ Template/LO/SubProcesses/cuts.f	2016-08-19 23:07:19 +0000
@@ -1,1317 +1,1750 @@
-      logical function pass_point(p)
-c************************************************************************
-c     This function is called from sample to see if it needs to 
-c     bother calculating the weight from all the different conficurations
-c     You can either just return true, or have it call passcuts
-c************************************************************************
-      implicit none
-c
-c     Arguments
-c
-      double precision p
-c
-c     External
-c
-      logical passcuts
-      external passcuts
-c-----
-c  Begin Code
-c-----
-      pass_point = .true.
-c      pass_point = passcuts(p)
-      end
-C 
-      LOGICAL FUNCTION PASSCUTS(P)
-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 ...
-C            ALL MOMENTA ARE IN THE REST FRAME!!
-C            COMMON/JETCUTS/   CUTS ON JETS
-C     OUTPUT:
-C            TRUE IF EVENTS PASSES ALL CUTS LISTED
-C**************************************************************************
-      IMPLICIT NONE
-c
-c     Constants
-c
-      include 'genps.inc'
-      include 'nexternal.inc'
-C
-C     ARGUMENTS
-C
-      REAL*8 P(0:3,nexternal)
-
-C
-C     LOCAL
-C
-      LOGICAL FIRSTTIME,FIRSTTIME2,pass_bw,notgood,good,foundheavy
-      LOGICAL DEBUG
-      integer i,j,njets,nheavyjets,nleptons,hardj1,hardj2
-      REAL*8 XVAR,ptmax1,ptmax2,htj,tmp,inclht
-      real*8 ptemp(0:3), ptemp2(0:3)
-      character*20 formstr
-C
-C     PARAMETERS
-C
-      real*8 PI
-      parameter( PI = 3.14159265358979323846d0 )
-C
-C     EXTERNAL
-C
-      REAL*8 R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,PtDot
-      logical cut_bw,setclscales
-      external R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,cut_bw,setclscales,PtDot
-C
-C     GLOBAL
-C
-      include 'run.inc'
-      include 'cuts.inc'
-
-      double precision ptjet(nexternal)
-      double precision ptheavyjet(nexternal)
-      double precision ptlepton(nexternal)
-      double precision temp
-
-C VARIABLES TO SPECIFY JETS
-      DOUBLE PRECISION PJET(NEXTERNAL,0:3)
-      DOUBLE PRECISION PTMIN
-      DOUBLE PRECISION PT1,PT2
-      INTEGER K,J1,J2
-
-C VARIABLES FOR KT CUT
-      DOUBLE PRECISION PTNOW,COSTH,PABS1,PABS2
-      DOUBLE PRECISION ETA1,ETA2,COSH_DETA,COS_DPHI,KT1SQ,KT2SQ, DPHI
-
-      double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
-      double precision emin(nincoming+1:nexternal)
-      double precision                    r2min(nincoming+1:nexternal,nincoming+1:nexternal)
-      double precision s_min(nexternal,nexternal)
-      double precision etmax(nincoming+1:nexternal),etamin(nincoming+1:nexternal)
-      double precision emax(nincoming+1:nexternal)
-      double precision r2max(nincoming+1:nexternal,nincoming+1:nexternal)
-      double precision s_max(nexternal,nexternal)
-      double precision ptll_min(nexternal,nexternal),ptll_max(nexternal,nexternal)
-      double precision inclHtmin,inclHtmax
-      common/to_cuts/  etmin, emin, etamax, r2min, s_min,
-     $     etmax, emax, etamin, r2max, s_max, ptll_min, ptll_max, inclHtmin,inclHtmax
-
-      double precision ptjmin4(4),ptjmax4(4),htjmin4(2:4),htjmax4(2:4)
-      logical jetor
-      common/to_jet_cuts/ ptjmin4,ptjmax4,htjmin4,htjmax4,jetor
-
-      double precision ptlmin4(4),ptlmax4(4)
-      common/to_lepton_cuts/ ptlmin4,ptlmax4
-
-c
-c     Special cuts
-c
-
-      integer        lbw(0:nexternal)  !Use of B.W.
-      common /to_BW/ lbw
-C
-C     SPECIAL CUTS
-C
-      LOGICAL  IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
-      LOGICAL  IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
-      LOGICAL  IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
-      logical  do_cuts(nexternal)
-      COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
-     . IS_A_ONIUM, do_cuts
-C
-C     Keep track of whether cuts already calculated for this event
-C
-      LOGICAL CUTSDONE,CUTSPASSED
-      COMMON/TO_CUTSDONE/CUTSDONE,CUTSPASSED
-      DATA CUTSDONE,CUTSPASSED/.FALSE.,.FALSE./
-
-C $B$ MW_NEW_DEF $E$ !this is a tag for MadWeight
-
-      double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
-      common/to_xqcuts/xqcutij,xqcuti
-
-c jet cluster algorithm
-      integer nQCD !,NJET,JET(nexternal)
-c      double precision plab(0:3, nexternal)
-      double precision pQCD(0:3,nexternal)!,PJET(0:3,nexternal)
-c      double precision rfj,sycut,palg,fastjetdmerge
-c      integer njet_eta
-c     Photon isolation
-      integer nph,nem,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)
-      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)
-
-      include 'coupl.inc'
-C
-C
-c
-      DATA FIRSTTIME,FIRSTTIME2/.TRUE.,.TRUE./
-
-c put momenta in common block for couplings.f
-      double precision pp(0:3,max_particles)
-      common /momenta_pp/pp
-
-      DATA DEBUG/.FALSE./
-
-C-----
-C  BEGIN CODE
-C-----
-
-
-
-      PASSCUTS=.TRUE.             !EVENT IS OK UNLESS OTHERWISE CHANGED
-      IF (FIRSTTIME) THEN
-         FIRSTTIME=.FALSE.
-c      Preparation for reweighting by setting up clustering by diagrams
-         call initcluster()
-c
-c
-         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'i8)'
-         write(*,formstr) 'Particle',(i,i=nincoming+1,nexternal)
-         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'f8.1)'
-         write(*,formstr) 'Et >',(etmin(i),i=nincoming+1,nexternal)
-         write(*,formstr) 'E >',(emin(i),i=nincoming+1,nexternal)
-         write(*,formstr) 'Eta <',(etamax(i),i=nincoming+1,nexternal)
-         write(*,formstr) 'xqcut: ',(xqcuti(i),i=nincoming+1,nexternal)
-         write(formstr,'(a,i2.2,a)')'(a,i2,a,',nexternal,'f8.1)'
-         do j=nincoming+1,nexternal-1
-            write(*,formstr) 'd R #',j,'  >',(-0.0,i=nincoming+1,j),
-     &           (r2min(i,j),i=j+1,nexternal)
-            do i=j+1,nexternal
-               r2min(i,j)=r2min(i,j)*dabs(r2min(i,j))    !Since r2 returns distance squared
-               r2max(i,j)=r2max(i,j)*dabs(r2max(i,j))
-            enddo
-         enddo
-         do j=nincoming+1,nexternal-1
-            write(*,formstr) 's min #',j,'>',
-     &           (s_min(i,j),i=nincoming+1,nexternal)
-         enddo
-         do j=nincoming+1,nexternal-1
-            write(*,formstr) 'xqcutij #',j,'>',
-     &           (xqcutij(i,j),i=nincoming+1,nexternal)
-         enddo
-
-cc
-cc     Set the strong coupling
-cc
-c         call set_ren_scale(P,scale)
-c
-cc     Check that the user funtions for setting the scales
-cc     have been edited if the choice of an event-by-event
-cc     scale choice has been made 
-c
-c         if(.not.fixed_ren_scale) then
-c            if(scale.eq.0d0) then
-c               write(6,*) 
-c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
-c               write(6,*) ' Dynamical renormalization scale choice '
-c               write(6,*) ' selected but user subroutine' 
-c               write(6,*) ' set_ren_scale not edited in file:setpara.f'
-c               write(6,*) ' Switching to a fixed_ren_scale choice'
-c               write(6,*) ' with scale=zmass'
-c               scale=91.2d0
-c               write(6,*) 'scale=',scale
-c               fixed_ren_scale=.true.
-c               call set_ren_scale(P,scale)
-c            endif
-c         endif
-         
-c         if(.not.fixed_fac_scale) then
-c            call set_fac_scale(P,q2fact)
-c            if(q2fact(1).eq.0d0.or.q2fact(2).eq.0d0) then
-c               write(6,*) 
-c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
-c               write(6,*) ' Dynamical renormalization scale choice '
-c               write(6,*) ' selected but user subroutine' 
-c               write(6,*) ' set_fac_scale not edited in file:setpara.f'
-c               write(6,*) ' Switching to a fixed_fac_scale choice'
-c               write(6,*) ' with q2fact(i)=zmass**2'
-c               fixed_fac_scale=.true.
-c               q2fact(1)=91.2d0**2
-c               q2fact(2)=91.2d0**2
-c               write(6,*) 'scales=',q2fact(1),q2fact(2)
-c            endif
-c         endif
-
-         if(fixed_ren_scale) then
-            G = SQRT(4d0*PI*ALPHAS(scale))
-            call update_as_param()
-         endif
-
-c     Put momenta in the common block to zero to start
-         do i=0,3
-            do j=1,max_particles
-               pp(i,j) = 0d0
-            enddo
-         enddo
-         
-      ENDIF ! IF FIRSTTIME
-
-      IF (CUTSDONE) THEN
-         PASSCUTS=CUTSPASSED
-         RETURN
-      ENDIF
-      CUTSDONE=.TRUE.
-c      CUTSPASSED=.FALSE.
-
-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
-      
-c
-c     Limit S_hat
-c
-c      if (x1*x2*stot .gt. 500**2) then
-c         passcuts=.false.
-c         return
-c      endif
-
-C $B$ DESACTIVATE_CUT $E$ !This is a tag for MadWeight
-
-      if(debug) write (*,*) '============================='
-      if(debug) write (*,*) ' EVENT STARTS TO BE CHECKED  '
-      if(debug) write (*,*) '============================='
-c     
-c     p_t min & max cuts
-c     
-      do i=nincoming+1,nexternal
-         if(debug) write (*,*) 'pt(',i,')=',pt(p(0,i)),'   ',etmin(i),
-     $        ':',etmax(i)
-         notgood=(pt(p(0,i)) .lt. etmin(i)).or.
-     &        (etmax(i).ge.0d0.and.pt(p(0,i)) .gt. etmax(i))
-         if (notgood) then
-            if(debug) write (*,*) i,' -> fails'
-            passcuts=.false.
-            return
-         endif
-      enddo
-c
-c    missing ET min & max cut + Invariant mass of leptons and neutrino 
-c    nb: missing Et defined as the vector sum over the neutrino's pt
-c
-c-- reset ptemp(0:3)
-      do j=0,3
-         ptemp(j)=0 ! for the neutrino
-         ptemp2(j)=0 ! for the leptons
-      enddo
-c-  sum over the momenta
-      do i=nincoming+1,nexternal
-         if(is_a_nu(i)) then            
-         if(debug) write (*,*) i,' -> neutrino '
-            do j=0,3
-               ptemp(j)=ptemp(j)+p(j,i)
-            enddo
-         elseif(is_a_l(i)) then            
-         if(debug) write (*,*) i,' -> lepton '
-            do j=0,3
-               ptemp2(j)=ptemp2(j)+p(j,i)
-            enddo
-         endif
-
-      enddo
-c-  check the et
-      if(debug.and.ptemp(0).eq.0d0) write (*,*) 'No et miss in event'
-      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Et miss =',pt(ptemp(0)),'   ',misset,':',missetmax
-      if(debug.and.ptemp2(0).eq.0d0) write (*,*) 'No leptons in event'
-      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Energy of leptons =',pt(ptemp2(0))
-      if(ptemp(0).gt.0d0) then
-         notgood=(pt(ptemp(0)) .lt. misset).or.
-     &        (missetmax.ge.0d0.and.pt(ptemp(0)) .gt. missetmax)
-         if (notgood) then
-            if(debug) write (*,*) ' missing et cut -> fails'
-            passcuts=.false.
-            return
-         endif
-      endif
-      if (mmnl.gt.0d0.or.mmnlmax.ge.0d0)then
-         if(dsqrt(SumDot(ptemp,ptemp2,1d0)).lt.mmnl.or.mmnlmax.ge.0d0.and.dsqrt(SumDot(ptemp, ptemp2,1d0)).gt.mmnlmax) then
-            if(debug) write (*,*) 'lepton invariant mass -> fails'
-            passcuts=.false.
-            return
-         endif
-      endif
-c
-c     pt cut on heavy particles
-c     gives min(pt) for (at least) one heavy particle
-c
-      if(ptheavy.gt.0d0)then
-         passcuts=.false.
-         foundheavy=.false.
-         do i=nincoming+1,nexternal
-            if(is_heavy(i)) then            
-               if(debug) write (*,*) i,' -> heavy '
-               foundheavy=.true.
-               if(pt(p(0,i)).gt.ptheavy) passcuts=.true.
-            endif
-         enddo
-         
-         if(.not.passcuts.and.foundheavy)then
-            if(debug) write (*,*) ' heavy particle cut -> fails'
-            return
-         else
-            passcuts=.true.
-         endif
-      endif
-c     
-c     E min & max cuts
-c     
-      do i=nincoming+1,nexternal
-         if(debug) write (*,*) 'p(0,',i,')=',p(0,i),'   ',emin(i),':',emax(i)
-         notgood=(p(0,i) .le. emin(i)).or.
-     &        (emax(i).ge.0d0 .and. p(0,i) .gt. emax(i))
-         if (notgood) then
-            if(debug) write (*,*) i,' -> fails'
-            passcuts=.false.
-            return
-         endif
-      enddo
-c     
-c     Rapidity  min & max cuts
-c     
-      do i=nincoming+1,nexternal
-         if(debug) write (*,*) 'abs(rap(',i,'))=',abs(rap(p(0,i))),'   ',etamin(i),':',etamax(i)
-         notgood=(etamax(i).ge.0.and.abs(rap(p(0,i))) .gt. etamax(i)).or.
-     &        (abs(rap(p(0,i))) .lt. etamin(i))
-         if (notgood) then
-            if(debug) write (*,*) i,' -> fails'
-            passcuts=.false.
-            return
-         endif
-      enddo
-c     
-c     DeltaR min & max cuts
-c     
-      do i=nincoming+1,nexternal-1
-         do j=i+1,nexternal
-            if(debug) write (*,*) 'r2(',i, ',' ,j,')=',dsqrt(r2(p(0,i),p(0,j)))
-            if(debug) write (*,*) dsqrt(r2min(j,i)),dsqrt(r2max(j,i))
-            if(r2min(j,i).gt.0.or.r2max(j,i).ge.0d0) then
-               tmp=r2(p(0,i),p(0,j))
-               notgood=(tmp .lt. r2min(j,i)).or.
-     $              (r2max(j,i).ge.0d0 .and. tmp .gt. r2max(j,i))
-               if (notgood) then
-                  if(debug) write (*,*) i,j,' -> fails'
-                  passcuts=.false.
-                  return
-               endif
-            endif
-         enddo
-      enddo
-
-
-c     s-channel min & max pt of sum of 4-momenta
-c     
-      do i=nincoming+1,nexternal-1
-         do j=i+1,nexternal
-            if(debug)write (*,*) 'ptll(',i,',',j,')=',PtDot(p(0,i),p(0,j))
-            if(debug)write (*,*) ptll_min(j,i),ptll_max(j,i)
-            if(ptll_min(j,i).gt.0.or.ptll_max(j,i).ge.0d0) then
-               tmp=PtDot(p(0,i),p(0,j))
-               notgood=(tmp .lt. ptll_min(j,i).or.
-     $              ptll_max(j,i).ge.0d0 .and. tmp.gt.ptll_max(j,i))
-               if (notgood) then
-                  if(debug) write (*,*) i,j,' -> fails'
-                  passcuts=.false.
-                  return
-               endif
-            endif
-         enddo
-      enddo
-
-
-
-
-c     
-c     s-channel min & max invariant mass cuts
-c     
-      do i=nincoming+1,nexternal-1
-         do j=i+1,nexternal
-            if(debug) write (*,*) 's(',i,',',j,')=',Sumdot(p(0,i),p(0,j),+1d0)
-            if(debug) write (*,*) s_min(j,i),s_max(j,i)
-            if(s_min(j,i).gt.0.or.s_max(j,i).ge.0d0) then
-               tmp=SumDot(p(0,i),p(0,j),+1d0)
-               if(s_min(j,i).le.s_max(j,i) .or. s_max(j,i).lt.0d0)then
-                  notgood=(tmp .lt. s_min(j,i).or.
-     $                 s_max(j,i).ge.0d0 .and. tmp .gt. s_max(j,i)) 
-                  if (notgood) then
-                     if(debug) write (*,*) i,j,' -> fails'
-                     passcuts=.false.
-                     return
-                  endif
-               else
-                  notgood=(tmp .lt. s_min(j,i).and.tmp .gt. s_max(j,i)) 
-                  if (notgood) then
-                     if(debug) write (*,*) i,j,' -> fails'
-                     passcuts=.false.
-                     return
-                  endif
-               endif
-            endif
-         enddo
-      enddo
-C     $B$DESACTIVATE_BW_CUT$B$ This is a Tag for MadWeight
-c     
-c     B.W. phase space cuts
-c     
-      pass_bw=cut_bw(p)
-c     JA 4/8/11 always check pass_bw
-      if ( pass_bw.and..not.CUTSPASSED) then
-         passcuts=.false.
-         if(debug) write (*,*) ' pass_bw -> fails'
-         return
-      endif
-C     $E$DESACTIVATE_BW_CUT$E$ This is a Tag for MadWeight
-        CUTSPASSED=.FALSE.
-C     
-C     maximal and minimal pt of the jets sorted by pt
-c     
-      njets=0
-      nheavyjets=0
-
-c- fill ptjet with the pt's of the jets.
-      do i=nincoming+1,nexternal
-         if(is_a_j(i)) then
-            njets=njets+1
-            ptjet(njets)=pt(p(0,i))
-         endif
-         if(is_a_b(i)) then
-            nheavyjets=nheavyjets+1
-            ptheavyjet(nheavyjets)=pt(p(0,i))
-         endif
-
-      enddo
-      if(debug) write (*,*) 'not yet ordered ',njets,'   ',ptjet
-
-C----------------------------------------------------------------------------
-C     DURHAM_KT CUT
-C----------------------------------------------------------------------------
-      IF(NJETS.GT.0 .AND.KT_DURHAM.GT.0D0) THEN
-C RESET JET MOMENTA
-      njets=0
-      DO I=1,NEXTERNAL
-        DO J=0,3
-          PJET(I,J) = 0E0
-        ENDDO
-      ENDDO
-
-      do i=nincoming+1,nexternal
-         if(is_a_j(i).and.do_cuts(i)) then
-           njets=njets+1
-           DO J=0,3
-             PJET(NJETS,J) = P(J,I)
-           ENDDO
-         endif
-      enddo
-
-C DURHAM KT SEPARATION CUT
-
-
-        PTMIN = EBEAM(1) + EBEAM(2)
-
-        DO I=1,NJETS
-
-C         PT WITH RESPECT TO Z AXIS FOR HADRONIC COLLISIONS
-          IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
-            PT1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2)
-            PTMIN = MIN( PTMIN, PT1 )
-          ENDIF
-
-          DO J=I+1,NJETS
-C           GET ANGLE BETWEEN JETS
-            PABS1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2 + PJET(I,3)**2)
-            PABS2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2 + PJET(J,3)**2)
-C           CHECK IF 3-MOMENTA DO NOT VANISH
-            IF(PABS1*PABS2 .NE. 0D0) THEN
-              COSTH = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) + PJET(I,3)*PJET(J,3) )/(PABS1*PABS2)
-            ELSE
-C           IF 3-MOMENTA VANISH, MAKE JET COSTH = 1D0 SO THAT JET MEASURE VANISHES
-              COSTH = 1D0
-            ENDIF
-C           GET PT AND ETA OF JETS
-            PT2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2)
-            ETA1 = 0.5D0*LOG( (PJET(I,0) + PJET(I,3)) / (PJET(I,0) - PJET(I,3)) )
-            ETA2 = 0.5D0*LOG( (PJET(J,0) + PJET(J,3)) / (PJET(J,0) - PJET(J,3)) )
-C           GET COSH OF DELTA ETA, COS OF DELTA PHI
-            COSH_DETA = DCOSH( ETA1 - ETA2 )
-            COS_DPHI = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) ) / (PT1*PT2)
-            DPHI = DACOS( COS_DPHI )
-            IF ( (LPP(1).EQ.0) .AND. (LPP(2).EQ.0)) THEN
-C             KT FOR E+E- COLLISION
-              PTNOW = DSQRT( 2D0*MIN(PJET(I,0)**2,PJET(J,0)**2)*( 1D0-COSTH ) )
-             ELSE
-C             HADRONIC KT, FASTJET DEFINITION
-              PTNOW = DSQRT( MIN(PT1**2,PT2**2)*( (ETA1 - ETA2 )**2 + DPHI**2 )/(D_PARAMETER**2) )
-            ENDIF
-
-            PTMIN = MIN( PTMIN, PTNOW )
-
-          ENDDO ! LOOP OVER NJET
-
-        ENDDO ! LOOP OVER NJET
-
-C CHECK COMPATIBILITY WITH CUT
-        IF( (PTMIN .LT. KT_DURHAM)) THEN
-          PASSCUTS = .FALSE.
-          RETURN
-        ENDIF
-      ENDIF ! IF NJETS.GT. 0 .AND. DO_KT_DURHAM
-
-C----------------------------------------------------------------------------
-C----------------------------------------------------------------------------
-
-
-
-c- check existance of jets if jet cuts are on
-      if(njets.lt.1.and.(htjmin.gt.0.or.ptj1min.gt.0).or.
-     $     njets.lt.2.and.ptj2min.gt.0.or.
-     $     njets.lt.3.and.ptj3min.gt.0.or.
-     $     njets.lt.4.and.ptj4min.gt.0)then
-         if(debug) write (*,*) i, ' too few jets -> fails'
-         passcuts=.false.
-         return
-      endif
-
-c - sort jet pts
-      do i=1,njets-1
-         do j=i+1,njets
-            if(ptjet(j).gt.ptjet(i)) then
-               temp=ptjet(i)
-               ptjet(i)=ptjet(j)
-               ptjet(j)=temp
-            endif
-         enddo
-      enddo
-      if(debug) write (*,*) 'ordered ',njets,'   ',ptjet
-c
-c     Use "and" or "or" prescriptions 
-c     
-      inclht=0
-
-      if(njets.gt.0) then
-
-       notgood=.not.jetor
-       if(debug) write (*,*) 'jetor :',jetor  
-       if(debug) write (*,*) '0',notgood   
-      
-      do i=1,min(njets,4)
-            if(debug) write (*,*) i,ptjet(i), '   ',ptjmin4(i),
-     $        ':',ptjmax4(i)
-         if(jetor) then     
-c---  if one of the jets does not pass, the event is rejected
-            notgood=notgood.or.(ptjmax4(i).ge.0d0 .and.
-     $           ptjet(i).gt.ptjmax4(i)) .or.
-     $           (ptjet(i).lt.ptjmin4(i))
-            if(debug) write (*,*) i,' notgood total:', notgood   
-         else
-c---  all cuts must fail to reject the event
-            notgood=notgood.and.(ptjmax4(i).ge.0d0 .and.
-     $              ptjet(i).gt.ptjmax4(i) .or.
-     $              (ptjet(i).lt.ptjmin4(i)))
-            if(debug) write (*,*) i,' notgood total:', notgood   
-         endif
-      enddo
-
-
-      if (notgood) then
-         if(debug) write (*,*) i, ' multiple pt -> fails'
-         passcuts=.false.
-         return
-      endif
-
-c---------------------------
-c      Ht cuts
-C---------------------------
-      htj=ptjet(1)
-
-      do i=2,njets
-         htj=htj+ptjet(i)
-         if(debug) write (*,*) i, 'htj ',htj
-         if(debug.and.i.le.4) write (*,*) 'htmin ',i,' ', htjmin4(i),':',htjmax4(i)
-         if(i.le.4)then
-            if(htj.lt.htjmin4(i) .or.
-     $        htjmax4(i).ge.0d0.and.htj.gt.htjmax4(i)) then
-            if(debug) write (*,*) i, ' ht -> fails'
-            passcuts=.false.
-            return
-            endif
-         endif
-      enddo
-
-      if(htj.lt.htjmin.or.htjmax.ge.0d0.and.htj.gt.htjmax)then
-         if(debug) write (*,*) i, ' htj -> fails'
-         passcuts=.false.
-         return
-      endif
-
-      inclht=htj
-
-      endif !if there are jets 
-
-      if(nheavyjets.gt.0) then
-         do i=1,nheavyjets
-            inclht=inclht+ptheavyjet(i)
-         enddo
-      endif !if there are heavyjets
-
-      if(inclht.lt.inclHtmin.or.
-     $     inclHtmax.ge.0d0.and.inclht.gt.inclHtmax)then
-         if(debug) write (*,*) ' inclhtmin=',inclHtmin,' -> fails'
-         passcuts=.false.
-         return
-      endif
-
-C     
-C     maximal and minimal pt of the leptons sorted by pt
-c     
-      nleptons=0
-
-      if(ptl1min.gt.0.or.ptl2min.gt.0.or.ptl3min.gt.0.or.ptl4min.gt.0.or.
-     $     ptl1max.ge.0d0.or.ptl2max.ge.0d0.or.
-     $     ptl3max.ge.0d0.or.ptl4max.ge.0d0) then
-
-c     - fill ptlepton with the pt's of the leptons.
-         do i=nincoming+1,nexternal
-            if(is_a_l(i)) then
-               nleptons=nleptons+1
-               ptlepton(nleptons)=pt(p(0,i))
-            endif
-         enddo
-         if(debug) write (*,*) 'not yet ordered ',njets,'   ',ptjet
-
-c     - check existance of leptons if lepton cuts are on
-         if(nleptons.lt.1.and.ptl1min.gt.0.or.
-     $        nleptons.lt.2.and.ptl2min.gt.0.or.
-     $        nleptons.lt.3.and.ptl3min.gt.0.or.
-     $        nleptons.lt.4.and.ptl4min.gt.0)then
-            if(debug) write (*,*) i, ' too few leptons -> fails'
-            passcuts=.false.
-            return
-         endif
-
-c     - sort lepton pts
-         do i=1,nleptons-1
-            do j=i+1,nleptons
-               if(ptlepton(j).gt.ptlepton(i)) then
-                  temp=ptlepton(i)
-                  ptlepton(i)=ptlepton(j)
-                  ptlepton(j)=temp
-               endif
-            enddo
-         enddo
-         if(debug) write (*,*) 'ordered ',nleptons,'   ',ptlepton
-
-         if(nleptons.gt.0) then
-
-            notgood = .false.
-            do i=1,min(nleptons,4)
-               if(debug) write (*,*) i,ptlepton(i), '   ',ptlmin4(i),':',ptlmax4(i)
-c---  if one of the leptons does not pass, the event is rejected
-               notgood=notgood.or.
-     $              (ptlmax4(i).ge.0d0.and.ptlepton(i).gt.ptlmax4(i)).or.
-     $              (ptlepton(i).lt.ptlmin4(i))
-               if(debug) write (*,*) i,' notgood total:', notgood   
-            enddo
-
-
-            if (notgood) then
-               if(debug) write (*,*) i, ' multiple pt -> fails'
-               passcuts=.false.
-               return
-            endif
-         endif
-      endif
-C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>     
-C     SPECIAL CUTS
-C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-C     REQUIRE AT LEAST ONE JET WITH PT>XPTJ
-         
-         IF(xptj.gt.0d0) THEN
-            xvar=0
-            do i=nincoming+1,nexternal
-               if(is_a_j(i)) xvar=max(xvar,pt(p(0,i)))
-            enddo
-            if (xvar .lt. xptj) then
-               passcuts=.false.
-               return
-            endif
-         ENDIF
-
-C     REQUIRE AT LEAST ONE PHOTON WITH PT>XPTA
-         
-         IF(xpta.gt.0d0) THEN
-            xvar=0
-            do i=nincoming+1,nexternal
-               if(is_a_a(i)) xvar=max(xvar,pt(p(0,i)))
-            enddo
-            if (xvar .lt. xpta) then
-               passcuts=.false.
-               return
-            endif
-         ENDIF
-
-C     REQUIRE AT LEAST ONE B  WITH PT>XPTB
-         
-         IF(xptb.gt.0d0) THEN
-            xvar=0
-            do i=nincoming+1,nexternal
-               if(is_a_b(i)) xvar=max(xvar,pt(p(0,i)))
-            enddo
-            if (xvar .lt. xptb) then
-               passcuts=.false.
-               return
-            endif
-         ENDIF
-
-C     REQUIRE AT LEAST ONE LEPTON  WITH PT>XPTL
-         
-         IF(xptl.gt.0d0) THEN
-            xvar=0
-            do i=nincoming+1,nexternal
-               if(is_a_l(i)) xvar=max(xvar,pt(p(0,i)))
-            enddo
-            if (xvar .lt. xptl) then
-               passcuts=.false.
-               if(debug) write (*,*) ' xptl -> fails'
-               return
-            endif
-         ENDIF
-C
-C     WBF CUTS: TWO TYPES
-C    
-C     FIRST TYPE:  implemented by FM
-C
-C     1. FIND THE 2 HARDEST JETS
-C     2. REQUIRE |RAP(J)|>XETAMIN
-C     3. REQUIRE RAP(J1)*ETA(J2)<0
-C
-C     SECOND TYPE : added by Simon de Visscher 1-08-2007
-C
-C     1. FIND THE 2 HARDEST JETS
-C     2. REQUIRE |RAP(J1)-RAP(J2)|>DELTAETA
-C     3. REQUIRE RAP(J1)*RAP(J2)<0
-C
-C
-         hardj1=0
-         hardj2=0
-         ptmax1=0d0
-         ptmax2=0d0
-
-C-- START IF AT LEAST ONE OF THE CUTS IS ACTIVATED
-         
-         IF(XETAMIN.GT.0D0.OR.DELTAETA.GT.0D0) THEN
-            
-C-- FIND THE HARDEST JETS
-
-            do i=nincoming+1,nexternal
-               if(is_a_j(i)) then
-c                  write (*,*) i,pt(p(0,i))
-                  if(pt(p(0,i)).gt.ptmax1) then
-                     hardj2=hardj1
-                     ptmax2=ptmax1
-                     hardj1=i
-                     ptmax1=pt(p(0,i))
-                  elseif(pt(p(0,i)).gt.ptmax2) then
-                     hardj2=i
-                     ptmax2=pt(p(0,i))
-                  endif
-c                  write (*,*) hardj1,hardj2,ptmax1,ptmax2
-               endif
-            enddo
-            
-            if (hardj2.eq.0) goto 21 ! bypass vbf cut since not enough jets
-
-C-- NOW APPLY THE CUT I            
-
-            if (abs(rap(p(0,hardj1))) .lt. xetamin
-     &       .or.abs(rap(p(0,hardj2))) .lt. xetamin
-     &       .or.rap(p(0,hardj1))*rap(p(0,hardj2)) .gt.0d0) then
-             passcuts=.false.
-             return
-            endif
-
-            
-C-- NOW APPLY THE CUT II
-            
-            if (abs(rap(p(0,hardj1))-rap(p(0,hardj2))) .lt. deltaeta) then
-             passcuts=.false.
-             return
-            endif
-         
-c            write (*,*) hardj1,hardj2,rap(p(0,hardj1)),rap(p(0,hardj2))
-         
-         ENDIF
-
-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 If we do not require a mimimum jet energy, there's no need to apply
-c jet clustering and all that.
-      if (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
-c the number of (light) QCD partons at the real-emission level (i.e. one
-c more than the Born).
-
-      nQCD=0
-      do j=nincoming+1,nexternal
-         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
-            enddo
-         endif
-      enddo
-
-        nph=0
-        do j=nincoming+1,nexternal
-          if (is_a_a(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_l(j)) then
-              nem=nem+1
-              do i=0,3
-                pem(i,nem)=p(i,j) ! Use C.o.M. frame momenta
-              enddo
-            endif
-          enddo
-        endif
-
-        alliso=.true.
-
-        j=0
-        dowhile(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
-
-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)))
-            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
-      endif
-
- 444    continue
-c End photon isolation
-
-
-C...Set couplings if event passed cuts
-
- 21   if(.not.fixed_ren_scale) then
-         call set_ren_scale(P,scale)
-         if(scale.gt.0) G = SQRT(4d0*PI*ALPHAS(scale))
-      endif
-
-      if(.not.fixed_fac_scale) then
-         call set_fac_scale(P,q2fact)
-      endif
-
-c
-c     Here we cluster event and reset factorization and renormalization
-c     scales on an event-by-event basis, as well as check xqcut for jets
-c
-c     Note the following condition is the first line of setclscales
-c      if(xqcut.gt.0d0.or.ickkw.gt.0.or.scale.eq.0.or.q2fact(1).eq.0)then
-c     Do not duplicate it since some variable are set for syscalc in the fct
-        if(.not.setclscales(p))then
-         if(debug) write (*,*) ' setclscales -> fails'
-         passcuts=.false.
-         return
-       endif
-c      endif
-
-c     Set couplings in model files
-      if(.not.fixed_ren_scale.or..not.fixed_couplings) then
-         if (.not.fixed_couplings)then
-            do i=0,3
-               do j=1,nexternal
-                  pp(i,j)=p(i,j)
-               enddo
-            enddo
-         endif
-         call update_as_param()
-      endif
-
-      IF (FIRSTTIME2) THEN
-        FIRSTTIME2=.FALSE.
-        write(6,*) 'alpha_s for scale ',scale,' is ', G**2/(16d0*atan(1d0))
-      ENDIF
-
-      if(debug) write (*,*) '============================='
-      if(debug) write (*,*) ' EVENT PASSED THE CUTS       '
-      if(debug) write (*,*) '============================='
-
-      CUTSPASSED=.TRUE.
-
-      RETURN
-      END
-
-
-C
-C     FUNCTION FOR ISOLATION
-C
-
-      function iso_getdrv40(p1,p2)
-      implicit none
-      real*8 iso_getdrv40,p1(0:3),p2(0:3)
-      real*8 iso_getdr
-c
-      iso_getdrv40=iso_getdr(p1(0),p1(1),p1(2),p1(3),
-     #                       p2(0),p2(1),p2(2),p2(3))
-      return
-      end
-
-
-      function iso_getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
-      implicit none
-      real*8 iso_getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
-     # iso_getpseudorap,iso_getdelphi
-c
-      deta=iso_getpseudorap(en1,ptx1,pty1,pl1)-
-     #     iso_getpseudorap(en2,ptx2,pty2,pl2)
-      dphi=iso_getdelphi(ptx1,pty1,ptx2,pty2)
-      iso_getdr=sqrt(dphi**2+deta**2)
-      return
-      end
-
-
-      function iso_getpseudorap(en,ptx,pty,pl)
-      implicit none
-      real*8 iso_getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
-      parameter (tiny=1.d-5)
-c
-      pt=sqrt(ptx**2+pty**2)
-      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
-        eta=sign(1.d0,pl)*1.d8
-      else
-        th=atan2(pt,pl)
-        eta=-log(tan(th/2.d0))
-      endif
-      iso_getpseudorap=eta
-      return
-      end
-
-
-      function iso_getdelphi(ptx1,pty1,ptx2,pty2)
-      implicit none
-      real*8 iso_getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
-      parameter (tiny=1.d-5)
-c
-      pt1=sqrt(ptx1**2+pty1**2)
-      pt2=sqrt(ptx2**2+pty2**2)
-      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
-        tmp=ptx1*ptx2+pty1*pty2
-        tmp=tmp/(pt1*pt2)
-        if(abs(tmp).gt.1.d0+tiny)then
-          write(*,*)'Cosine larger than 1'
-          stop
-        elseif(abs(tmp).ge.1.d0)then
-          tmp=sign(1.d0,tmp)
-        endif
-        tmp=acos(tmp)
-      else
-        tmp=1.d8
-      endif
-      iso_getdelphi=tmp
-      return
-      end
-
-      function chi_gamma_iso(dr,R0,xn,epsgamma,pTgamma)
-c Eq.(3.4) of Phys.Lett. B429 (1998) 369-374 [hep-ph/9801442]
-      implicit none
-      real*8 chi_gamma_iso,dr,R0,xn,epsgamma,pTgamma
-      real*8 tmp,axn
-c
-      axn=abs(xn)
-      tmp=epsgamma*pTgamma
-      if(axn.ne.0.d0)then
-        tmp=tmp*( (1-cos(dr))/(1-cos(R0)) )**axn
-      endif
-      chi_gamma_iso=tmp
-      return
-      end
-
-
-*
-* $Id: sortzv.F,v 1.1.1.1 1996/02/15 17:49:50 mclareni Exp $
-*
-* $Log: sortzv.F,v $
-* Revision 1.1.1.1  1996/02/15 17:49:50  mclareni
-* Kernlib
-*
-*
-c$$$#include "kerngen/pilot.h"
-      SUBROUTINE SORTZV (A,INDEX,N1,MODE,NWAY,NSORT)
-C
-C CERN PROGLIB# M101    SORTZV          .VERSION KERNFOR  3.15  820113
-C ORIG. 02/10/75
-C
-      DIMENSION A(N1),INDEX(N1)
-C
-C
-      N = N1
-      IF (N.LE.0)            RETURN
-      IF (NSORT.NE.0) GO TO 2
-      DO 1 I=1,N
-    1 INDEX(I)=I
-C
-    2 IF (N.EQ.1)            RETURN
-      IF (MODE)    10,20,30
-   10 CALL SORTTI (A,INDEX,N)
-      GO TO 40
-C
-   20 CALL SORTTC(A,INDEX,N)
-      GO TO 40
-C
-   30 CALL SORTTF (A,INDEX,N)
-C
-   40 IF (NWAY.EQ.0) GO TO 50
-      N2 = N/2
-      DO 41 I=1,N2
-      ISWAP = INDEX(I)
-      K = N+1-I
-      INDEX(I) = INDEX(K)
-   41 INDEX(K) = ISWAP
-   50 RETURN
-      END
-*     ========================================
-      SUBROUTINE SORTTF (A,INDEX,N1)
-C
-      DIMENSION A(N1),INDEX(N1)
-C
-      N = N1
-      DO 3 I1=2,N
-      I3 = I1
-      I33 = INDEX(I3)
-      AI = A(I33)
-    1 I2 = I3/2
-      IF (I2) 3,3,2
-    2 I22 = INDEX(I2)
-      IF (AI.LE.A (I22)) GO TO 3
-      INDEX (I3) = I22
-      I3 = I2
-      GO TO 1
-    3 INDEX (I3) = I33
-    4 I3 = INDEX (N)
-      INDEX (N) = INDEX (1)
-      AI = A(I3)
-      N = N-1
-      IF (N-1) 12,12,5
-    5 I1 = 1
-    6 I2 = I1 + I1
-      IF (I2.LE.N) I22= INDEX(I2)
-      IF (I2-N) 7,9,11
-    7 I222 = INDEX (I2+1)
-      IF (A(I22)-A(I222)) 8,9,9
-    8 I2 = I2+1
-      I22 = I222
-    9 IF (AI-A(I22)) 10,11,11
-   10 INDEX(I1) = I22
-      I1 = I2
-      GO TO 6
-   11 INDEX (I1) = I3
-      GO TO 4
-   12 INDEX (1) = I3
-      RETURN
-      END
-*     ========================================
-      SUBROUTINE SORTTI (A,INDEX,N1)
-C
-      INTEGER A,AI
-      DIMENSION A(N1),INDEX(N1)
-C
-      N = N1
-      DO 3 I1=2,N
-      I3 = I1
-      I33 = INDEX(I3)
-      AI = A(I33)
-    1 I2 = I3/2
-      IF (I2) 3,3,2
-    2 I22 = INDEX(I2)
-      IF (AI.LE.A (I22)) GO TO 3
-      INDEX (I3) = I22
-      I3 = I2
-      GO TO 1
-    3 INDEX (I3) = I33
-    4 I3 = INDEX (N)
-      INDEX (N) = INDEX (1)
-      AI = A(I3)
-      N = N-1
-      IF (N-1) 12,12,5
-    5 I1 = 1
-    6 I2 = I1 + I1
-      IF (I2.LE.N) I22= INDEX(I2)
-      IF (I2-N) 7,9,11
-    7 I222 = INDEX (I2+1)
-      IF (A(I22)-A(I222)) 8,9,9
-    8 I2 = I2+1
-      I22 = I222
-    9 IF (AI-A(I22)) 10,11,11
-   10 INDEX(I1) = I22
-      I1 = I2
-      GO TO 6
-   11 INDEX (I1) = I3
-      GO TO 4
-   12 INDEX (1) = I3
-      RETURN
-      END
-*     ========================================
-      SUBROUTINE SORTTC (A,INDEX,N1)
-C
-      INTEGER A,AI
-      DIMENSION A(N1),INDEX(N1)
-C
-      N = N1
-      DO 3 I1=2,N
-      I3 = I1
-      I33 = INDEX(I3)
-      AI = A(I33)
-    1 I2 = I3/2
-      IF (I2) 3,3,2
-    2 I22 = INDEX(I2)
-      IF(ICMPCH(AI,A(I22)))3,3,21
-   21 INDEX (I3) = I22
-      I3 = I2
-      GO TO 1
-    3 INDEX (I3) = I33
-    4 I3 = INDEX (N)
-      INDEX (N) = INDEX (1)
-      AI = A(I3)
-      N = N-1
-      IF (N-1) 12,12,5
-    5 I1 = 1
-    6 I2 = I1 + I1
-      IF (I2.LE.N) I22= INDEX(I2)
-      IF (I2-N) 7,9,11
-    7 I222 = INDEX (I2+1)
-      IF (ICMPCH(A(I22),A(I222))) 8,9,9
-    8 I2 = I2+1
-      I22 = I222
-    9 IF (ICMPCH(AI,A(I22))) 10,11,11
-   10 INDEX(I1) = I22
-      I1 = I2
-      GO TO 6
-   11 INDEX (I1) = I3
-      GO TO 4
-   12 INDEX (1) = I3
-      RETURN
-      END
-*     ========================================
-      FUNCTION ICMPCH(IC1,IC2)
-C     FUNCTION TO COMPARE TWO 4 CHARACTER EBCDIC STRINGS - IC1,IC2
-C     ICMPCH=-1 IF HEX VALUE OF IC1 IS LESS THAN IC2
-C     ICMPCH=0  IF HEX VALUES OF IC1 AND IC2 ARE THE SAME
-C     ICMPCH=+1 IF HEX VALUES OF IC1 IS GREATER THAN IC2
-      I1=IC1
-      I2=IC2
-      IF(I1.GE.0.AND.I2.GE.0)GOTO 40
-      IF(I1.GE.0)GOTO 60
-      IF(I2.GE.0)GOTO 80
-      I1=-I1
-      I2=-I2
-      IF(I1-I2)80,70,60
- 40   IF(I1-I2)60,70,80
- 60   ICMPCH=-1
-      RETURN
- 70   ICMPCH=0
-      RETURN
- 80   ICMPCH=1
-      RETURN
-      END
-
+      logical function pass_point(p)
+c************************************************************************
+c     This function is called from sample to see if it needs to 
+c     bother calculating the weight from all the different conficurations
+c     You can either just return true, or have it call passcuts
+c************************************************************************
+      implicit none
+c
+c     Arguments
+c
+      double precision p
+c
+c     External
+c
+      logical passcuts
+      external passcuts
+c-----
+c  Begin Code
+c-----
+      pass_point = .true.
+c      pass_point = passcuts(p)
+      end
+C 
+      LOGICAL FUNCTION PASSCUTS(P)
+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 ...
+C            ALL MOMENTA ARE IN THE REST FRAME!!
+C            COMMON/JETCUTS/   CUTS ON JETS
+C     OUTPUT:
+C            TRUE IF EVENTS PASSES ALL CUTS LISTED
+C**************************************************************************
+      IMPLICIT NONE
+c
+c     Constants
+c
+      include 'genps.inc'
+      include 'nexternal.inc'
+C
+C     ARGUMENTS
+C
+      REAL*8 P(0:3,nexternal)
+
+C
+C     LOCAL
+C
+      LOGICAL FIRSTTIME,FIRSTTIME2,pass_bw,notgood,good,foundheavy
+      LOGICAL DEBUG
+      integer i,j,njets,nheavyjets,nleptons,hardj1,hardj2
+      REAL*8 XVAR,ptmax1,ptmax2,htj,tmp,inclht
+      real*8 ptemp(0:3), ptemp2(0:3)
+      character*20 formstr
+C
+C     PARAMETERS
+C
+      real*8 PI
+      parameter( PI = 3.14159265358979323846d0 )
+C
+C     EXTERNAL
+C
+      REAL*8 R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,PtDot
+      logical cut_bw,setclscales
+      external R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,cut_bw,setclscales,PtDot
+C
+C     GLOBAL
+C
+      include 'run.inc'
+      include 'cuts.inc'
+
+      double precision ptjet(nexternal)
+      double precision ptheavyjet(nexternal)
+      double precision ptlepton(nexternal)
+      double precision temp
+
+C     VARIABLES TO SPECIFY JETS
+      DOUBLE PRECISION PJET(NEXTERNAL,0:3)
+      DOUBLE PRECISION PTMIN
+      DOUBLE PRECISION PT1,PT2
+
+C     INTEGERS FOR COUNTING.
+      INTEGER K,L,J1,J2
+
+C VARIABLES FOR KT CUT
+      DOUBLE PRECISION PTNOW,COSTH,PABS1,PABS2
+      DOUBLE PRECISION ETA1,ETA2,COSH_DETA,COS_DPHI,KT1SQ,KT2SQ, DPHI
+
+      double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
+      double precision emin(nincoming+1:nexternal)
+      double precision                    r2min(nincoming+1:nexternal,nincoming+1:nexternal)
+      double precision s_min(nexternal,nexternal)
+      double precision etmax(nincoming+1:nexternal),etamin(nincoming+1:nexternal)
+      double precision emax(nincoming+1:nexternal)
+      double precision r2max(nincoming+1:nexternal,nincoming+1:nexternal)
+      double precision s_max(nexternal,nexternal)
+      double precision ptll_min(nexternal,nexternal),ptll_max(nexternal,nexternal)
+      double precision inclHtmin,inclHtmax
+      common/to_cuts/  etmin, emin, etamax, r2min, s_min,
+     $     etmax, emax, etamin, r2max, s_max, ptll_min, ptll_max, inclHtmin,inclHtmax
+
+      double precision ptjmin4(4),ptjmax4(4),htjmin4(2:4),htjmax4(2:4)
+      logical jetor
+      common/to_jet_cuts/ ptjmin4,ptjmax4,htjmin4,htjmax4,jetor
+
+      double precision ptlmin4(4),ptlmax4(4)
+      common/to_lepton_cuts/ ptlmin4,ptlmax4
+
+c
+c     Special cuts
+c
+
+      integer        lbw(0:nexternal)  !Use of B.W.
+      common /to_BW/ lbw
+C
+C     SPECIAL CUTS
+C
+      LOGICAL  IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
+      LOGICAL  IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
+      LOGICAL  IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
+      logical  do_cuts(nexternal)
+      COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
+     . IS_A_ONIUM, do_cuts
+
+C
+C     MERGING SCALE CUT
+C
+C     Retrieve which external particles undergo the ktdurham and ptlund cuts.
+      LOGICAL  is_pdg_for_merging_cut(NEXTERNAL)
+      logical from_decay(-(nexternal+3):nexternal)
+      COMMON /TO_MERGE_CUTS/is_pdg_for_merging_cut, from_decay
+
+C
+C     ADDITIONAL VARIABLES FOR PTLUND CUT
+C
+      INTEGER NMASSLESS
+      DOUBLE PRECISION PINC(NINCOMING,0:3)
+      DOUBLE PRECISION PRADTEMP(0:3), PRECTEMP(0:3), PEMTTEMP(0:3)
+      DOUBLE PRECISION PTMINSAVE, RHOPYTHIA
+      EXTERNAL RHOPYTHIA
+C
+C     FLAVOUR INFORMATION NECESSARY TO RECONSTRUCT PTLUND
+C
+      INTEGER JETFLAVOUR(NEXTERNAL), INCFLAVOUR(NINCOMING)
+      include 'maxamps.inc'
+      integer idup(nexternal,maxproc,maxsproc)
+      integer mothup(2,nexternal)
+      integer icolup(2,nexternal,maxflow,maxsproc)
+      include 'leshouche.inc'
+
+C
+C     Keep track of whether cuts already calculated for this event
+C
+      LOGICAL CUTSDONE,CUTSPASSED
+      COMMON/TO_CUTSDONE/CUTSDONE,CUTSPASSED
+      DATA CUTSDONE,CUTSPASSED/.FALSE.,.FALSE./
+
+C $B$ MW_NEW_DEF $E$ !this is a tag for MadWeight
+
+      double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
+      common/to_xqcuts/xqcutij,xqcuti
+
+c jet cluster algorithm
+      integer nQCD !,NJET,JET(nexternal)
+c      double precision plab(0:3, nexternal)
+      double precision pQCD(0:3,nexternal)!,PJET(0:3,nexternal)
+c      double precision rfj,sycut,palg,fastjetdmerge
+c      integer njet_eta
+c     Photon isolation
+      integer nph,nem,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)
+      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)
+
+      include 'coupl.inc'
+
+C
+C
+c
+      DATA FIRSTTIME,FIRSTTIME2/.TRUE.,.TRUE./
+
+c put momenta in common block for couplings.f
+      double precision pp(0:3,max_particles)
+      common /momenta_pp/pp
+
+      DATA DEBUG/.FALSE./
+
+C-----
+C  BEGIN CODE
+C-----
+
+
+
+      PASSCUTS=.TRUE.             !EVENT IS OK UNLESS OTHERWISE CHANGED
+      IF (FIRSTTIME) THEN
+         FIRSTTIME=.FALSE.
+c      Preparation for reweighting by setting up clustering by diagrams
+         call initcluster()
+c
+c
+         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'i8)'
+         write(*,formstr) 'Particle',(i,i=nincoming+1,nexternal)
+         write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'f8.1)'
+         write(*,formstr) 'Et >',(etmin(i),i=nincoming+1,nexternal)
+         write(*,formstr) 'E >',(emin(i),i=nincoming+1,nexternal)
+         write(*,formstr) 'Eta <',(etamax(i),i=nincoming+1,nexternal)
+         write(*,formstr) 'xqcut: ',(xqcuti(i),i=nincoming+1,nexternal)
+         write(formstr,'(a,i2.2,a)')'(a,i2,a,',nexternal,'f8.1)'
+         do j=nincoming+1,nexternal-1
+            write(*,formstr) 'd R #',j,'  >',(-0.0,i=nincoming+1,j),
+     &           (r2min(i,j),i=j+1,nexternal)
+            do i=j+1,nexternal
+               r2min(i,j)=r2min(i,j)*dabs(r2min(i,j))    !Since r2 returns distance squared
+               r2max(i,j)=r2max(i,j)*dabs(r2max(i,j))
+            enddo
+         enddo
+         do j=nincoming+1,nexternal-1
+            write(*,formstr) 's min #',j,'>',
+     &           (s_min(i,j),i=nincoming+1,nexternal)
+         enddo
+         do j=nincoming+1,nexternal-1
+            write(*,formstr) 'xqcutij #',j,'>',
+     &           (xqcutij(i,j),i=nincoming+1,nexternal)
+         enddo
+
+cc
+cc     Set the strong coupling
+cc
+c         call set_ren_scale(P,scale)
+c
+cc     Check that the user funtions for setting the scales
+cc     have been edited if the choice of an event-by-event
+cc     scale choice has been made 
+c
+c         if(.not.fixed_ren_scale) then
+c            if(scale.eq.0d0) then
+c               write(6,*) 
+c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
+c               write(6,*) ' Dynamical renormalization scale choice '
+c               write(6,*) ' selected but user subroutine' 
+c               write(6,*) ' set_ren_scale not edited in file:setpara.f'
+c               write(6,*) ' Switching to a fixed_ren_scale choice'
+c               write(6,*) ' with scale=zmass'
+c               scale=91.2d0
+c               write(6,*) 'scale=',scale
+c               fixed_ren_scale=.true.
+c               call set_ren_scale(P,scale)
+c            endif
+c         endif
+         
+c         if(.not.fixed_fac_scale) then
+c            call set_fac_scale(P,q2fact)
+c            if(q2fact(1).eq.0d0.or.q2fact(2).eq.0d0) then
+c               write(6,*) 
+c               write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
+c               write(6,*) ' Dynamical renormalization scale choice '
+c               write(6,*) ' selected but user subroutine' 
+c               write(6,*) ' set_fac_scale not edited in file:setpara.f'
+c               write(6,*) ' Switching to a fixed_fac_scale choice'
+c               write(6,*) ' with q2fact(i)=zmass**2'
+c               fixed_fac_scale=.true.
+c               q2fact(1)=91.2d0**2
+c               q2fact(2)=91.2d0**2
+c               write(6,*) 'scales=',q2fact(1),q2fact(2)
+c            endif
+c         endif
+
+         if(fixed_ren_scale) then
+            G = SQRT(4d0*PI*ALPHAS(scale))
+            call update_as_param()
+         endif
+
+c     Put momenta in the common block to zero to start
+         do i=0,3
+            do j=1,max_particles
+               pp(i,j) = 0d0
+            enddo
+         enddo
+         
+      ENDIF ! IF FIRSTTIME
+
+      IF (CUTSDONE) THEN
+         PASSCUTS=CUTSPASSED
+         RETURN
+      ENDIF
+      CUTSDONE=.TRUE.
+c      CUTSPASSED=.FALSE.
+
+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
+      
+c
+c     Limit S_hat
+c
+c      if (x1*x2*stot .gt. 500**2) then
+c         passcuts=.false.
+c         return
+c      endif
+
+C $B$ DESACTIVATE_CUT $E$ !This is a tag for MadWeight
+
+      if(debug) write (*,*) '============================='
+      if(debug) write (*,*) ' EVENT STARTS TO BE CHECKED  '
+      if(debug) write (*,*) '============================='
+c     
+c     p_t min & max cuts
+c     
+      do i=nincoming+1,nexternal
+         if(debug) write (*,*) 'pt(',i,')=',pt(p(0,i)),'   ',etmin(i),
+     $        ':',etmax(i)
+         notgood=(pt(p(0,i)) .lt. etmin(i)).or.
+     &        (etmax(i).ge.0d0.and.pt(p(0,i)) .gt. etmax(i))
+         if (notgood) then
+            if(debug) write (*,*) i,' -> fails'
+            passcuts=.false.
+            return
+         endif
+      enddo
+c
+c    missing ET min & max cut + Invariant mass of leptons and neutrino 
+c    nb: missing Et defined as the vector sum over the neutrino's pt
+c
+c-- reset ptemp(0:3)
+      do j=0,3
+         ptemp(j)=0 ! for the neutrino
+         ptemp2(j)=0 ! for the leptons
+      enddo
+c-  sum over the momenta
+      do i=nincoming+1,nexternal
+         if(is_a_nu(i)) then            
+         if(debug) write (*,*) i,' -> neutrino '
+            do j=0,3
+               ptemp(j)=ptemp(j)+p(j,i)
+            enddo
+         elseif(is_a_l(i)) then            
+         if(debug) write (*,*) i,' -> lepton '
+            do j=0,3
+               ptemp2(j)=ptemp2(j)+p(j,i)
+            enddo
+         endif
+
+      enddo
+c-  check the et
+      if(debug.and.ptemp(0).eq.0d0) write (*,*) 'No et miss in event'
+      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Et miss =',pt(ptemp(0)),'   ',misset,':',missetmax
+      if(debug.and.ptemp2(0).eq.0d0) write (*,*) 'No leptons in event'
+      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Energy of leptons =',pt(ptemp2(0))
+      if(ptemp(0).gt.0d0) then
+         notgood=(pt(ptemp(0)) .lt. misset).or.
+     &        (missetmax.ge.0d0.and.pt(ptemp(0)) .gt. missetmax)
+         if (notgood) then
+            if(debug) write (*,*) ' missing et cut -> fails'
+            passcuts=.false.
+            return
+         endif
+      endif
+      if (mmnl.gt.0d0.or.mmnlmax.ge.0d0)then
+         if(dsqrt(SumDot(ptemp,ptemp2,1d0)).lt.mmnl.or.mmnlmax.ge.0d0.and.dsqrt(SumDot(ptemp, ptemp2,1d0)).gt.mmnlmax) then
+            if(debug) write (*,*) 'lepton invariant mass -> fails'
+            passcuts=.false.
+            return
+         endif
+      endif
+c
+c     pt cut on heavy particles
+c     gives min(pt) for (at least) one heavy particle
+c
+      if(ptheavy.gt.0d0)then
+         passcuts=.false.
+         foundheavy=.false.
+         do i=nincoming+1,nexternal
+            if(is_heavy(i)) then            
+               if(debug) write (*,*) i,' -> heavy '
+               foundheavy=.true.
+               if(pt(p(0,i)).gt.ptheavy) passcuts=.true.
+            endif
+         enddo
+         
+         if(.not.passcuts.and.foundheavy)then
+            if(debug) write (*,*) ' heavy particle cut -> fails'
+            return
+         else
+            passcuts=.true.
+         endif
+      endif
+c     
+c     E min & max cuts
+c     
+      do i=nincoming+1,nexternal
+         if(debug) write (*,*) 'p(0,',i,')=',p(0,i),'   ',emin(i),':',emax(i)
+         notgood=(p(0,i) .le. emin(i)).or.
+     &        (emax(i).ge.0d0 .and. p(0,i) .gt. emax(i))
+         if (notgood) then
+            if(debug) write (*,*) i,' -> fails'
+            passcuts=.false.
+            return
+         endif
+      enddo
+c     
+c     Rapidity  min & max cuts
+c     
+      do i=nincoming+1,nexternal
+         if(debug) write (*,*) 'abs(rap(',i,'))=',abs(rap(p(0,i))),'   ',etamin(i),':',etamax(i)
+         notgood=(etamax(i).ge.0.and.abs(rap(p(0,i))) .gt. etamax(i)).or.
+     &        (abs(rap(p(0,i))) .lt. etamin(i))
+         if (notgood) then
+            if(debug) write (*,*) i,' -> fails'
+            passcuts=.false.
+            return
+         endif
+      enddo
+c     
+c     DeltaR min & max cuts
+c     
+      do i=nincoming+1,nexternal-1
+         do j=i+1,nexternal
+            if(debug) write (*,*) 'r2(',i, ',' ,j,')=',dsqrt(r2(p(0,i),p(0,j)))
+            if(debug) write (*,*) dsqrt(r2min(j,i)),dsqrt(r2max(j,i))
+            if(r2min(j,i).gt.0.or.r2max(j,i).ge.0d0) then
+               tmp=r2(p(0,i),p(0,j))
+               notgood=(tmp .lt. r2min(j,i)).or.
+     $              (r2max(j,i).ge.0d0 .and. tmp .gt. r2max(j,i))
+               if (notgood) then
+                  if(debug) write (*,*) i,j,' -> fails'
+                  passcuts=.false.
+                  return
+               endif
+            endif
+         enddo
+      enddo
+
+
+c     s-channel min & max pt of sum of 4-momenta
+c     
+      do i=nincoming+1,nexternal-1
+         do j=i+1,nexternal
+            if(debug)write (*,*) 'ptll(',i,',',j,')=',PtDot(p(0,i),p(0,j))
+            if(debug)write (*,*) ptll_min(j,i),ptll_max(j,i)
+            if(ptll_min(j,i).gt.0.or.ptll_max(j,i).ge.0d0) then
+               tmp=PtDot(p(0,i),p(0,j))
+               notgood=(tmp .lt. ptll_min(j,i).or.
+     $              ptll_max(j,i).ge.0d0 .and. tmp.gt.ptll_max(j,i))
+               if (notgood) then
+                  if(debug) write (*,*) i,j,' -> fails'
+                  passcuts=.false.
+                  return
+               endif
+            endif
+         enddo
+      enddo
+
+
+
+
+c     
+c     s-channel min & max invariant mass cuts
+c     
+      do i=nincoming+1,nexternal-1
+         do j=i+1,nexternal
+            if(debug) write (*,*) 's(',i,',',j,')=',Sumdot(p(0,i),p(0,j),+1d0)
+            if(debug) write (*,*) s_min(j,i),s_max(j,i)
+            if(s_min(j,i).gt.0.or.s_max(j,i).ge.0d0) then
+               tmp=SumDot(p(0,i),p(0,j),+1d0)
+               if(s_min(j,i).le.s_max(j,i) .or. s_max(j,i).lt.0d0)then
+                  notgood=(tmp .lt. s_min(j,i).or.
+     $                 s_max(j,i).ge.0d0 .and. tmp .gt. s_max(j,i)) 
+                  if (notgood) then
+                     if(debug) write (*,*) i,j,' -> fails'
+                     passcuts=.false.
+                     return
+                  endif
+               else
+                  notgood=(tmp .lt. s_min(j,i).and.tmp .gt. s_max(j,i)) 
+                  if (notgood) then
+                     if(debug) write (*,*) i,j,' -> fails'
+                     passcuts=.false.
+                     return
+                  endif
+               endif
+            endif
+         enddo
+      enddo
+C     $B$DESACTIVATE_BW_CUT$B$ This is a Tag for MadWeight
+c     
+c     B.W. phase space cuts
+c     
+      pass_bw=cut_bw(p)
+c     JA 4/8/11 always check pass_bw
+      if ( pass_bw.and..not.CUTSPASSED) then
+         passcuts=.false.
+         if(debug) write (*,*) ' pass_bw -> fails'
+         return
+      endif
+C     $E$DESACTIVATE_BW_CUT$E$ This is a Tag for MadWeight
+        CUTSPASSED=.FALSE.
+C     
+C     maximal and minimal pt of the jets sorted by pt
+c     
+      njets=0
+      nheavyjets=0
+
+c- fill ptjet with the pt's of the jets.
+      do i=nincoming+1,nexternal
+         if(is_a_j(i)) then
+            njets=njets+1
+            ptjet(njets)=pt(p(0,i))
+         endif
+         if(is_a_b(i)) then
+            nheavyjets=nheavyjets+1
+            ptheavyjet(nheavyjets)=pt(p(0,i))
+         endif
+
+      enddo
+      if(debug) write (*,*) 'not yet ordered ',njets,'   ',ptjet
+
+C----------------------------------------------------------------------------
+C     DURHAM_KT CUT
+C----------------------------------------------------------------------------
+
+      IF ( KT_DURHAM .GT. 0D0) THEN
+
+C       RESET JET MOMENTA
+        njets=0
+        DO I=1,NEXTERNAL
+          DO J=0,3
+            PJET(I,J) = 0D0
+          ENDDO
+        ENDDO
+
+        do i=nincoming+1,nexternal
+           if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) ) then
+             njets=njets+1
+             DO J=0,3
+               PJET(NJETS,J) = P(J,I)
+             ENDDO
+           endif
+        enddo
+
+C       COUNT NUMBER OF MASSLESS OUTGOING PARTICLES, SINCE WE DO NOT WANT
+C       TO APPLY A CUT FOR JUST A SINGLE MASSIVE PARTICLE IN THE FINAL STATE.
+        NMASSLESS = 0
+        DO I=NINCOMING+1,NEXTERNAL
+          if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and. 
+     &                          DSQRT(DOT( P(0,I), P(0,I))) .LT. 5D0) THEN
+            NMASSLESS = NMASSLESS + 1
+          ENDIF
+        ENDDO
+
+C       DURHAM KT SEPARATION CUT
+        IF(NJETS.GT.0 .AND. NMASSLESS .GT. 0) THEN
+
+        PTMIN = EBEAM(1) + EBEAM(2)
+
+        DO I=1,NJETS
+
+C         PT WITH RESPECT TO Z AXIS FOR HADRONIC COLLISIONS
+          IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
+            PT1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2)
+            PTMIN = MIN( PTMIN, PT1 )
+          ENDIF
+
+          DO J=I+1,NJETS
+C           GET ANGLE BETWEEN JETS
+            PABS1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2 + PJET(I,3)**2)
+            PABS2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2 + PJET(J,3)**2)
+C           CHECK IF 3-MOMENTA DO NOT VANISH
+            IF(PABS1*PABS2 .NE. 0D0) THEN
+              COSTH = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) + PJET(I,3)*PJET(J,3) )/(PABS1*PABS2)
+            ELSE
+C           IF 3-MOMENTA VANISH, MAKE JET COSTH = 1D0 SO THAT JET MEASURE VANISHES
+              COSTH = 1D0
+            ENDIF
+C           GET PT AND ETA OF JETS
+            PT2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2)
+            ETA1 = 0.5D0*LOG( (PJET(I,0) + PJET(I,3)) / (PJET(I,0) - PJET(I,3)) )
+            ETA2 = 0.5D0*LOG( (PJET(J,0) + PJET(J,3)) / (PJET(J,0) - PJET(J,3)) )
+C           GET COSH OF DELTA ETA, COS OF DELTA PHI
+            COSH_DETA = DCOSH( ETA1 - ETA2 )
+            COS_DPHI = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) ) / (PT1*PT2)
+            DPHI = DACOS( COS_DPHI )
+            IF ( (LPP(1).EQ.0) .AND. (LPP(2).EQ.0)) THEN
+C             KT FOR E+E- COLLISION
+              PTNOW = DSQRT( 2D0*MIN(PJET(I,0)**2,PJET(J,0)**2)*( 1D0-COSTH ) )
+             ELSE
+C             HADRONIC KT, FASTJET DEFINITION
+              PTNOW = DSQRT( MIN(PT1**2,PT2**2)*( (ETA1 - ETA2 )**2 + DPHI**2 )/(D_PARAMETER**2) )
+            ENDIF
+
+            PTMIN = MIN( PTMIN, PTNOW )
+
+          ENDDO ! LOOP OVER NJET
+
+        ENDDO ! LOOP OVER NJET
+
+C       CHECK COMPATIBILITY WITH CUT
+        IF( (PTMIN .LT. KT_DURHAM)) THEN
+          PASSCUTS = .FALSE.
+          RETURN
+        ENDIF
+
+        ENDIF ! IF NJETS.GT. 0
+
+      ENDIF ! KT_DURHAM .GT. 0D0
+
+C----------------------------------------------------------------------------
+C     PTLUND CUT
+C----------------------------------------------------------------------------
+
+      IF(PT_LUND .GT. 0D0 ) THEN
+
+C       Reset jet momenta
+        NJETS=0
+        DO I=1,NEXTERNAL
+          JETFLAVOUR(I) = 0
+          DO J=0,3
+            PJET(I,J) = 0D0
+          ENDDO
+        ENDDO
+
+C       Fill incoming particle momenta
+        DO I=1,NINCOMING
+          INCFLAVOUR(I) = IDUP(I,1,1)
+          DO J=0,3
+            PINC(I,J) = P(J,I)
+          ENDDO
+        ENDDO
+
+C       Fill final jet momenta
+        DO I=NINCOMING+1,NEXTERNAL
+          if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) ) then
+            NJETS=NJETS+1
+            JETFLAVOUR(NJETS) = IDUP(I,1,1)
+            DO J=0,3
+              PJET(NJETS,J) = P(J,I)
+            ENDDO
+          ENDIF
+        ENDDO
+
+C       PROCESS WITH EXACTLY TWO MASSLESS OUTGOING PARTICLES IS SPECIAL
+C       BECAUSE AN ENERGY-SHARING VARIABLE LIKE "Z" DOES NOT MAKE SENSE.
+C       IN THIS CASE, ONLY APPLY MINIMAL pT W.R.T BEAM CUT.
+C       THIS CUT WILL ONLY APPLY TO THE TWO-MASSLESS PARTICLE STATE.
+        NMASSLESS = 0
+        DO I=NINCOMING+1,NEXTERNAL
+          if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and. 
+     &                          DSQRT(DOT( P(0,I), P(0,I))) .LT. 5D0) THEN
+            NMASSLESS = NMASSLESS + 1
+          ENDIF
+        ENDDO
+        IF (NMASSLESS .EQ. 2 .AND. NJETS .EQ. 2 .AND.
+     &                                     NEXTERNAL-NINCOMING .EQ. 2) THEN
+          PTMINSAVE = EBEAM(1) + EBEAM(2)
+          DO I=NINCOMING+1,NEXTERNAL
+            if( .not. from_decay(I) ) then  
+              PTMINSAVE = MIN(PTMINSAVE, PT(p(0,i)))
+            ENDIF
+          ENDDO
+C         CHECK COMPATIBILITY WITH CUT
+          IF ( ((LPP(1).NE.0) .OR. (LPP(2).NE.0)) .AND.
+     &                                         PTMINSAVE .LT. PT_LUND) THEN
+            PASSCUTS = .FALSE.
+            RETURN
+          ENDIF
+C         RESET NJETS TO AVOID FURTHER MERGING SCALE CUT.
+          NJETS=0
+        ENDIF
+
+C       PYTHIA PT SEPARATION CUT
+        IF(NJETS.GT.0 .AND. NMASSLESS .GT. 0) THEN
+
+        PTMINSAVE = EBEAM(1) + EBEAM(2)
+
+        DO I=1,NJETS
+
+          PTMIN = EBEAM(1) + EBEAM(2)
+          PTMINSAVE = MIN(PTMIN,PTMINSAVE)
+
+C         Compute pythia ISR separation between i-jet and incoming.
+C         Only SM-like emissions off the beam are possible.
+          IF ( ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) .AND. 
+     &                               ABS(JETFLAVOUR(I)) .LT. 30 ) THEN
+C           Check separation to first incoming particle
+            DO L=0,3
+              PRADTEMP(L) = PINC(1,L)
+              PEMTTEMP(L) = PJET(I,L)
+              PRECTEMP(L) = PINC(2,L)
+            ENDDO
+            PT1 = RHOPYTHIA(PRADTEMP, PEMTTEMP, PRECTEMP, INCFLAVOUR(1),
+     &                                            JETFLAVOUR(I), -1, -1)
+            PTMIN = MIN( PTMIN, PT1 )
+C           Check separation to second incoming particle
+            DO L=0,3
+              PRADTEMP(L) = PINC(2,L)
+              PEMTTEMP(L) = PJET(I,L)
+              PRECTEMP(L) = PINC(1,L)
+            ENDDO
+            PT2 = RHOPYTHIA(PRADTEMP, PEMTTEMP, PRECTEMP, INCFLAVOUR(2),
+     &                                            JETFLAVOUR(I), -1, -1)
+            PTMIN = MIN( PTMIN, PT2 )
+          ENDIF
+
+C         Compute pythia FSR separation between two jets,
+C         without any knowledge of colour connections
+          DO J=1,NJETS
+            DO K=1,NJETS
+              IF ( I .NE. J .AND. I .NE. K .AND. J .NE. K ) THEN
+
+C               Check separation between final partons i and j, with k as spectator
+                DO L=0,3
+                  PRADTEMP(L) = PJET(J,L)
+                  PEMTTEMP(L) = PJET(I,L)
+                  PRECTEMP(L) = PJET(K,L)
+                ENDDO
+
+                TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
+     &                               JETFLAVOUR(J), JETFLAVOUR(I), 1, 1)
+C               Only SM-like emissions off the beam are possible, no additional
+C               BSM particles will be produced as as shower emissions.
+                IF ( ABS(JETFLAVOUR(I)) .LT. 30 ) THEN
+                  PTMIN = MIN(PTMIN, TEMP);
+                ENDIF
+
+                TEMP = RHOPYTHIA( PEMTTEMP, PRADTEMP, PRECTEMP,
+     &                               JETFLAVOUR(I), JETFLAVOUR(J), 1, 1)
+C               Only SM-like emissions off the beam are possible, no additional
+C               BSM particles will be produced as as shower emissions.
+                IF ( ABS(JETFLAVOUR(J)) .LT. 30 ) THEN
+                  PTMIN = MIN(PTMIN, TEMP);
+                ENDIF
+
+              ENDIF
+
+            ENDDO ! LOOP OVER NJET
+          ENDDO ! LOOP OVER NJET
+
+C         Compute pythia FSR separation between two jets, with initial spectator
+          IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
+            DO J=1,NJETS
+
+C             BSM particles can only be radiators, and will not be produced
+C             as shower emissions.
+              IF ( ABS(JETFLAVOUR(I)) .GT. 1000000 ) THEN
+                EXIT
+              ENDIF
+
+C             Allow both initial partons as recoiler
+              IF ( I .NE. J ) THEN
+
+C               Check with first initial as recoiler
+                DO L=0,3
+                  PRADTEMP(L) = PJET(J,L)
+                  PEMTTEMP(L) = PJET(I,L)
+                  PRECTEMP(L) = PINC(1,L)
+                ENDDO
+                TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP, 
+     &                             JETFLAVOUR(J), JETFLAVOUR(I), 1, -1);
+                IF ( LPP(1) .NE. 0 ) THEN
+                  PTMIN = MIN(PTMIN, TEMP);
+                ENDIF
+                DO L=0,3
+                  PRADTEMP(L) = PJET(J,L)
+                  PEMTTEMP(L) = PJET(I,L)
+                  PRECTEMP(L) = PINC(2,L)
+                ENDDO
+                TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
+     &                             JETFLAVOUR(J), JETFLAVOUR(I), 1, -1);
+                IF ( LPP(2) .NE. 0 ) THEN
+                  PTMIN = MIN(PTMIN, TEMP);
+                ENDIF
+              ENDIF
+            ENDDO ! LOOP OVER NJET
+          ENDIF
+
+          PTMINSAVE = MIN(PTMIN,PTMINSAVE)
+
+        ENDDO ! LOOP OVER NJET
+
+C       CHECK COMPATIBILITY WITH CUT
+        IF (PTMINSAVE .LT. PT_LUND) THEN
+          PASSCUTS = .FALSE.
+          RETURN
+        ENDIF
+ 
+        ENDIF ! IF NJETS.GT. 0 
+
+      ENDIF ! PT_LUND .GT. 0D0 
+
+C----------------------------------------------------------------------------
+C----------------------------------------------------------------------------
+
+
+
+
+c- check existance of jets if jet cuts are on
+      if(njets.lt.1.and.(htjmin.gt.0.or.ptj1min.gt.0).or.
+     $     njets.lt.2.and.ptj2min.gt.0.or.
+     $     njets.lt.3.and.ptj3min.gt.0.or.
+     $     njets.lt.4.and.ptj4min.gt.0)then
+         if(debug) write (*,*) i, ' too few jets -> fails'
+         passcuts=.false.
+         return
+      endif
+
+c - sort jet pts
+      do i=1,njets-1
+         do j=i+1,njets
+            if(ptjet(j).gt.ptjet(i)) then
+               temp=ptjet(i)
+               ptjet(i)=ptjet(j)
+               ptjet(j)=temp
+            endif
+         enddo
+      enddo
+      if(debug) write (*,*) 'ordered ',njets,'   ',ptjet
+c
+c     Use "and" or "or" prescriptions 
+c     
+      inclht=0
+
+      if(njets.gt.0) then
+
+       notgood=.not.jetor
+       if(debug) write (*,*) 'jetor :',jetor  
+       if(debug) write (*,*) '0',notgood   
+      
+      do i=1,min(njets,4)
+            if(debug) write (*,*) i,ptjet(i), '   ',ptjmin4(i),
+     $        ':',ptjmax4(i)
+         if(jetor) then     
+c---  if one of the jets does not pass, the event is rejected
+            notgood=notgood.or.(ptjmax4(i).ge.0d0 .and.
+     $           ptjet(i).gt.ptjmax4(i)) .or.
+     $           (ptjet(i).lt.ptjmin4(i))
+            if(debug) write (*,*) i,' notgood total:', notgood   
+         else
+c---  all cuts must fail to reject the event
+            notgood=notgood.and.(ptjmax4(i).ge.0d0 .and.
+     $              ptjet(i).gt.ptjmax4(i) .or.
+     $              (ptjet(i).lt.ptjmin4(i)))
+            if(debug) write (*,*) i,' notgood total:', notgood   
+         endif
+      enddo
+
+
+      if (notgood) then
+         if(debug) write (*,*) i, ' multiple pt -> fails'
+         passcuts=.false.
+         return
+      endif
+
+c---------------------------
+c      Ht cuts
+C---------------------------
+      htj=ptjet(1)
+
+      do i=2,njets
+         htj=htj+ptjet(i)
+         if(debug) write (*,*) i, 'htj ',htj
+         if(debug.and.i.le.4) write (*,*) 'htmin ',i,' ', htjmin4(i),':',htjmax4(i)
+         if(i.le.4)then
+            if(htj.lt.htjmin4(i) .or.
+     $        htjmax4(i).ge.0d0.and.htj.gt.htjmax4(i)) then
+            if(debug) write (*,*) i, ' ht -> fails'
+            passcuts=.false.
+            return
+            endif
+         endif
+      enddo
+
+      if(htj.lt.htjmin.or.htjmax.ge.0d0.and.htj.gt.htjmax)then
+         if(debug) write (*,*) i, ' htj -> fails'
+         passcuts=.false.
+         return
+      endif
+
+      inclht=htj
+
+      endif !if there are jets 
+
+      if(nheavyjets.gt.0) then
+         do i=1,nheavyjets
+            inclht=inclht+ptheavyjet(i)
+         enddo
+      endif !if there are heavyjets
+
+      if(inclht.lt.inclHtmin.or.
+     $     inclHtmax.ge.0d0.and.inclht.gt.inclHtmax)then
+         if(debug) write (*,*) ' inclhtmin=',inclHtmin,' -> fails'
+         passcuts=.false.
+         return
+      endif
+
+C     
+C     maximal and minimal pt of the leptons sorted by pt
+c     
+      nleptons=0
+
+      if(ptl1min.gt.0.or.ptl2min.gt.0.or.ptl3min.gt.0.or.ptl4min.gt.0.or.
+     $     ptl1max.ge.0d0.or.ptl2max.ge.0d0.or.
+     $     ptl3max.ge.0d0.or.ptl4max.ge.0d0) then
+
+c     - fill ptlepton with the pt's of the leptons.
+         do i=nincoming+1,nexternal
+            if(is_a_l(i)) then
+               nleptons=nleptons+1
+               ptlepton(nleptons)=pt(p(0,i))
+            endif
+         enddo
+         if(debug) write (*,*) 'not yet ordered ',njets,'   ',ptjet
+
+c     - check existance of leptons if lepton cuts are on
+         if(nleptons.lt.1.and.ptl1min.gt.0.or.
+     $        nleptons.lt.2.and.ptl2min.gt.0.or.
+     $        nleptons.lt.3.and.ptl3min.gt.0.or.
+     $        nleptons.lt.4.and.ptl4min.gt.0)then
+            if(debug) write (*,*) i, ' too few leptons -> fails'
+            passcuts=.false.
+            return
+         endif
+
+c     - sort lepton pts
+         do i=1,nleptons-1
+            do j=i+1,nleptons
+               if(ptlepton(j).gt.ptlepton(i)) then
+                  temp=ptlepton(i)
+                  ptlepton(i)=ptlepton(j)
+                  ptlepton(j)=temp
+               endif
+            enddo
+         enddo
+         if(debug) write (*,*) 'ordered ',nleptons,'   ',ptlepton
+
+         if(nleptons.gt.0) then
+
+            notgood = .false.
+            do i=1,min(nleptons,4)
+               if(debug) write (*,*) i,ptlepton(i), '   ',ptlmin4(i),':',ptlmax4(i)
+c---  if one of the leptons does not pass, the event is rejected
+               notgood=notgood.or.
+     $              (ptlmax4(i).ge.0d0.and.ptlepton(i).gt.ptlmax4(i)).or.
+     $              (ptlepton(i).lt.ptlmin4(i))
+               if(debug) write (*,*) i,' notgood total:', notgood   
+            enddo
+
+
+            if (notgood) then
+               if(debug) write (*,*) i, ' multiple pt -> fails'
+               passcuts=.false.
+               return
+            endif
+         endif
+      endif
+C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>     
+C     SPECIAL CUTS
+C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+C     REQUIRE AT LEAST ONE JET WITH PT>XPTJ
+         
+         IF(xptj.gt.0d0) THEN
+            xvar=0
+            do i=nincoming+1,nexternal
+               if(is_a_j(i)) xvar=max(xvar,pt(p(0,i)))
+            enddo
+            if (xvar .lt. xptj) then
+               passcuts=.false.
+               return
+            endif
+         ENDIF
+
+C     REQUIRE AT LEAST ONE PHOTON WITH PT>XPTA
+         
+         IF(xpta.gt.0d0) THEN
+            xvar=0
+            do i=nincoming+1,nexternal
+               if(is_a_a(i)) xvar=max(xvar,pt(p(0,i)))
+            enddo
+            if (xvar .lt. xpta) then
+               passcuts=.false.
+               return
+            endif
+         ENDIF
+
+C     REQUIRE AT LEAST ONE B  WITH PT>XPTB
+         
+         IF(xptb.gt.0d0) THEN
+            xvar=0
+            do i=nincoming+1,nexternal
+               if(is_a_b(i)) xvar=max(xvar,pt(p(0,i)))
+            enddo
+            if (xvar .lt. xptb) then
+               passcuts=.false.
+               return
+            endif
+         ENDIF
+
+C     REQUIRE AT LEAST ONE LEPTON  WITH PT>XPTL
+         
+         IF(xptl.gt.0d0) THEN
+            xvar=0
+            do i=nincoming+1,nexternal
+               if(is_a_l(i)) xvar=max(xvar,pt(p(0,i)))
+            enddo
+            if (xvar .lt. xptl) then
+               passcuts=.false.
+               if(debug) write (*,*) ' xptl -> fails'
+               return
+            endif
+         ENDIF
+C
+C     WBF CUTS: TWO TYPES
+C    
+C     FIRST TYPE:  implemented by FM
+C
+C     1. FIND THE 2 HARDEST JETS
+C     2. REQUIRE |RAP(J)|>XETAMIN
+C     3. REQUIRE RAP(J1)*ETA(J2)<0
+C
+C     SECOND TYPE : added by Simon de Visscher 1-08-2007
+C
+C     1. FIND THE 2 HARDEST JETS
+C     2. REQUIRE |RAP(J1)-RAP(J2)|>DELTAETA
+C     3. REQUIRE RAP(J1)*RAP(J2)<0
+C
+C
+         hardj1=0
+         hardj2=0
+         ptmax1=0d0
+         ptmax2=0d0
+
+C-- START IF AT LEAST ONE OF THE CUTS IS ACTIVATED
+         
+         IF(XETAMIN.GT.0D0.OR.DELTAETA.GT.0D0) THEN
+            
+C-- FIND THE HARDEST JETS
+
+            do i=nincoming+1,nexternal
+               if(is_a_j(i)) then
+c                  write (*,*) i,pt(p(0,i))
+                  if(pt(p(0,i)).gt.ptmax1) then
+                     hardj2=hardj1
+                     ptmax2=ptmax1
+                     hardj1=i
+                     ptmax1=pt(p(0,i))
+                  elseif(pt(p(0,i)).gt.ptmax2) then
+                     hardj2=i
+                     ptmax2=pt(p(0,i))
+                  endif
+c                  write (*,*) hardj1,hardj2,ptmax1,ptmax2
+               endif
+            enddo
+            
+            if (hardj2.eq.0) goto 21 ! bypass vbf cut since not enough jets
+
+C-- NOW APPLY THE CUT I            
+
+            if (abs(rap(p(0,hardj1))) .lt. xetamin
+     &       .or.abs(rap(p(0,hardj2))) .lt. xetamin
+     &       .or.rap(p(0,hardj1))*rap(p(0,hardj2)) .gt.0d0) then
+             passcuts=.false.
+             return
+            endif
+
+            
+C-- NOW APPLY THE CUT II
+            
+            if (abs(rap(p(0,hardj1))-rap(p(0,hardj2))) .lt. deltaeta) then
+             passcuts=.false.
+             return
+            endif
+         
+c            write (*,*) hardj1,hardj2,rap(p(0,hardj1)),rap(p(0,hardj2))
+         
+         ENDIF
+
+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 If we do not require a mimimum jet energy, there's no need to apply
+c jet clustering and all that.
+      if (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
+c the number of (light) QCD partons at the real-emission level (i.e. one
+c more than the Born).
+
+      nQCD=0
+      do j=nincoming+1,nexternal
+         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
+            enddo
+         endif
+      enddo
+
+        nph=0
+        do j=nincoming+1,nexternal
+          if (is_a_a(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_l(j)) then
+              nem=nem+1
+              do i=0,3
+                pem(i,nem)=p(i,j) ! Use C.o.M. frame momenta
+              enddo
+            endif
+          enddo
+        endif
+
+        alliso=.true.
+
+        j=0
+        dowhile(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
+
+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)))
+            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
+      endif
+
+ 444    continue
+c End photon isolation
+
+
+C...Set couplings if event passed cuts
+
+ 21   if(.not.fixed_ren_scale) then
+         call set_ren_scale(P,scale)
+         if(scale.gt.0) G = SQRT(4d0*PI*ALPHAS(scale))
+      endif
+
+      if(.not.fixed_fac_scale) then
+         call set_fac_scale(P,q2fact)
+      endif
+
+c
+c     Here we cluster event and reset factorization and renormalization
+c     scales on an event-by-event basis, as well as check xqcut for jets
+c
+c     Note the following condition is the first line of setclscales
+c      if(xqcut.gt.0d0.or.ickkw.gt.0.or.scale.eq.0.or.q2fact(1).eq.0)then
+c     Do not duplicate it since some variable are set for syscalc in the fct
+        if(.not.setclscales(p))then
+           cutsdone=.false.
+           cutspassed=.false.
+           passcuts = .false.
+           if(debug) write (*,*) 'setclscales -> fails'
+           return
+       endif
+c      endif
+
+c     Set couplings in model files
+      if(.not.fixed_ren_scale.or..not.fixed_couplings) then
+         if (.not.fixed_couplings)then
+            do i=0,3
+               do j=1,nexternal
+                  pp(i,j)=p(i,j)
+               enddo
+            enddo
+         endif
+         call update_as_param()
+      endif
+
+      IF (FIRSTTIME2) THEN
+        FIRSTTIME2=.FALSE.
+        write(6,*) 'alpha_s for scale ',scale,' is ', G**2/(16d0*atan(1d0))
+      ENDIF
+
+      if(debug) write (*,*) '============================='
+      if(debug) write (*,*) ' EVENT PASSED THE CUTS       '
+      if(debug) write (*,*) '============================='
+
+      CUTSPASSED=.TRUE.
+
+      RETURN
+      END
+
+
+C
+C     FUNCTION FOR ISOLATION
+C
+
+      function iso_getdrv40(p1,p2)
+      implicit none
+      real*8 iso_getdrv40,p1(0:3),p2(0:3)
+      real*8 iso_getdr
+c
+      iso_getdrv40=iso_getdr(p1(0),p1(1),p1(2),p1(3),
+     #                       p2(0),p2(1),p2(2),p2(3))
+      return
+      end
+
+
+      function iso_getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
+      implicit none
+      real*8 iso_getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
+     # iso_getpseudorap,iso_getdelphi
+c
+      deta=iso_getpseudorap(en1,ptx1,pty1,pl1)-
+     #     iso_getpseudorap(en2,ptx2,pty2,pl2)
+      dphi=iso_getdelphi(ptx1,pty1,ptx2,pty2)
+      iso_getdr=sqrt(dphi**2+deta**2)
+      return
+      end
+
+
+      function iso_getpseudorap(en,ptx,pty,pl)
+      implicit none
+      real*8 iso_getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
+      parameter (tiny=1.d-5)
+c
+      pt=sqrt(ptx**2+pty**2)
+      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
+        eta=sign(1.d0,pl)*1.d8
+      else
+        th=atan2(pt,pl)
+        eta=-log(tan(th/2.d0))
+      endif
+      iso_getpseudorap=eta
+      return
+      end
+
+
+      function iso_getdelphi(ptx1,pty1,ptx2,pty2)
+      implicit none
+      real*8 iso_getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
+      parameter (tiny=1.d-5)
+c
+      pt1=sqrt(ptx1**2+pty1**2)
+      pt2=sqrt(ptx2**2+pty2**2)
+      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
+        tmp=ptx1*ptx2+pty1*pty2
+        tmp=tmp/(pt1*pt2)
+        if(abs(tmp).gt.1.d0+tiny)then
+          write(*,*)'Cosine larger than 1'
+          stop
+        elseif(abs(tmp).ge.1.d0)then
+          tmp=sign(1.d0,tmp)
+        endif
+        tmp=acos(tmp)
+      else
+        tmp=1.d8
+      endif
+      iso_getdelphi=tmp
+      return
+      end
+
+      function chi_gamma_iso(dr,R0,xn,epsgamma,pTgamma)
+c Eq.(3.4) of Phys.Lett. B429 (1998) 369-374 [hep-ph/9801442]
+      implicit none
+      real*8 chi_gamma_iso,dr,R0,xn,epsgamma,pTgamma
+      real*8 tmp,axn
+c
+      axn=abs(xn)
+      tmp=epsgamma*pTgamma
+      if(axn.ne.0.d0)then
+        tmp=tmp*( (1-cos(dr))/(1-cos(R0)) )**axn
+      endif
+      chi_gamma_iso=tmp
+      return
+      end
+
+
+*
+* $Id: sortzv.F,v 1.1.1.1 1996/02/15 17:49:50 mclareni Exp $
+*
+* $Log: sortzv.F,v $
+* Revision 1.1.1.1  1996/02/15 17:49:50  mclareni
+* Kernlib
+*
+*
+c$$$#include "kerngen/pilot.h"
+      SUBROUTINE SORTZV (A,INDEX,N1,MODE,NWAY,NSORT)
+C
+C CERN PROGLIB# M101    SORTZV          .VERSION KERNFOR  3.15  820113
+C ORIG. 02/10/75
+C
+      DIMENSION A(N1),INDEX(N1)
+C
+C
+      N = N1
+      IF (N.LE.0)            RETURN
+      IF (NSORT.NE.0) GO TO 2
+      DO 1 I=1,N
+    1 INDEX(I)=I
+C
+    2 IF (N.EQ.1)            RETURN
+      IF (MODE)    10,20,30
+   10 CALL SORTTI (A,INDEX,N)
+      GO TO 40
+C
+   20 CALL SORTTC(A,INDEX,N)
+      GO TO 40
+C
+   30 CALL SORTTF (A,INDEX,N)
+C
+   40 IF (NWAY.EQ.0) GO TO 50
+      N2 = N/2
+      DO 41 I=1,N2
+      ISWAP = INDEX(I)
+      K = N+1-I
+      INDEX(I) = INDEX(K)
+   41 INDEX(K) = ISWAP
+   50 RETURN
+      END
+*     ========================================
+      SUBROUTINE SORTTF (A,INDEX,N1)
+C
+      DIMENSION A(N1),INDEX(N1)
+C
+      N = N1
+      DO 3 I1=2,N
+      I3 = I1
+      I33 = INDEX(I3)
+      AI = A(I33)
+    1 I2 = I3/2
+      IF (I2) 3,3,2
+    2 I22 = INDEX(I2)
+      IF (AI.LE.A (I22)) GO TO 3
+      INDEX (I3) = I22
+      I3 = I2
+      GO TO 1
+    3 INDEX (I3) = I33
+    4 I3 = INDEX (N)
+      INDEX (N) = INDEX (1)
+      AI = A(I3)
+      N = N-1
+      IF (N-1) 12,12,5
+    5 I1 = 1
+    6 I2 = I1 + I1
+      IF (I2.LE.N) I22= INDEX(I2)
+      IF (I2-N) 7,9,11
+    7 I222 = INDEX (I2+1)
+      IF (A(I22)-A(I222)) 8,9,9
+    8 I2 = I2+1
+      I22 = I222
+    9 IF (AI-A(I22)) 10,11,11
+   10 INDEX(I1) = I22
+      I1 = I2
+      GO TO 6
+   11 INDEX (I1) = I3
+      GO TO 4
+   12 INDEX (1) = I3
+      RETURN
+      END
+*     ========================================
+      SUBROUTINE SORTTI (A,INDEX,N1)
+C
+      INTEGER A,AI
+      DIMENSION A(N1),INDEX(N1)
+C
+      N = N1
+      DO 3 I1=2,N
+      I3 = I1
+      I33 = INDEX(I3)
+      AI = A(I33)
+    1 I2 = I3/2
+      IF (I2) 3,3,2
+    2 I22 = INDEX(I2)
+      IF (AI.LE.A (I22)) GO TO 3
+      INDEX (I3) = I22
+      I3 = I2
+      GO TO 1
+    3 INDEX (I3) = I33
+    4 I3 = INDEX (N)
+      INDEX (N) = INDEX (1)
+      AI = A(I3)
+      N = N-1
+      IF (N-1) 12,12,5
+    5 I1 = 1
+    6 I2 = I1 + I1
+      IF (I2.LE.N) I22= INDEX(I2)
+      IF (I2-N) 7,9,11
+    7 I222 = INDEX (I2+1)
+      IF (A(I22)-A(I222)) 8,9,9
+    8 I2 = I2+1
+      I22 = I222
+    9 IF (AI-A(I22)) 10,11,11
+   10 INDEX(I1) = I22
+      I1 = I2
+      GO TO 6
+   11 INDEX (I1) = I3
+      GO TO 4
+   12 INDEX (1) = I3
+      RETURN
+      END
+*     ========================================
+      SUBROUTINE SORTTC (A,INDEX,N1)
+C
+      INTEGER A,AI
+      DIMENSION A(N1),INDEX(N1)
+C
+      N = N1
+      DO 3 I1=2,N
+      I3 = I1
+      I33 = INDEX(I3)
+      AI = A(I33)
+    1 I2 = I3/2
+      IF (I2) 3,3,2
+    2 I22 = INDEX(I2)
+      IF(ICMPCH(AI,A(I22)))3,3,21
+   21 INDEX (I3) = I22
+      I3 = I2
+      GO TO 1
+    3 INDEX (I3) = I33
+    4 I3 = INDEX (N)
+      INDEX (N) = INDEX (1)
+      AI = A(I3)
+      N = N-1
+      IF (N-1) 12,12,5
+    5 I1 = 1
+    6 I2 = I1 + I1
+      IF (I2.LE.N) I22= INDEX(I2)
+      IF (I2-N) 7,9,11
+    7 I222 = INDEX (I2+1)
+      IF (ICMPCH(A(I22),A(I222))) 8,9,9
+    8 I2 = I2+1
+      I22 = I222
+    9 IF (ICMPCH(AI,A(I22))) 10,11,11
+   10 INDEX(I1) = I22
+      I1 = I2
+      GO TO 6
+   11 INDEX (I1) = I3
+      GO TO 4
+   12 INDEX (1) = I3
+      RETURN
+      END
+*     ========================================
+      FUNCTION ICMPCH(IC1,IC2)
+C     FUNCTION TO COMPARE TWO 4 CHARACTER EBCDIC STRINGS - IC1,IC2
+C     ICMPCH=-1 IF HEX VALUE OF IC1 IS LESS THAN IC2
+C     ICMPCH=0  IF HEX VALUES OF IC1 AND IC2 ARE THE SAME
+C     ICMPCH=+1 IF HEX VALUES OF IC1 IS GREATER THAN IC2
+      I1=IC1
+      I2=IC2
+      IF(I1.GE.0.AND.I2.GE.0)GOTO 40
+      IF(I1.GE.0)GOTO 60
+      IF(I2.GE.0)GOTO 80
+      I1=-I1
+      I2=-I2
+      IF(I1-I2)80,70,60
+ 40   IF(I1-I2)60,70,80
+ 60   ICMPCH=-1
+      RETURN
+ 70   ICMPCH=0
+      RETURN
+ 80   ICMPCH=1
+      RETURN
+      END
+
+c************************************************************************
+c     Returns pTLund (i.e. the Pythia8 evolution variable) between two 
+c     particles with momenta prad and pemt with momentum prec as spectator
+c************************************************************************
+
+      DOUBLE PRECISION FUNCTION RHOPYTHIA(PRAD, PEMT, PREC, FLAVRAD,
+     &                                        FLAVEMT, RADTYPE, RECTYPE)
+
+       IMPLICIT NONE
+c
+c     Arguments
+c
+      DOUBLE PRECISION PRAD(0:3),PEMT(0:3), PREC(0:3)
+      INTEGER FLAVRAD, FLAVEMT, RADTYPE,RECTYPE
+c
+c     Local
+c
+      DOUBLE PRECISION Q(0:3),SUM(0:3), qBR(0:3), qAR(0:3)
+      DOUBLE PRECISION Q2, m2Rad, m2Emt, m2Dip, qBR2, qAR2, x1, x2, z
+      DOUBLE PRECISION m2RadAft, m2EmtAft, m2Rec, m2RadBef, m2ar, rescale
+      DOUBLE PRECISION TEMP, lambda13, k1, k3, m2Final
+      DOUBLE PRECISION m0u, m0d, m0c, m0s, m0t, m0b, m0w, m0z, m0x
+      DOUBLE PRECISION PRECAFT(0:3)
+      INTEGER emtsign
+      INTEGER idRadBef
+      LOGICAL allowed
+
+c-----
+c  Begin Code
+c-----
+
+C     Set masses. Currently no way of getting those?
+      m0u = 0.0
+      m0d = 0.0
+      m0c = 1.5
+      m0s = 0.0
+      m0t = 172.5
+      m0w = 80.4
+      m0z = 91.188
+      m0x = 400.0
+
+C     Store recoiler momentum (since FI splittings require recoiler
+C     rescaling)
+      PRECAFT(0) = PREC(0)
+      PRECAFT(1) = PREC(1)
+      PRECAFT(2) = PREC(2)
+      PRECAFT(3) = PREC(3)
+C     Get sign of emitted momentum
+      emtsign = 1
+      if(radtype .eq. -1) emtsign = -1
+
+C     Get virtuality
+      Q(0) = pRad(0) + emtsign*pEmt(0)
+      Q(1) = pRad(1) + emtsign*pEmt(1)
+      Q(2) = pRad(2) + emtsign*pEmt(2)
+      Q(3) = pRad(3) + emtsign*pEmt(3)
+      Q2   = emtsign * ( Q(0)**2 - Q(1)**2 - Q(2)**2 - Q(3)**2 );
+
+C     Reset allowed
+      allowed = .true.
+
+C     Splitting not possible for negative virtuality.
+      if ( Q2 .lt. 0.0 ) allowed = .false.
+
+C     Try to reconstruct flavour of radiator before emission.
+      idRadBef = 0
+C     gluon radiation: idBef = idAft
+      if (abs(flavEmt) .eq. 21 .or. abs(flavEmt) .eq. 22 ) idRadBef=flavRad
+C     final state gluon splitting: idBef = 21
+      if (radtype .eq.  1 .and. flavEmt .eq. -flavRad) idRadBef=21
+C     final state quark -> gluon conversion
+      if (radtype .eq.  1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. 21) idRadBef=flavEmt
+C     initial state gluon splitting: idBef = -idEmt
+      if (radtype .eq. -1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. 21) idRadBef=-flavEmt
+C     initial state gluon -> quark conversion
+      if (radtype .eq. -1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. flavEmt) idRadBef=21
+C     W-boson radiation
+      if (flavEmt .eq.  24) idRadBef = flavRad+1
+      if (flavEmt .eq. -24) idRadBef = flavRad-1
+
+C     Set particle masses.
+      m2RadAft = 0.0
+      m2EmtAft = 0.0
+      m2Rec    = 0.0
+      m2RadBef = 0.0
+
+      m2RadAft = pRad(0)*pRad(0)-pRad(1)*pRad(1)-pRad(2)*pRad(2)-pRad(3)*pRad(3)
+      m2EmtAft = pEmt(0)*pEmt(0)-pEmt(1)*pEmt(1)-pEmt(2)*pEmt(2)-pEmt(3)*pEmt(3)
+      m2Rec    = pRec(0)*pRec(0)-pRec(1)*pRec(1)-pRec(2)*pRec(2)-pRec(3)*pRec(3)
+
+      if (m2RadAft .lt. 1d-4) m2RadAft = 0.0
+      if (m2EmtAft .lt. 1d-4) m2EmtAft = 0.0
+      if (m2Rec    .lt. 1d-4) m2Rec    = 0.0
+
+      if (abs(flavRad) .ne. 21 .and. abs(flavRad) .ne. 22 .and.
+     &    abs(flavEmt) .ne. 24 .and.
+     &    abs(flavRad) .ne. abs(flavEmt)) then
+        m2RadBef = m2RadAft
+      else if (abs(flavEmt) .eq. 24) then
+        if (idRadBef .ne. 0) then
+          if( abs(idRadBef) .eq. 4 ) m2RadBef       = m0c**2
+          if( abs(idRadBef) .eq. 5 ) m2RadBef       = m0b**2
+          if( abs(idRadBef) .eq. 6 ) m2RadBef       = m0t**2
+          if( abs(idRadBef) .eq. 9000001 ) m2RadBef = m0x**2
+        endif
+      else if (radtype .eq. -1) then
+        if (abs(flavRad) .eq. 21 .and. abs(flavEmt) .eq. 21) m2RadBef = m2EmtAft
+      endif
+
+C     Calculate dipole mass for final-state radiation.
+      m2Final = 0.0
+      m2Final = m2Final + (pRad(0) + pRec(0) + pEmt(0))**2
+      m2Final = m2Final - (pRad(1) + pRec(1) + pEmt(1))**2
+      m2Final = m2Final - (pRad(2) + pRec(2) + pEmt(2))**2
+      m2Final = m2Final - (pRad(3) + pRec(3) + pEmt(3))**2
+
+C     Final state splitting not possible for negative dipole mass.
+      if ( radtype .eq. 1 .and. m2Final .lt. 0.0 ) allowed = .false.
+
+C     Rescale recoiler for final-intial splittings.
+      rescale = 1.0
+      if (radtype .eq. 1 .and. rectype .eq. -1) then
+        m2ar = m2Final - 2.0*Q2 + 2.0*m2RadBef
+        rescale = (1.0 - (Q2 - m2RadBef) / (m2ar-m2RadBef))
+     &           /(1.0 + (Q2 - m2RadBef) / (m2ar-m2RadBef))
+        pRecAft(0) = pRecAft(0)*rescale
+        pRecAft(1) = pRecAft(1)*rescale
+        pRecAft(2) = pRecAft(2)*rescale
+        pRecAft(3) = pRecAft(3)*rescale
+      endif
+
+C     Final-initial splitting not possible for negative rescaling.
+      if ( rescale .lt. 0.0 ) allowed = .false.
+
+C     Construct dipole momentum for FSR.
+      sum(0) = pRad(0) + pRecAft(0) + pEmt(0)
+      sum(1) = pRad(1) + pRecAft(1) + pEmt(1)
+      sum(2) = pRad(2) + pRecAft(2) + pEmt(2)
+      sum(3) = pRad(3) + pRecAft(3) + pEmt(3)
+      m2Dip = sum(0)**2 - sum(1)**2 - sum(2)**2 - sum(3)**2
+
+C     Construct 2->3 variables for FSR
+      x1     = 2. * ( sum(0)*pRad(0) - sum(1)*pRad(1)
+     &             -  sum(2)*pRad(2) - sum(3)*pRad(3) ) / m2Dip
+      x2     = 2. * ( sum(0)*pRecAft(0) - sum(1)*pRecAft(1)
+     &             -  sum(2)*pRecAft(2) - sum(3)*pRecAft(3) ) / m2Dip
+
+C     Final state splitting not possible for ill-defined
+C     3-body-variables.
+      if ( radtype .eq. 1 .and. x1 .lt. 0.0 ) allowed = .false.
+      if ( radtype .eq. 1 .and. x1 .gt. 1.0 ) allowed = .false.
+      if ( radtype .eq. 1 .and. x2 .lt. 0.0 ) allowed = .false.
+      if ( radtype .eq. 1 .and. x2 .gt. 1.0 ) allowed = .false.
+
+C     Auxiliary variables for massive FSR
+      lambda13 = DSQRT( (Q2 - m2RadAft - m2EmtAft )**2 - 4.0 * m2RadAft*m2EmtAft)
+      k1 = ( Q2 - lambda13 + (m2EmtAft - m2RadAft ) ) / ( 2.0 * Q2)
+      k3 = ( Q2 - lambda13 - (m2EmtAft - m2RadAft ) ) / ( 2.0 * Q2)
+
+C     Construct momenta of dipole before/after splitting for ISR
+      qBR(0) = pRad(0) + pRec(0) - pEmt(0)
+      qBR(1) = pRad(1) + pRec(1) - pEmt(1)
+      qBR(2) = pRad(2) + pRec(2) - pEmt(2)
+      qBR(3) = pRad(3) + pRec(3) - pEmt(3)
+      qBR2   = qBR(0)**2 - qBR(1)**2 - qBR(2)**2 - qBR(3)**2
+
+      qAR(0) = pRad(0) + pRec(0)
+      qAR(1) = pRad(1) + pRec(1)
+      qAR(2) = pRad(2) + pRec(2)
+      qAR(3) = pRad(3) + pRec(3)
+      qAR2   = qAR(0)**2 - qAR(1)**2 - qAR(2)**2 - qAR(3)**2
+
+C     Calculate z of splitting, different for FSR and ISR
+      z = 1.0 / (1.0 - k1 -k3) * ( x1 / (2.0-x2) - k3)
+      if(radtype .eq. -1 ) z = qBR2 / qAR2;
+
+C     Splitting not possible for ill-defined energy sharing.
+      if ( z .lt. 0.0 .or. z .gt. 1.0 ) allowed = .false.
+
+C     pT^2 = separation * virtuality (corrected with mass for FSR)
+      if (radtype .eq.  1) temp = z*(1-z)*(Q2 - m2RadBef)
+      if (radtype .eq. -1) temp = (1-z)*Q2
+
+C     Check threshold in ISR
+      if (radtype .ne. 1) then
+        if ((abs(flavRad) .eq. 4 .or. abs(flavEmt) .eq. 4)
+     &    .and. dsqrt(temp) .le. 2.0*m0c**2 ) temp = (1.-z)*(Q2+m0c**2)
+        if ((abs(flavRad) .eq. 5 .or. abs(flavEmt) .eq. 5)
+     &    .and. dsqrt(temp) .le. 2.0*m0b**2 ) temp = (1.-z)*(Q2+m0b**2)
+      endif
+
+C     Kinematically impossible splittings should not be included in the
+C     pT definition!
+      if( .not. allowed) temp = 1d15
+
+      if(temp .lt. 0.0) temp = 0.0
+
+C     Return pT
+      rhoPythia = dsqrt(temp);
+
+      RETURN
+      END

=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2014-11-06 01:54:56 +0000
+++ Template/LO/SubProcesses/genps.f	2016-08-19 23:07:19 +0000
@@ -1159,7 +1159,7 @@
 
 C     IF THERE IS NO S CHANNEL POLE USE BELOW:
 
-      TAUMIN = SMIN/S
+      TAUMIN = 0d0 !SMIN/S !keep scale fix
       TAUMAX = 1D0
       TAU    = (TAUMAX-TAUMIN)*X(1)+TAUMIN
       SJACOBI=  sjacobi*(TAUMAX-TAUMIN)

=== modified file 'Template/LO/SubProcesses/makefile'
--- Template/LO/SubProcesses/makefile	2016-01-28 14:24:29 +0000
+++ Template/LO/SubProcesses/makefile	2016-08-19 23:07:19 +0000
@@ -34,10 +34,10 @@
 # Binaries
 
 $(PROG): $(PROCESS) auto_dsig.o $(LIBS)
-	$(FC) $(LDFLAGS) -o $(PROG) $(PROCESS) $(LINKLIBS)
+	$(FC) -o $(PROG) $(PROCESS) $(LINKLIBS) $(LDFLAGS)
 
 gensym: $(SYMMETRY) configs.inc $(LIBDIR)libmodel.$(libext) $(LIBDIR)libgeneric.$(libext)
-	$(FC) $(LDFLAGS) -o gensym $(SYMMETRY) -L../../lib/ -lmodel -lgeneric
+	$(FC) -o gensym $(SYMMETRY) -L../../lib/ -lmodel -lgeneric $(LDFLAGS)
 
 $(LIBDIR)libmodel.$(libext): ../../Cards/param_card.dat
 	cd ../../Source/MODEL; make

=== modified file 'Template/LO/SubProcesses/myamp.f'
--- Template/LO/SubProcesses/myamp.f	2016-03-18 17:36:28 +0000
+++ Template/LO/SubProcesses/myamp.f	2016-08-19 23:07:19 +0000
@@ -211,7 +211,8 @@
 c           Here we set if the BW is "on-shell" for LesHouches
 c
             onshell = (abs(xmass - prmass(i,iconfig)) .lt.
-     $           bwcutoff*prwidth(i,iconfig))
+     $           bwcutoff*prwidth(i,iconfig).and.
+     $           prwidth(i,iconfig)/prmass(i,iconfig).lt.0.1d0)
             if(onshell)then
 c     Remove on-shell forbidden s-channels (gForceBW=2) (JA 2/10/11)
               if(gForceBW(i,iconfig).eq.2) then
@@ -605,7 +606,7 @@
 c-----------------------
          if (swidth(i) .eq. 0d0) then
             if (xo.lt.smin/stot)then
-                xo = -1d0*smin/stot
+                xo = 1d0*smin/stot
             endif
             swidth(i) = xo
             spole(i)= -2.0d0    ! 1/s pole

=== modified file 'Template/LO/SubProcesses/refine.sh'
--- Template/LO/SubProcesses/refine.sh	2015-05-22 18:00:40 +0000
+++ Template/LO/SubProcesses/refine.sh	2016-08-19 23:07:19 +0000
@@ -8,7 +8,12 @@
 if [[ -e MadLoop5_resources.tar.gz && ! -e MadLoop5_resources ]]; then
 tar -xzf MadLoop5_resources.tar.gz
 fi
-k=%(name)s_app.log
+keeplog=%(keeplog)s
+if [ "$keeplog" = true ] ; then
+    k=%(name)s_app.log
+else
+    k=/dev/null
+fi
 script=%(script_name)s                         
 
 grid_directory=%(base_directory)s
@@ -24,9 +29,10 @@
           	 fi
           fi	
      fi
-     
      cd $j
-     rm -f $k
+     if [ "$keeplog" = true ] ; then
+	 rm -f $k
+     fi
      rm -f moffset.dat >& /dev/null
       echo   %(offset)s  > moffset.dat
      if  [[ -e ftn26 ]]; then
@@ -50,17 +56,27 @@
      # filesystem problem (executable not found)
      for((try=1;try<=16;try+=1)); 
      do
-         ../madevent 2>&1 >> $k <input_sg.txt | tee -a $k;
-     status_code=${PIPESTATUS[0]};
-         if [ -s $k ]
-         then
-             break
-         else
-             echo $try > fail.log 
-         fi
+	 if [ "$keeplog" = true ] ; then
+             ../madevent 2>&1 >> $k <input_sg.txt | tee -a $k;
+	     status_code=${PIPESTATUS[0]};
+             if [ -s $k ]
+             then
+		 break
+             else
+		 echo $try > fail.log 
+             fi
+	 else
+	     ../madevent 2>&1 >> log.txt <input_sg.txt | tee -a log.txt;
+	     status_code=${PIPESTATUS[0]};
+	     if [ -s log.txt ]
+             then
+		 rm log.txt
+                 break
+             else
+                  echo $try > fail.log
+             fi
+	 fi
      done
-     echo "" >> $k; echo "ls status:" >> $k; ls >> $k
-     cat $k >> log.txt
      if [[ $status_code -ne 0 ]]; then 
 	 rm results.dat
 	 echo "ERROR DETECTED"
@@ -69,5 +85,16 @@
      if [[ -e ftn26 ]]; then
          cp ftn26 ftn25
      fi
+
+     if [ "$keeplog" = true ] ; then
+	 echo "" >> $k; echo "ls status:" >> $k; ls >> $k	 
+     else
+	 rm ftn26 > /dev/null
+     fi
+
+
+
+
+
      cd ../
 

=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2016-02-11 10:49:17 +0000
+++ Template/LO/SubProcesses/reweight.f	2016-08-19 23:07:19 +0000
@@ -485,7 +485,6 @@
 C   ARGUMENTS 
 C   
       DOUBLE PRECISION P(0:3,NEXTERNAL)
-
 C   global variables
 C     Present process number
       INTEGER IMIRROR,IPROC
@@ -557,8 +556,8 @@
 c     First time, cluster according to this config and store jets
 c     (following times, only accept configurations if the same partons
 c      are flagged as jets)
+      chclusold=chcluster
       if(njetstore(iconfig).eq.-1)then
-         chclusold=chcluster
          chcluster=.true.
       endif
  100  clustered = cluster(p(0,1))
@@ -844,7 +843,7 @@
      $           write(*,*) 'Bad clustering, jets fail. Reclustering ',
      $           iconfig
             chcluster=.true.
-c            goto 100 ! not
+            goto 100
          endif
       endif
       
@@ -1214,11 +1213,14 @@
       endif
 
 
-      if(.not.setclscales(p)) then ! recluster to have the correct iqjets
-        write(*,*) "Fail to cluster the events from the rewgt function"
-        stop
+      if(.not.setclscales(p)) then ! assign the correct id information.
+         write(*,*) "Fail to cluster the events from the rewgt function"
+         stop 1
+c        rewgt = 0d0
+        return
       endif
 
+
 c     Store pdf information for systematics studies (initial)
 c     need to be done after      setclscales since that one clean the syscalc value
       if(use_syst)then

=== modified file 'Template/LO/SubProcesses/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2016-05-11 10:34:14 +0000
+++ Template/LO/SubProcesses/setcuts.f	2016-08-19 23:07:19 +0000
@@ -81,6 +81,11 @@
       logical  do_cuts(nexternal)
       COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
      . IS_A_ONIUM, do_cuts 
+
+c     Store which external particles undergo the ktdurham and ptlund cuts.
+      LOGICAL  is_pdg_for_merging_cut(NEXTERNAL)
+      COMMON /TO_MERGE_CUTS/is_pdg_for_merging_cut, from_decay
+
 c
 c
 c     reading parameters
@@ -107,6 +112,11 @@
      $     ,ptll_max_save,etmin_save,etmax_save,emin_save,emax_save
      $     ,etamin_save,etamax_save
 c
+c     count the number of j/bjet/photon/lepton
+c
+      integer nb_j, nb_b, nb_a, nb_l, nb_nocut
+      double precision smin_p, smin_m ! local variable to compute smin
+c
 c     setup masses for the final-state particles
 c
       include 'pmass.inc'
@@ -133,6 +143,8 @@
             fixed_fac_scale=.true.
             use_syst=.false.
          endif
+
+
 c
 c     set ptj and s_min if xqcut and ktscheme = 1, to improve
 c     integration speed, and set drjj and drjl to 0.
@@ -220,62 +232,24 @@
          if (abs(idup(i,1,iproc)).eq.14) is_a_nu(i)=.true.  ! no cuts on vm vm~
          if (abs(idup(i,1,iproc)).eq.16) is_a_nu(i)=.true.  ! no cuts on vt vt~
          if (pmass(i).gt.10d0)     is_heavy(i)=.true. ! heavy fs particle
-c-onium
-         if (idup(i,1,iproc).eq.441)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.10441)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.100441)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.443)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.10443)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.20443)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.100443)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.30443)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.9000443)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.9010443)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.9020443)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.445)   is_a_onium(i)=.true. ! charmonium
-         if (idup(i,1,iproc).eq.9000445)   is_a_onium(i)=.true. ! charmonium
-
-         if (idup(i,1,iproc).eq.551)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.10551)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.100551)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.110551)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.200551)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.210551)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.10553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.20553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.30553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.100553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.110553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.120553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.130553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.200553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.210553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.220553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.300553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.9000553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.9010553)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.555)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.10555)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.20555)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.100555)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.110555)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.200555)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.557)   is_a_onium(i)=.true. ! bottomonium
-         if (idup(i,1,iproc).eq.100557)   is_a_onium(i)=.true. ! bottomonium
-
-         if (idup(i,1,iproc).eq.541)   is_a_onium(i)=.true. ! Bc
-         if (idup(i,1,iproc).eq.10541)   is_a_onium(i)=.true. ! Bc
-         if (idup(i,1,iproc).eq.543)   is_a_onium(i)=.true. ! Bc
-         if (idup(i,1,iproc).eq.10543)   is_a_onium(i)=.true. ! Bc
-         if (idup(i,1,iproc).eq.20543)   is_a_onium(i)=.true. ! Bc
-         if (idup(i,1,iproc).eq.545)   is_a_onium(i)=.true. ! Bc
+
+c        Remember mergeable particles
+         do j=1,pdgs_for_merging_cut(0)
+           if (    pdgs_for_merging_cut(j) .ne. 0
+     $       .and. abs(idup(i,1,iproc)) .eq.pdgs_for_merging_cut(j) ) then
+             is_pdg_for_merging_cut(i)=.true.
+             exit
+            endif
+         enddo
+
       enddo
 
+
+
+
 c
 c     et and eta cuts
 c
-         Smin = 0d0
       do i=nincoming+1,nexternal
          etmin(i)  = 0d0
          etmax(i)  = -1
@@ -291,7 +265,6 @@
 c        JET
             if(is_a_j(i))then
                  etmin(i)=ptj
-                 SMIN = SMIN + max(ptj,ej)
                  etmax(i)=ptjmax
                  emin(i)=ej
                  emax(i)=ejmax
@@ -301,7 +274,6 @@
 c        LEPTON
             if(is_a_l(i))then
                  etmin(i)=ptl
-                 SMIN = SMIN + max(ptl,el)
                  etmax(i)=ptlmax
                  emin(i)=el
                  emax(i)=elmax
@@ -311,7 +283,6 @@
 c        BJET
             if(is_a_b(i))then
                  etmin(i)=ptb
-                 SMIN = SMIN + max(ptb,eb)
                  etmax(i)=ptbmax
                  emin(i)=eb
                  emax(i)=ebmax
@@ -321,7 +292,6 @@
 c        PHOTON
             if(is_a_a(i))then
                  etmin(i) = max(pta, ptgmin)
-                 SMIN = SMIN + max(etmin(i),ea)
                  etmax(i)=ptamax
                  emin(i)=ea
                  emax(i)=eamax
@@ -329,14 +299,16 @@
                  etamin(i)=etaamin
             endif
 c        QUARKONIUM
-            if(is_a_onium(i))then
-                 etmin(i)=ptonium
-                 SMIN = SMIN + ptonium
-                 etamax(i)=etaonium
-            endif
+c            if(is_a_onium(i))then
+c                 etmin(i)=ptonium
+c                 SMIN = SMIN + max(0d0,ptonium**2
+c                 etamax(i)=etaonium
+c            endif
          endif
       enddo
-      SMIN = SMIN **2
+
+
+
 c
 c     delta r cut
 c
@@ -473,35 +445,166 @@
 
 
 c
-c     Compute Smin (for efficiency
+c     Compute Smin (for efficiency purpose)
 c
-      do i=nincoming+1,nexternal-1
-      do j=nincoming+1,nexternal-1
-         if(j.lt.i)then
-            s_min(i,j) = max(s_min(j,i),s_min(i,j))
-         else
-            smin=0.0d0**2
-
-            if(do_cuts(i).and.do_cuts(j)) then
-               if(is_a_j(i).and.is_a_j(j)) s_min(j,i)=mmjj*dabs(mmjj)
-               if(is_a_a(i).and.is_a_a(j)) s_min(j,i)=mmaa*dabs(mmaa)
-               if( is_a_b(i).and.is_a_b(j) ) s_min(j,i)=mmbb*dabs(mmbb)
-               if((is_a_l(i).and.is_a_l(j)).and.
-     &            (abs(idup(i,1,iproc)).eq.abs(idup(j,1,iproc))).and.
-     &            (idup(i,1,iproc)*idup(j,1,iproc).lt.0))
-     &            s_min(j,i)=mmll*dabs(mmll)  !only on l+l- pairs (same flavour)
-
-               if(is_a_j(i).and.is_a_j(j)) s_max(j,i)=mmjjmax*dabs(mmjjmax)
-               if(is_a_a(i).and.is_a_a(j)) s_max(j,i)=mmaamax*dabs(mmaamax)
-               if( is_a_b(i).and.is_a_b(j) ) s_max(j,i)=mmbbmax*dabs(mmbbmax)
-               if((is_a_l(i).and.is_a_l(j)).and.
-     &            (abs(idup(i,1,iproc)).eq.abs(idup(j,1,iproc))).and.
-     &            (idup(i,1,iproc)*idup(j,1,iproc).lt.0))
-     &            s_max(j,i)=mmllmax*dabs(mmllmax)  !only on l+l- pairs (same flavour)
-            endif
-         endif
-      enddo
-      enddo
+      smin = 0d0
+c     check for the jet
+      smin_m = 0d0
+      smin_p = 0d0
+      nb_j = 0
+      nb_nocut = 0
+      do i=nincoming+1,nexternal
+          if (is_a_j(i)) then
+              nb_j = nb_j + 1
+              smin_m = smin_m - pmass(i)**2
+              if (do_cuts(i))then
+                 if (nb_j.eq.1) then
+                  smin_p = smin_p + max(ej,ptj,xptj,0d0,
+     &                                    max(ptj1min,ptj2min,ptj3min,ptj4min))
+                 elseif(nb_j.eq.2) then
+                  smin_p = max(ht2min, smin_p + max(ej,ptj,0d0,
+     &                                            max(ptj2min,ptj3min,ptj4min)))
+                 elseif(nb_j.eq.3) then
+                  smin_p = max(ht3min, smin_p + max(ej,ptj,0d0,
+     &                                                    max(ptj3min,ptj4min)))
+                 elseif(nb_j.eq.4) then
+                  smin_p = max(ht4min, smin_p + max(ej,ptj,ptj4min,0d0))
+                 else
+                  smin_p = smin_p + max(ej,ptj,0d0)
+                 endif
+              else
+                 nb_nocut = nb_nocut + 1
+                 if (nb_j.eq.1) then
+                  smin_p = smin_p + max(0d0,ptj1min,ptj2min,ptj3min,ptj4min)
+                 elseif(nb_j.eq.2) then
+                  smin_p = max(ht2min, smin_p + max(ptj2min,ptj3min,ptj4min,0d0))
+                 elseif(nb_j.eq.3) then
+                  smin_p = max(ht3min, smin_p + max(0d0,ptj3min,ptj4min))
+                 elseif(nb_j.eq.4) then
+                  smin_p = max(ht4min, smin_p + max(ptj4min,0d0))
+                 endif
+              endif
+          endif
+        enddo
+      if (nb_j.gt.0)then
+         if ((nb_j-nb_nocut).gt.0)then
+            smin_m = smin_m + (nb_j-nb_nocut)*(nb_j-nb_nocut-1)/2d0*mmjj**2
+         endif
+        smin = smin + max(smin_p**2, smin_m, htjmin**2)
+      endif
+c     check for the bjet
+      smin_m = 0d0
+      smin_p = 0d0
+      nb_b = 0
+      nb_nocut = 0
+      do i=nincoming+1,nexternal
+          if (is_a_b(i).and.do_cuts(i)) then
+              nb_b = nb_b + 1
+              smin_m = smin_m - pmass(i)**2
+              if (do_cuts(i)) then
+                 if (nb_b.eq.1) then
+                    smin_p = smin_p + max(eb,ptb,xptb,0d0)
+                 else
+                    smin_p = smin_p + max(eb,ptb,0d0)
+                 endif
+              else
+                 nb_nocut = nb_nocut +1
+                 if (nb_b.eq.1) then
+                    smin_p = smin_p + max(xptb,0d0)
+                 endif
+              endif
+          endif
+      enddo
+      if (nb_b.gt.0)then
+         if ((nb_b-nb_nocut).gt.0) then
+            smin_m = smin_m + (nb_b-nb_nocut)*(nb_b-nb_nocut-1)/2d0*mmbb**2
+         endif
+         smin = smin + max(smin_p**2, smin_m, (ihtmin**2-htjmin**2))
+      endif
+c     check for the photon
+      smin_m = 0d0
+      smin_p = 0d0
+      nb_a = 0
+      nb_nocut = 0
+      do i=nincoming+1,nexternal
+         if (is_a_a(i))then
+            nb_a = nb_a + 1
+            if (do_cuts(i))then
+               if (ptgmin.eq.0d0) then
+                  if (nb_a.eq.1) then
+                     smin_p = smin_p + max(ea,pta,xpta,0d0)
+                  else
+                     smin_p = smin_p + max(ea,pta,0d0)
+                  endif
+               endif
+            else
+               nb_nocut = nb_nocut + 1
+               if(ptgmin.eq.0d0)then
+                  if (nb_a.eq.1) then
+                     smin_p = smin_p + max(xpta,0d0)
+                  endif
+               endif
+            endif
+         endif
+      enddo
+      if (nb_a.gt.0) then
+         if ((nb_a-nb_nocut).gt.0)then
+            smin_m = (nb_a-nb_nocut)*(nb_a-nb_nocut-1)/2d0*mmaa**2
+         endif
+        smin = smin + max(smin_p**2, smin_m)
+      endif
+c     check for lepton
+      smin_m = 0d0
+      smin_p = 0d0
+      nb_l = 0
+      nb_nocut = 0
+      do i=nincoming+1,nexternal-1
+         if (is_a_l(i))then
+            nb_l = nb_l + 1
+            smin_m = smin_m - pmass(i)**2
+            if (do_cuts(i))then 
+               if (nb_l.eq.1) then
+                  smin_p = smin_p + max(el,ptl,xptl,0d0,
+     &                 max(ptl1min,ptl2min,ptl3min,ptl4min))
+               elseif(nb_l.eq.2) then
+                  smin_p = smin_p + max(el,ptl,0d0,max(ptl2min,ptl3min,ptl4min))
+               elseif(nb_l.eq.3) then
+                  smin_p = smin_p + max(el,ptl,0d0, max(ptl3min,ptl4min))
+               elseif(nb_l.eq.4) then
+                  smin_p = smin_p + max(el,ptl,ptl4min,0d0)
+               else
+                  smin_p = smin_p + max(el,ptl,0d0)
+               endif
+            else
+               nb_nocut = nb_nocut + 1
+               if (nb_l.eq.1) then
+                  smin_p = smin_p + max(xptl,0d0,
+     &                 max(ptl1min,ptl2min,ptl3min,ptl4min))
+               elseif(nb_l.eq.2) then
+                  smin_p = smin_p + max(0d0,ptl2min,ptl3min,ptl4min)
+               elseif(nb_l.eq.3) then
+                  smin_p = smin_p + max(0d0,ptl3min,ptl4min)
+               elseif(nb_l.eq.4) then
+                  smin_p = smin_p + max(ptl4min,0d0)
+               endif
+            endif
+         endif
+      enddo
+      if (nb_l.gt.0)then
+         if ((nb_l-nb_nocut).gt.0)then
+            smin_m = smin_m + (nb_l-nb_nocut)*((nb_l-nb_nocut)-1)/2d0*mmll**2
+         endif
+        smin = smin + max(smin_p**2, smin_m, mmnl**2, ptllmin**2, misset**2)
+      endif
+c     ensure symmetry of s_min(i,j)
+      do i=nincoming+1,nexternal-1
+        do j=nincoming+1,nexternal-1
+           if(j.lt.i)then
+               s_min(i,j) = max(s_min(j,i),s_min(i,j))
+            endif
+        enddo
+      enddo
+      write(*,*) "Define smin to", smin
 
 c Check that results are consistent among all the grouped subprocesses
 

=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2015-02-18 01:33:22 +0000
+++ Template/LO/SubProcesses/unwgt.f	2016-08-19 23:07:19 +0000
@@ -138,6 +138,38 @@
       nw = 0
       maxwgt = 0d0
       end
+C**************************************************************************
+C      HELPING ROUTINE FOR PERFORMING THE Z BOOST OF THE EVENTS
+C**************************************************************************
+      double precision function get_betaz(pin,pout)
+C     compute the boost for the requested transformation
+      implicit none
+      double precision pin(0:3), pout(0:3)
+      double precision denom
+
+      denom = pin(0)*pout(0) - pin(3)*pout(3)
+      if (denom.ne.0d0) then
+         get_betaz = (pin(3) * pout(0) - pout(3) * pin(0)) / denom
+      else
+         get_betaz = 0d0
+      endif
+      return
+      end
+
+      subroutine zboost_with_beta(pin, beta, pout)
+c     apply the boost
+      implicit none
+      double precision pin(0:3), pout(0:3)
+      double precision beta, gamma
+
+      gamma = 1d0/DSQRT(1-beta**2)
+      pout(0) = gamma * pin(0) - gamma * beta * pin(3)
+      pout(1) = pin(1)
+      pout(2) = pin(2)
+      pout(3) = - gamma * beta * pin(0) + gamma * pin(3)
+      return
+      end
+
 
       SUBROUTINE unwgt(px,wgt,numproc)
 C**************************************************************************
@@ -436,7 +468,9 @@
       integer ip, np, ic, nc, jpart(7,-nexternal+3:2*nexternal-3)
       integer ida(2),ito(-nexternal+3:nexternal),ns,nres,ires,icloop
       integer iseed
-      double precision pboost(0:3),pb(0:4,-nexternal+3:2*nexternal-3),eta
+      double precision pboost(0:3),pb(0:4,-nexternal+3:2*nexternal-3)
+      double precision beta, get_betaz
+      double precision ebi(0:3), ebo(0:3)
       double precision ptcltmp(nexternal), pdum(0:3)
 
       integer idup(nexternal,maxproc,maxsproc)
@@ -564,10 +598,21 @@
       if (nincoming.eq.2) then
          if (xbk(1) .gt. 0d0 .and. xbk(1) .lt. 1d0 .and.
      $       xbk(2) .gt. 0d0 .and. xbk(2) .lt. 1d0) then
-            eta = sqrt(xbk(1)*ebeam(1)/
-     $                 (xbk(2)*ebeam(2)))
-            pboost(0)=p(0,1)*(eta + 1d0/eta)
-            pboost(3)=p(0,1)*(eta - 1d0/eta)
+            pboost(0) = -1d0 ! special flag
+            ! construct the beam momenta in each frame and compute the related (z)boost
+            ebi(0) = p(0,1)/xbk(1)
+            ebi(1) = 0
+            ebi(2) = 0
+            ebi(3) = DSQRT(ebi(0)**2-m1**2)
+            ebo(0) = ebeam(1)
+            ebo(1) = 0
+            ebo(2) = 0
+            ebo(3) = DSQRT(ebo(0)**2-m1**2)
+            beta = get_betaz(ebi, ebo)
+c           eta = sqrt(xbk(1)*ebeam(1)/
+c     $                 (xbk(2)*ebeam(2)))
+c            pboost(0)=p(0,1)*(eta + 1d0/eta)
+c            pboost(3)=p(0,1)*(eta - 1d0/eta)
          else if (xbk(1) .gt. 0d0 .and. xbk(1) .lt. 1d0 .and.
      $       xbk(2) .eq. 1d0) then
             pboost(0)=xbk(1)*ebeam(1)+ebeam(2)
@@ -587,12 +632,18 @@
      $           xbk(1),xbk(2)
          endif
       endif
-      do j=1,nexternal
-         call boostx(p(0,j),pboost,pb(0,isym(j,jsym)))
-c      Add mass information in pb(4)
-         pb(4,isym(j,jsym))=pmass(j)
-      enddo
-
+      if (pboost(0).ne.-1d0) then
+         do j=1,nexternal
+            call boostx(p(0,j),pboost,pb(0,isym(j,jsym)))
+            ! Add mass information in pb(4)
+            pb(4,isym(j,jsym))=pmass(j)
+         enddo
+      else
+         do j=1,nexternal
+            call zboost_with_beta(p(0,j),beta,pb(0,isym(j,jsym)))
+            pb(4,isym(j,jsym))=pmass(j)
+         enddo
+      endif
 c
 c     Add info on resonant mothers
 c

=== modified file 'Template/LO/bin/internal/run_delphes3'
--- Template/LO/bin/internal/run_delphes3	2013-02-25 00:52:53 +0000
+++ Template/LO/bin/internal/run_delphes3	2016-08-19 23:07:19 +0000
@@ -9,6 +9,7 @@
 run=$2
 tag=$3
 cross=$4 
+file=$5
 
 main=`pwd`
 
@@ -22,17 +23,34 @@
   exit
 fi
 
+if [ ${file: -3} == ".gz" ]
+then
+  if [ ${file: -7} == ".hep.gz" ]
+  then
+      gunzip --stdout $file | $delphesdir/DelphesSTDHEP ../Cards/delphes_card.dat  ${run}/${tag}_delphes_events.root      
+  else
+      gunzip --stdout $file | $delphesdir/DelphesHepMC ../Cards/delphes_card.dat  ${run}/${tag}_delphes_events.root      
+  fi
+else
+  if [ ${file: -4} == ".hep" ]
+  then
+      $delphesdir/DelphesSTDHEP ../Cards/delphes_card.dat  ${run}/${tag}_delphes_events.root $file
+  else
+      $delphesdir/DelphesHepMC ../Cards/delphes_card.dat  ${run}/${tag}_delphes_events.root $file      
+  fi
+fi
+
 echo $$ >> ../myprocid
 
-# Set delphes path in delphes_card.dat
-
-$delphesdir/DelphesSTDHEP ../Cards/delphes_card.dat delphes.root pythia_events.hep
-$delphesdir/root2lhco delphes.root delphes_events.lhco
-
-if [ -e delphes_events.lhco ]; then
-# write the delphes banner
-    sed -e "s/^/#/g" ${run}/${run}_${tag}_banner.txt > ${run}/${tag}_delphes_events.lhco
-    echo "##  Integrated weight (pb)  : ${cross}" >> ${run}/${tag}_delphes_events.lhco
-    cat delphes_events.lhco >> ${run}/${tag}_delphes_events.lhco
-    rm -f delphes_events.lhco
-fi
+# Uncomment the following to have the LHCO file:
+#
+#
+#$delphesdir/root2lhco ${run}/${tag}_delphes.root delphes_events.lhco
+#
+#if [ -e delphes_events.lhco ]; then
+## write the delphes banner
+#    sed -e "s/^/#/g" ${run}/${run}_${tag}_banner.txt > ${run}/${tag}_delphes_events.lhco
+#    echo "##  Integrated weight (pb)  : ${cross}" >> ${run}/${tag}_delphes_events.lhco
+#    cat delphes_events.lhco >> ${run}/${tag}_delphes_events.lhco
+#    rm -f delphes_events.lhco
+#fi

=== modified file 'Template/LO/bin/internal/store4grid'
--- Template/LO/bin/internal/store4grid	2012-09-05 07:27:45 +0000
+++ Template/LO/bin/internal/store4grid	2016-08-19 23:07:19 +0000
@@ -32,14 +32,11 @@
 if [[ -e default_results.dat ]]; then
     rm -f default_results.dat
 fi    
-cp -f $1_results.dat default_results.dat
 for i in `cat subproc.mg` ; do
     cd $i
-    echo $i
     if [[ -e default_results.dat ]]; then
 	rm -f default_results.dat
     fi    
-    cp -f $1_results.dat default_results.dat
     for k in G* ; do
 	if [[ ! -d $k ]]; then
 	    continue

=== modified file 'Template/LO/bin/madevent'
--- Template/LO/bin/madevent	2016-02-23 11:25:14 +0000
+++ Template/LO/bin/madevent	2016-08-19 23:07:19 +0000
@@ -205,7 +205,6 @@
     if cmd_line.history[-1] not in ['EOF','quit','exit']:
         cmd_line.results.store_result()
 except:
-    print error
     pass
 
 

=== modified file 'Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_h.f'
--- Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_h.f	2016-02-23 10:48:46 +0000
+++ Template/NLO/MCatNLO/PY8Analyzer/py8an_HwU_pp_h.f	2016-08-19 23:07:19 +0000
@@ -159,7 +159,7 @@
       common/cww/ww
 c
       if(nnn.eq.0)ww(1)=1d0
-      do i=,nnn
+      do i=1,nnn
          ww(i)=xww(i)
       enddo
 c

=== modified file 'Template/NLO/Source/make_opts.inc'
--- Template/NLO/Source/make_opts.inc	2016-03-11 11:20:49 +0000
+++ Template/NLO/Source/make_opts.inc	2016-08-19 23:07:19 +0000
@@ -55,7 +55,7 @@
 
 UNAME := $(shell uname -s)
 ifeq ($(origin LDFLAGS), undefined)
-     LDFLAGS=$(STDLIB)
+     LDFLAGS=$(STDLIB) $(MACFLAG)
 endif
 # Options: dynamic, lhapdf
 
@@ -89,10 +89,10 @@
 
 # Option lhapdf
 
-ifdef lhapdf
+ifneq ($(lhapdf),)
   CXXFLAGS += $(shell $(lhapdf) --cppflags)
   alfas_functions=alfas_functions_lhapdf
-  llhapdf+=-lLHAPDF 
+  llhapdf = $(shell $(lhapdf) --libs)
   reweight_xsec_events_pdf_dummy=
 else
   alfas_functions=alfas_functions
@@ -101,7 +101,7 @@
 endif
 
 # Option APPLGrid
-ifdef applgrid
+ifneq ($(applgrid),)
   APPLLIBS=$(shell applgrid-config --ldcflags) $(shell amcfast-config --ldflags)
   applgrid_interface=appl_interface.o
 else

=== modified file 'Template/NLO/Source/setrun.f'
--- Template/NLO/Source/setrun.f	2016-03-21 08:45:42 +0000
+++ Template/NLO/Source/setrun.f	2016-08-19 23:07:19 +0000
@@ -204,6 +204,10 @@
 c fill the nmemPDF(i) array with the number of PDF error set. This we
 c get from LHAPDF.
          call numberPDFm(1,nmemPDF(1))
+         if (nmemPDF(1).eq.1) then
+            nmemPDF(1)=0
+            lpdfvar(1)=0
+         endif
       else
          nmemPDF(1)=0
       endif

=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
--- Template/NLO/SubProcesses/BinothLHA.f	2016-03-02 09:31:23 +0000
+++ Template/NLO/SubProcesses/BinothLHA.f	2016-08-19 23:07:19 +0000
@@ -85,6 +85,8 @@
 c itself again to perform stability check to make sure no unstable EPS
 c splips unnoticed.
          CALL FORCE_STABILITY_CHECK(.TRUE.)
+         CALL COLLIER_COMPUTE_UV_POLES(.FALSE.)
+         CALL COLLIER_COMPUTE_IR_POLES(.FALSE.)
          firsttime_run = .false.
       endif
       if (firsttime) then
@@ -102,6 +104,7 @@
 c Just set the accuracy found to a positive value as it is not specified
 c once the initial pole check is performed.
          if (mc_hel.eq.0) then
+
             call sloopmatrix_thres(p,virt_wgts,tolerance,accuracies
      $           ,ret_code)
             prec_found = accuracies(0)            
@@ -164,18 +167,19 @@
 c MadLoop initialization PS points.
       cpol=.false.
       if ((firsttime .or. mc_hel.eq.0) .and. mod(ret_code,100)/10.ne.3
-     $     .and. mod(ret_code,100)/10.ne.4) then 
+     $     .and. mod(ret_code,100)/10.ne.4) then
          call getpoles(p,QES2,madfks_double,madfks_single,fksprefact)
          avgPoleRes(1)=(single+madfks_single)/2.0d0
          avgPoleRes(2)=(double+madfks_double)/2.0d0
          PoleDiff(1)=dabs(single - madfks_single)
          PoleDiff(2)=dabs(double - madfks_double)
          if ((dabs(avgPoleRes(1))+dabs(avgPoleRes(2))).ne.0d0) then
-            cpol = .not. (((PoleDiff(1)+PoleDiff(2))/
+            cpol = .not. ((((PoleDiff(1)+PoleDiff(2))/
      $           (dabs(avgPoleRes(1))+dabs(avgPoleRes(2)))) .lt.
-     $           tolerance*10d0)
+     $           tolerance*10d0).or.(mod(ret_code,10).eq.7))
          else
-            cpol = .not.(PoleDiff(1)+PoleDiff(2).lt.tolerance*10d0)
+            cpol = .not.((PoleDiff(1)+PoleDiff(2).lt.tolerance*10d0).or.
+     $                   (mod(ret_code,10).eq.7))
          endif
          if (tolerance.lt.0.0d0) then
             cpol = .false.

=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2016-03-02 09:31:23 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2016-08-19 23:07:19 +0000
@@ -116,6 +116,8 @@
 c Make sure that stability checks are always used by MadLoop, even for
 c initialization
       CALL FORCE_STABILITY_CHECK(.TRUE.)
+      CALL COLLIER_COMPUTE_UV_POLES(.TRUE.)
+      CALL COLLIER_COMPUTE_IR_POLES(.TRUE.)
 
 200   continue
           calculatedborn = .false.

=== modified file 'Template/NLO/SubProcesses/madfks_plot.f'
--- Template/NLO/SubProcesses/madfks_plot.f	2016-02-25 08:21:12 +0000
+++ Template/NLO/SubProcesses/madfks_plot.f	2016-08-19 23:07:19 +0000
@@ -66,6 +66,10 @@
                call initpdfsetbynamem(nn,lhaPDFsetname(nn))
                if (lpdfvar(nn)) then
                   call numberPDFm(nn,nmemPDF(nn))
+                  if (nmemPDF(nn).eq.1) then
+                     nmemPDF(nn)=0
+                     lpdfvar(nn)=.false.
+                  endif
                else
                   nmemPDF(nn)=0
                endif

=== modified file 'Template/NLO/SubProcesses/makefile_fks_dir'
--- Template/NLO/SubProcesses/makefile_fks_dir	2016-03-21 08:45:42 +0000
+++ Template/NLO/SubProcesses/makefile_fks_dir	2016-08-19 23:07:19 +0000
@@ -9,11 +9,11 @@
 ifdef fastjet_config
     FASTJET_CONFIG = $(fastjet_config)
     CXXFLAGS += $(shell $(FASTJET_CONFIG) --cxxflags)
-    FJLIBS += $(shell $(FASTJET_CONFIG) --libs --plugins ) -lstdc++
+    FJLIBS += $(shell $(FASTJET_CONFIG) --libs --plugins ) $(STDLIB)
     fastjetfortran_madfks=fastjetfortran_madfks_full.o fastjet_wrapper.o
 else
     fastjetfortran_madfks=fastjetfortran_madfks_core.o fjcore.o fastjet_wrapper.o
-    FJLIBS += -lstdc++
+    FJLIBS += $(STDLIB)
 endif
 
 NLOLIBS=-L$(LIBDIR) $(libmadloop) $(libcuttools) $(libOLP)
@@ -111,30 +111,30 @@
 	strip test_MC
 
 check_poles: $(POLES) $(libmadloop)
-	$(FC) $(LDFLAGS) -o check_poles $(POLES) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS)
+	$(FC) -o check_poles $(POLES) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(LDFLAGS)
 	strip check_poles
 
 madevent_vegas: $(VEGAS) $(libmadloop) makefile $(LIBS)
-	$(FC) $(LDFLAGS) -o madevent_vegas $(VEGAS) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS)
+	$(FC) -o madevent_vegas $(VEGAS) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(LDFLAGS)
 	rm handling_lhe_events.o
 	strip madevent_vegas
 
 madevent_mintMC: $(MINTMC) $(libmadloop) makefile $(LIBS)
-	$(FC) $(LDFLAGS) -o madevent_mintMC $(MINTMC) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS)
+	$(FC) -o madevent_mintMC $(MINTMC) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(LDFLAGS)
 	rm handling_lhe_events.o
 	strip madevent_mintMC
 
 madevent_mintFO: $(MINTFO) $(libmadloop) makefile $(LIBS)
-	$(FC) $(LDFLAGS) -o madevent_mintFO $(MINTFO) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS)
+	$(FC) -o madevent_mintFO $(MINTFO) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(LDFLAGS)
 	rm handling_lhe_events.o
 	strip madevent_mintFO
 
 madevent_reweight: $(BORN_REWEIGHT) $(libmadloop) makefile $(LIBS)
-	$(FC) $(LDFLAGS) -o madevent_reweight $(BORN_REWEIGHT) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS)
+	$(FC) -o madevent_reweight $(BORN_REWEIGHT) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(LDFLAGS)
 	rm handling_lhe_events.o
 
 reweight_xsec_events: $(RWGFILES) makefile $(LIBS)
-	$(FC) $(LDFLAGS) -o reweight_xsec_events $(RWGFILES) $(LINKLIBS) $(FJLIBS)
+	$(FC) -o reweight_xsec_events $(RWGFILES) $(LINKLIBS) $(FJLIBS) $(LDFLAGS)
 	rm handling_lhe_events.o
 	strip reweight_xsec_events
 

=== modified file 'Template/NLO/SubProcesses/makefile_loop.inc'
--- Template/NLO/SubProcesses/makefile_loop.inc	2016-03-02 03:03:20 +0000
+++ Template/NLO/SubProcesses/makefile_loop.inc	2016-08-19 23:07:19 +0000
@@ -24,6 +24,7 @@
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_coef_construction_*.f)) \
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard TIR_interface.f)) \
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard GOLEM_interface.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard COLLIER_interface.f)) \
 		 $(patsubst %(dotf)s,%(doto)s,$(wildcard compute_color_flows.f))
 
 # This is the core of madloop computationally wise, so make sure to turn optimizations on and bound checks off.

=== modified file 'Template/NLO/SubProcesses/reweight0.inc'
--- Template/NLO/SubProcesses/reweight0.inc	2016-02-23 10:48:46 +0000
+++ Template/NLO/SubProcesses/reweight0.inc	2016-08-19 23:07:19 +0000
@@ -71,7 +71,7 @@
 c New format to allow for multiple PDF sets and scales (both functional form
 c and normal)
       integer    maxPDFsets,   maxdynscales
-      parameter (maxPDFsets=10,maxdynscales=10)
+      parameter (maxPDFsets=25,maxdynscales=10)
       integer lhaPDFid(0:maxPDFsets),nmemPDF(maxPDFsets)
      $     ,dyn_scale(0:maxdynscales)
       logical lscalevar(maxdynscales),lpdfvar(maxPDFsets)

=== modified file 'Template/NLO/SubProcesses/reweight_xsec_events.f'
--- Template/NLO/SubProcesses/reweight_xsec_events.f	2016-03-30 23:53:52 +0000
+++ Template/NLO/SubProcesses/reweight_xsec_events.f	2016-08-19 23:07:19 +0000
@@ -87,6 +87,10 @@
                call initpdfsetbynamem(nn,lhaPDFsetname(nn))
                if (lpdfvar(nn)) then
                   call numberPDFm(nn,nmemPDF(nn))
+                  if (nmemPDF(nn).eq.1) then
+                     nmemPDF(nn)=0
+                     lpdfvar(nn)=.false.
+                  endif
                else
                   nmemPDF(nn)=0
                endif

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_test_MC.f'
--- Template/NLO/SubProcesses/symmetry_fks_test_MC.f	2016-03-16 14:21:42 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_test_MC.f	2016-08-19 23:07:19 +0000
@@ -13,6 +13,7 @@
       include 'nFKSconfigs.inc'
       include 'fks_info.inc'
       include 'run.inc'
+      include 'cuts.inc'
       
       double precision ZERO,one
       parameter       (ZERO = 0d0)
@@ -145,6 +146,9 @@
       integer fks_conf_number,fks_loop_min,fks_loop_max,fks_loop
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
+      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
       
       character*10 MonteCarlo
       common/cMonteCarloType/MonteCarlo
@@ -222,8 +226,14 @@
 c Note that tests are always performed at fixed energy with Bjorken x=1.
       totmass = 0.0d0
       include 'pmass.inc' ! make sure to set the masses after the model has been included
-      do i=1,nexternal
-        totmass = totmass + pmass(i)
+      do i=nincoming+1,nexternal
+         if (is_a_j(i) .and. i.ne.nexternal) then
+            totmass = totmass + max(ptj,pmass(i))
+         elseif ((is_a_lp(i).or.is_a_lm(i)) .and. i.ne.nexternal) then
+            totmass = totmass + max(mll/2d0,mll_sf/2d0,ptl,pmass(i))
+         else
+            totmass = totmass + pmass(i)
+         endif
       enddo
       if (lpp(1).ne.0) ebeam(1)=max(ebeam(1)/20d0,totmass)
       if (lpp(2).ne.0) ebeam(2)=max(ebeam(2)/20d0,totmass)
@@ -394,7 +404,6 @@
          wgt=1d0
          call generate_momenta(ndim,iconfig,wgt,x,p)
          calculatedBorn=.false.
-
          call set_cms_stuff(0)
          calculatedBorn=.false.
 

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_test_ME.f'
--- Template/NLO/SubProcesses/symmetry_fks_test_ME.f	2016-03-16 14:21:42 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_test_ME.f	2016-08-19 23:07:19 +0000
@@ -13,6 +13,7 @@
       include 'nFKSconfigs.inc'
       include 'fks_info.inc'
       include 'run.inc'
+      include 'cuts.inc'
       
       double precision ZERO,one
       parameter       (ZERO = 0d0)
@@ -134,6 +135,9 @@
       integer fks_conf_number,fks_loop_min,fks_loop_max,fks_loop
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
+      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
       
 c      integer icomp
 c
@@ -170,8 +174,14 @@
 c Note that tests are always performed at fixed energy with Bjorken x=1.
       totmass = 0.0d0
       include 'pmass.inc' ! make sure to set the masses after the model has been included
-      do i=1,nexternal
-        totmass = totmass + pmass(i)
+      do i=nincoming+1,nexternal
+         if (is_a_j(i) .and. i.ne.nexternal) then
+            totmass = totmass + max(ptj,pmass(i))
+         elseif ((is_a_lp(i).or.is_a_lm(i)) .and. i.ne.nexternal) then
+            totmass = totmass + max(mll/2d0,mll_sf/2d0,ptl,pmass(i))
+         else
+            totmass = totmass + pmass(i)
+         endif
       enddo
       if (lpp(1).ne.0) ebeam(1)=max(ebeam(1)/20d0,totmass)
       if (lpp(2).ne.0) ebeam(2)=max(ebeam(2)/20d0,totmass)

=== modified file 'Template/NLO/SubProcesses/write_event.f'
--- Template/NLO/SubProcesses/write_event.f	2016-04-16 20:30:38 +0000
+++ Template/NLO/SubProcesses/write_event.f	2016-08-19 23:07:19 +0000
@@ -229,12 +229,6 @@
       logical firsttime
       data firsttime/.true./
 c
-      if (firsttime) then
-         call write_header_init
-         firsttime=.false.
-      endif
-      ievent=66
-
       if (ickkw.eq.4) then
          scale = sqrt(muF12_current)
       elseif (ickkw.eq.-1) then
@@ -245,7 +239,16 @@
 
       aqcd=g**2/(4d0*pi)
       aqed=gal(1)**2/(4d0*pi)
-
+c
+c 'write_header_init' should be called after 'aqcd' has been set,
+c because it includes a call to 'setrun', which resets the value of
+c alpha_s to the one in the param_card.dat (without any running).
+      if (firsttime) then
+         call write_header_init
+         firsttime=.false.
+      endif
+      ievent=66
+c
       if(AddInfoLHE)then
         if(.not.doreweight)then
            write(buff,201)'#aMCatNLO',iSorH_lhe,ifks_lhe(nFKSprocess)

=== modified file 'Template/NLO/Utilities/VetoPrefactors/resum_reweighter.py'
--- Template/NLO/Utilities/VetoPrefactors/resum_reweighter.py	2015-03-06 09:47:53 +0000
+++ Template/NLO/Utilities/VetoPrefactors/resum_reweighter.py	2016-08-19 23:07:19 +0000
@@ -115,7 +115,7 @@
 outputEvtFile.write(evtFile.banner)
 # List the channel present in the process output
 channel_list = [ os.path.basename(chan) for chan in \
-                 glob.glob(pjoin(proc_path, 'SubProcesses','P*')) ]
+                 misc.glob('P*', pjoin(proc_path, 'SubProcesses')) ]
 
 # Now scan over events
 # For each encountered channel, store the corresponding process runner

=== modified file 'Template/NLO/Utilities/VetoPrefactors/virt_reweighter.py'
--- Template/NLO/Utilities/VetoPrefactors/virt_reweighter.py	2015-03-06 09:47:53 +0000
+++ Template/NLO/Utilities/VetoPrefactors/virt_reweighter.py	2016-08-19 23:07:19 +0000
@@ -183,7 +183,7 @@
 
 # List the channel present in the process output
 channel_list = [ os.path.basename(chan) for chan in \
-                 glob.glob(pjoin(proc_path, 'SubProcesses','P*')) ]
+                 misc.glob('P*', pjoin(proc_path, 'SubProcesses'))) ]
 
 # Now scan over events
 # For each encountered channel, store the corresponding process runner
@@ -269,8 +269,8 @@
 
         if ChannelNameToSetup!="":
             virtFolders = [ os.path.basename(virt) for virt in \
-                glob.glob(pjoin(proc_path, 'SubProcesses',
-                chanFolder,'V*%s*'%ChannelNameToSetup)) ]
+                misc.glob('V*%s*'%ChannelNameToSetup, pjoin(proc_path, 'SubProcesses',
+                chanFolder)) ]
             if len(virtFolders)!=1:
                 if len(virtFolders)==0:
                     logging.error("No virtual folder V*%s* found in %s."\

=== added file 'Template/NLO/Utilities/thrbdec2_usage.f'
--- Template/NLO/Utilities/thrbdec2_usage.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/Utilities/thrbdec2_usage.f	2016-08-19 23:07:19 +0000
@@ -0,0 +1,430 @@
+c From madfks_plot.ttlep5.f (used in 1407.2763)
+      subroutine outfun(pp,ybst_til_tolab,www,itype)
+.....
+c masses
+      double precision pmass(nexternal)
+      common/to_mass/pmass
+      real*8 rndec(10)
+      common/crndec/rndec
+      real*8 wgt1,wgt2,wwwd
+      real*8 ptop(5),ptbr(5),xpbq(5),xplp(5),xpnu(5),xpbb(5),
+     # xplm(5),xpnb(5)
+      integer maxnum,imode,iunwgt,itqrkseed,itbarseed,ntr
+      parameter (maxnum=10000000)
+      parameter (imode=1)
+c iunwgt=0 for weighted decays, iunwgt=1 for unweighted decays
+      parameter (iunwgt=1)
+      logical xext
+      parameter (xext=.true.)
+      double precision wmass
+
+.....
+
+c Stable tops
+c Assume t is particle #3, tbar is particle #4; light parton is #5
+      do i=1,4
+        ptop(i)=pplab(mod(i,4),3)
+        ptbr(i)=pplab(mod(i,4),4)
+      enddo
+      ptop(5)=pmass(3)
+      ptbr(5)=pmass(4)
+c
+      wmass=MW
+      itqrkseed=max(1,int(maxnum*rndec(1)))
+      itbarseed=max(1,int(maxnum*rndec(2)))
+      call thrbdec2(xext,imode,iunwgt,itqrkseed,ntr,
+     #              rndec(1),rndec(2),rndec(3),rndec(4),rndec(5),
+     #              wmass,ptop,xpbq,xplp,xpnu,wgt1)
+      call thrbdec2(xext,imode,iunwgt,itbarseed,ntr,
+     #              rndec(6),rndec(7),rndec(8),rndec(9),rndec(10),
+     #              wmass,ptbr,xpbb,xplm,xpnb,wgt2)
+      wwwd=wgt1*wgt2
+
+
+c
+c
+c
+c Package for top decay; the b and the leptons are taken to be massless
+c This is an upgrade of the package written for Wt production, in that it
+c gives one the possibility of using an external random number generator
+c in the case of weighted events
+c
+      subroutine thrbdec2(xext,imode,iunwgt,iseed,ntr,rnd1,rnd2,rnd3,
+     #                    rnd4,rnd5,xmw,xtq,xbq,xel,xnu,wgt)
+c Inputs:
+c   xext=.true. -> uses random numbers rnd1..rnd5, and ignores iseed
+c   imode=0 -> phase space only
+c   imode=1 -> top decay matrix element * phase space
+c   iunwgt=0 -> weighted events
+c   iunwgt=1 -> unweighted events
+c   iseed: random number seed
+c   xmw: W mass
+c   xtq: top four momentum
+c Outputs:
+c   xbq,xel,xnu: b quark, charged lepton, neutrino four momenta
+c   wgt: event weight; set to one if iunwgt=1
+c   ntr: the number of trials in the unweighted-event mode
+c The four momenta are arrays mom(i), with i=1,..,5, corresponding
+c to (px,py,pz,E,m)
+c
+c The matrix element squared for top decay, averaged over top spin, is
+c   |M|^2=2*g_W^4*|V_tb|^2*(top.el)*(b.nu)/( (Q2-mw^2)^2+(mw*wwidth)^2 )
+c (where Q2=2*el.nu is the W invariant mass), and the three-body phase-space
+c   dPhi_3=dQ2/(2*pi)*dPhi_2(t;W,b)*dPhi_2(W;el,nu)
+c   dPhi_2(t;W,b)=1/(8*(2*pi)**2)*(mt**2-Q2)/mt**2
+c   dPhi_2(W;el,nu)=1/(8*(2*pi)**2)
+c with the two two-body phase spaces computed in the t and W rest frames
+c respectively. In the weighted-event mode, the factor g_W^4*|V_tb|^2
+c is NOT included in wgt.
+c Since the phase-space variables are generated within this routine,
+c and only a seed or random numbers are given in input, the weight wgt 
+c includes the jacobian relevant to the definition of these variables
+      implicit none
+      logical xext
+      integer imode,iunwgt,iseed,ntr
+      real*8 rnd1,rnd2,rnd3,rnd4,rnd5,xmw,wgt,xtq(5),
+     # xbq(5),xel(5),xnu(5)
+      real*8 pi,pi2,wwidth,tiny,xmt,xmt2,xmw2,xmwlow2,xmwupp2,phspfact,
+     # bwmdpl,bwmdmn,bwfmpl,bwfmmn,bwdelf,thbdrandom,cth1,phi1,cth2,
+     # phi2,qw2,qw,drk2,thrbdec_dot,dpk1,decwgt,thrbdec_ibw,evwgt,
+     # bound,rat,rrnd,xx(5),xtq0(5),xbq0(5),xel0(5),xnu0(5)
+      integer i,ilocseed,maxtry