Merge lp:~maddevelopers/mg5amcnlo/3.0.2 into lp:mg5amcnlo

[Olivier Mattelaer]

This is mainly to start the discussion about releasing 3.0.2.
Are we waiting for something specific to release this?

[Olivier Mattelaer (olivier-mattelaer)]

First question,

Is it normal that the file tests/acceptance_tests/test_export_fks.py was removed?

Cheers,

Olivier

[marco zaro (marco-zaro)]

Hi Olivier,
I see that
 tests/unit_tests/iolibs/test_export_fks.py
exists, I don’t remember why and when it has been moved from/to acceptance_tests.

Ciao,

Marco

> On 27 Mar 2020, at 21:54, Olivier Mattelaer <email address hidden> wrote:
>
> First question,
>
> Is it normal that the file tests/acceptance_tests/test_export_fks.py was removed?
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.

[Olivier Mattelaer (olivier-mattelaer)]

Ok, so should we release this?

Olivier

[marco zaro (marco-zaro)]

Hi,
I am fine with releasing 3.0.2.
One detail: the new polyfit functionality included in the virt_tricks most of the time significantly increases the disk/ram usage, without having a major impact on run time. I tried it for a couple of processes with eternal photons (in the branch for the alpha0 scheme) and in all cases I had to disable this function via FKS_params.dat, specially when many jobs were submitted on the cluster.
So, I’d suggest either to clearly spell that there may be downsides on the memory/disk consumption, or to have the functionality disabled by default (sorry Rik for all your work…).

What does everybody think here?

Cheers,

marco

> On 6 Apr 2020, at 16:54, Olivier Mattelaer <email address hidden> wrote:
>
> Ok, so should we release this?
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.

[lp:~maddevelopers/mg5amcnlo/3.0.2]

980. By olivier-mattelaer on 2020-04-08

fixing tests

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

I do not know too much about that feature but presented like that this seems indeed that deactivating this by default make sense.

Otherwise, we have issue with the loop in this version:

======================================================================
FAIL: test_short_ML5_sm_vs_stored_ML4 (tests.parallel_tests.test_ML5.ML5Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/omattelaer/Documents/workspace/3.0.2/tests/parallel_tests/test_ML5.py", line 282, in test_short_ML5_sm_vs_stored_ML4
    chosen_runner='ML4')
  File "/Users/omattelaer/Documents/workspace/3.0.2/tests/parallel_tests/test_ML5.py", line 234, in compare_processes
    my_comp.assert_processes(self, tolerance)
  File "/Users/omattelaer/Documents/workspace/3.0.2/tests/parallel_tests/loop_me_comparator.py", line 1526, in assert_processes
    test_object.assertEqual(fail_str, "Failed for processes:")
AssertionError: 'Failed for processes:u u~ > d d~ 4.8307135319e+00 -3.4038975349e+005.7713195566e+00 d g > d g 1.1009861074e+01 2.7752500246e+00 5.9735543592e-01 g g > d d~ 9.9155829338e-01 -7.2430527735e+001.3172228006e+00 g g > t t~ 4.8445035533e+01 4.0080931226e+01 9.4481931273e-02 d~ d > g a -1.2494442497e+00-5.3667497832e+006.2230725293e-01 d~ d > g z -1.1932658581e+00-5.3105713960e+006.3305789752e-01 ' != 'Failed for processes:'

======================================================================
FAIL: test_short_ML5_sm_vs_stored_ML5 (tests.parallel_tests.test_ML5.ML5Test)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/omattelaer/Documents/workspace/3.0.2/tests/parallel_tests/test_ML5.py", line 274, in test_short_ML5_sm_vs_stored_ML5
    chosen_runner='ML5_opt')
  File "/Users/omattelaer/Documents/workspace/3.0.2/tests/parallel_tests/test_ML5.py", line 234, in compare_processes
    my_comp.assert_processes(self, tolerance)
  File "/Users/omattelaer/Documents/workspace/3.0.2/tests/parallel_tests/loop_me_comparator.py", line 1526, in assert_processes
    test_object.assertEqual(fail_str, "Failed for processes:")
AssertionError: 'Failed for processes:u u~ > d d~ 6.7303698120e+00 -1.5042412548e+001.5756617880e+00 d g > d g 1.9336116755e-01 -8.0412499551e+001.0492772446e+00 d~ u~ > d~ u~ 2.7814963527e+01 1.9580352460e+01 1.7374314097e-01 g u~ > g u~ 1.9336107775e-01 -8.0412499551e+001.0492772211e+00 g g > d d~ -1.2158385476e+00-9.4504496097e+007.7202218247e-01 g g > t t~ 4.6281795703e+01 3.7816621085e+01 1.0065795459e-01 g g > g g -1.0185218190e+01-1.8419829702e+012.8787267000e-01 d~ d > g a -2.0657523817e+00-6.1830579151e+004.9913931649e-01 u~ u > g z -2.1407021072e+00-6.2580076405e+004.9023072078e-01 d u~ > w- g -2.0691553635e+00-6.1864608969e+004.9872782401e-01 ' != 'Failed for processes:'

----------------------------------------------------------------------
Ran 2 tests in 1231.280s

FAILED (failures=2)
test_short_ML5_sm_vs_stored_ML4 test_short_ML5_sm_vs_stored_ML5

only the finite part is not in agreement...
Is this expected? ...

[Olivier Mattelaer (olivier-mattelaer)]

On a related note,
The python3 version of 3.0.2 seems to be ready:
lp:~maddevelopers/mg5amcnlo/3.0.2.py3

It obviously has the same issue with the MadLoop issue.
But otherwise it seems working as expected. So we can test it as well and move to a 4 release scheme in parallel.

I do not get any negative feedback on the python3 release (for the 2.x version) so I hope that we can quickly move --in a month or two-- to supporting only those python3 version (and simplify our lives). Remember that those python3 version still supports python2.7.

Cheers,

Olivier

[Richard Ruiz (rruiz)]

Hi Olivier,

Any feedback from users on how mg5mc@py3 is functioning? If not, I can ask
a few more people to stress-test.

best

On Thu, Apr 9, 2020 at 10:55 AM Olivier Mattelaer <
<email address hidden>> wrote:

> On a related note,
> The python3 version of 3.0.2 seems to be ready:
> lp:~maddevelopers/mg5amcnlo/3.0.2.py3
>
> It obviously has the same issue with the MadLoop issue.
> But otherwise it seems working as expected. So we can test it as well and
> move to a 4 release scheme in parallel.
>
> I do not get any negative feedback on the python3 release (for the 2.x
> version) so I hope that we can quickly move --in a month or two-- to
> supporting only those python3 version (and simplify our lives). Remember
> that those python3 version still supports python2.7.
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>

--
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Richard Ruiz, PhD
Fonds Speciaux de Recherche (FSR) Fellow,
Centre for Cosmology, Particle Physics and Phenomenology (CP3),
Université catholique de Louvain
Escargot: 2, Chemin du Cyclotron - Box L7.01.05,
B-1348 Louvain-la-Neuve, Belgium, Earth
Web: https://cp3.irmp.ucl.ac.be/Members/rruiz
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

[marco zaro (marco-zaro)]

Hi all,
any update on this merge?
Cheers,
Marco

On Thu, Apr 9, 2020 at 1:18 PM Richard Ruiz <email address hidden>
wrote:

> Hi Olivier,
>
> Any feedback from users on how mg5mc@py3 is functioning? If not, I can ask
> a few more people to stress-test.
>
> best
>
>
> On Thu, Apr 9, 2020 at 10:55 AM Olivier Mattelaer <
> <email address hidden>> wrote:
>
> > On a related note,
> > The python3 version of 3.0.2 seems to be ready:
> > lp:~maddevelopers/mg5amcnlo/3.0.2.py3
> >
> > It obviously has the same issue with the MadLoop issue.
> > But otherwise it seems working as expected. So we can test it as well and
> > move to a 4 release scheme in parallel.
> >
> > I do not get any negative feedback on the python3 release (for the 2.x
> > version) so I hope that we can quickly move --in a month or two-- to
> > supporting only those python3 version (and simplify our lives). Remember
> > that those python3 version still supports python2.7.
> >
> > Cheers,
> >
> > Olivier
> > --
> > https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> > Your team MadDevelopers is subscribed to branch
> > lp:~maddevelopers/mg5amcnlo/3.0.2.
> >
>
>
> --
> -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
> Richard Ruiz, PhD
> Fonds Speciaux de Recherche (FSR) Fellow,
> Centre for Cosmology, Particle Physics and Phenomenology (CP3),
> Université catholique de Louvain
> Escargot: 2, Chemin du Cyclotron - Box L7.01.05,
> B-1348 Louvain-la-Neuve, Belgium, Earth
> Web: https://cp3.irmp.ucl.ac.be/Members/rruiz
> -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
>
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadTeam is requested to review the proposed merge of
> lp:~maddevelopers/mg5amcnlo/3.0.2 into lp:mg5amcnlo/3.0.
>

[lp:~maddevelopers/mg5amcnlo/3.0.2]

981. By olivier-mattelaer on 2020-05-11

merge up to 2.7.2

[Olivier Mattelaer (olivier-mattelaer)]

I have a meeting with Valentin to discuss this merging this afternoon.
Otherwise I have tracked the issue to revision 974
which performs the merge with this branch:
lp:~maddevelopers/mg5amcnlo/3.0.2_loopparams

[lp:~maddevelopers/mg5amcnlo/3.0.2]

982. By olivier-mattelaer on 2020-05-11

fix issue with the non update of loop parameter

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

I have pushed a fix to the issue (thanks Valentin for the discussion on this).
Please take a look if you agree with such modification.

I will move forward on releasing this version now.

Cheers,

Olivier

[marco zaro (marco-zaro)]

Hi Olivier,
thanks for the fix.
how about the polyfit business?
Cheers,

Marco

On Wed, Apr 8, 2020 at 10:29 AM Marco Zaro <email address hidden> wrote:

> Hi,
> I am fine with releasing 3.0.2.
> One detail: the new polyfit functionality included in the virt_tricks most
> of the time significantly increases the disk/ram usage, without having a
> major impact on run time. I tried it for a couple of processes with eternal
> photons (in the branch for the alpha0 scheme) and in all cases I had to
> disable this function via FKS_params.dat, specially when many jobs were
> submitted on the cluster.
> So, I’d suggest either to clearly spell that there may be downsides on the
> memory/disk consumption, or to have the functionality disabled by default
> (sorry Rik for all your work…).
>
> What does everybody think here?
>
> Cheers,
>
> marco
>
> > On 6 Apr 2020, at 16:54, Olivier Mattelaer <email address hidden>
> wrote:
> >
> > Ok, so should we release this?
> >
> > Olivier
> > --
> > https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> > Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>
>

[Olivier Mattelaer (olivier-mattelaer)]

Should we disable it by default?

Cheers,

Olivier

> On 11 May 2020, at 19:06, marco zaro <email address hidden> wrote:
>
> Hi Olivier,
> thanks for the fix.
> how about the polyfit business?
> Cheers,
>
> Marco
>
>
> On Wed, Apr 8, 2020 at 10:29 AM Marco Zaro <email address hidden> wrote:
>
>> Hi,
>> I am fine with releasing 3.0.2.
>> One detail: the new polyfit functionality included in the virt_tricks most
>> of the time significantly increases the disk/ram usage, without having a
>> major impact on run time. I tried it for a couple of processes with eternal
>> photons (in the branch for the alpha0 scheme) and in all cases I had to
>> disable this function via FKS_params.dat, specially when many jobs were
>> submitted on the cluster.
>> So, I’d suggest either to clearly spell that there may be downsides on the
>> memory/disk consumption, or to have the functionality disabled by default
>> (sorry Rik for all your work…).
>>
>> What does everybody think here?
>>
>> Cheers,
>>
>> marco
>>
>>> On 6 Apr 2020, at 16:54, Olivier Mattelaer <email address hidden>
>> wrote:
>>>
>>> Ok, so should we release this?
>>>
>>> Olivier
>>> --
>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
>>> Your team MadDevelopers is subscribed to branch
>> lp:~maddevelopers/mg5amcnlo/3.0.2.
>>
>>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.

[marco zaro (marco-zaro)]

This is what I would suggest, ie to set

158 #UsePolyVirtual

159 .False.

inside Template/NLO/Cards/FKS_params.dat

cheers,

Marco Zaro

On Mon, May 11, 2020 at 8:04 PM Olivier Mattelaer <
<email address hidden>> wrote:

> Should we disable it by default?
>
> Cheers,
>
> Olivier
>
> > On 11 May 2020, at 19:06, marco zaro <email address hidden> wrote:
> >
> > Hi Olivier,
> > thanks for the fix.
> > how about the polyfit business?
> > Cheers,
> >
> > Marco
> >
> >
> > On Wed, Apr 8, 2020 at 10:29 AM Marco Zaro <email address hidden> wrote:
> >
> >> Hi,
> >> I am fine with releasing 3.0.2.
> >> One detail: the new polyfit functionality included in the virt_tricks
> most
> >> of the time significantly increases the disk/ram usage, without having a
> >> major impact on run time. I tried it for a couple of processes with
> eternal
> >> photons (in the branch for the alpha0 scheme) and in all cases I had to
> >> disable this function via FKS_params.dat, specially when many jobs were
> >> submitted on the cluster.
> >> So, I’d suggest either to clearly spell that there may be downsides on
> the
> >> memory/disk consumption, or to have the functionality disabled by
> default
> >> (sorry Rik for all your work…).
> >>
> >> What does everybody think here?
> >>
> >> Cheers,
> >>
> >> marco
> >>
> >>> On 6 Apr 2020, at 16:54, Olivier Mattelaer <
> <email address hidden>>
> >> wrote:
> >>>
> >>> Ok, so should we release this?
> >>>
> >>> Olivier
> >>> --
> >>>
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> >>> Your team MadDevelopers is subscribed to branch
> >> lp:~maddevelopers/mg5amcnlo/3.0.2.
> >>
> >>
> >
> > --
> > https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> > Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>

[Olivier Mattelaer (olivier-mattelaer)]

Rikkert can you comment on that?

Olivier

[Rikkert Frederix (frederix)]

Hi,

It's kind of a pity that this PolyVirtual doesn't work for you. Let's keep it turned off by default as Marco suggests -- at least for this release.

Best,
Rikkert

[marco zaro (marco-zaro)]

Hi Rikkert,
thanks for your work on this. As I told you, I had crashes on the louvain cluster, due to the large memory and space usage. I guess that if you have few jobs that can be run on multicore, then you do get improvements, but i think that on a cluster it can create problems…

cheers,

marco

> On 13 May 2020, at 17:42, Rikkert Frederix <email address hidden> wrote:
>
> Hi,
>
> It's kind of a pity that this PolyVirtual doesn't work for you. Let's keep it turned off by default as Marco suggests -- at least for this release.
>
> Best,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.

[Valentin Hirschi (valentin-hirschi)]

Very naively I would expect the source of large memory consumption as to be
stemming from either a memory leak of PolyFit unnecessarily keeping all
training sample fed in memory.
Rik, do you know which one it is? If the latter, isn't there any easy way
of truncating the internal "history" of PolyFit so as to effectively cap
its memory footprint?

On Wed, May 13, 2020 at 5:45 PM marco zaro <email address hidden> wrote:

> Hi Rikkert,
> thanks for your work on this. As I told you, I had crashes on the louvain
> cluster, due to the large memory and space usage. I guess that if you have
> few jobs that can be run on multicore, then you do get improvements, but i
> think that on a cluster it can create problems…
>
> cheers,
>
> marco
>
>
> > On 13 May 2020, at 17:42, Rikkert Frederix <email address hidden>
> wrote:
> >
> > Hi,
> >
> > It's kind of a pity that this PolyVirtual doesn't work for you. Let's
> keep it turned off by default as Marco suggests -- at least for this
> release.
> >
> > Best,
> > Rikkert
> >
> > --
> > https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> > Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>

--
Valentin

[Rikkert Frederix (frederix)]

Hi Valentin,

PolyVirtual keeps all events (i.e., the random numbers that corresponds to the generated momenta) for which the virtual have been computed in memory to get the best possible polynomial fit for the virtuals. It also writes them all to disk so that they can be re-used in the next integration step.

I had hoped that in practice there are few enough PS points for which the virtuals need to be computed (i.e., well below several million) that this would not become an issue for most processes (in particular since the virt-tricks should work slightly better with the PolyVirtuals than without). Apparently this is not good enough when running on a cluster requiring high accuracy.

Best,
Rikkert

[Valentin Hirschi (valentin-hirschi)]

Hi Rik,

Ok and I guess truncating the history to keep only a fixed number of past
events (may not impact the quality of the fit so much if this number is
large enough) would be hard to implement?

On Wed, May 13, 2020 at 6:25 PM Rikkert Frederix <email address hidden>
wrote:

> Hi Valentin,
>
> PolyVirtual keeps all events (i.e., the random numbers that corresponds to
> the generated momenta) for which the virtual have been computed in memory
> to get the best possible polynomial fit for the virtuals. It also writes
> them all to disk so that they can be re-used in the next integration step.
>
> I had hoped that in practice there are few enough PS points for which the
> virtuals need to be computed (i.e., well below several million) that this
> would not become an issue for most processes (in particular since the
> virt-tricks should work slightly better with the PolyVirtuals than
> without). Apparently this is not good enough when running on a cluster
> requiring high accuracy.
>
> Best,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>

--
Valentin

[Olivier Mattelaer (olivier-mattelaer)]

This looks like a very painfull version to release:
Looks like we have a new (MadLoop?) issue:

The (parralel) test crashing is the following one:

ERROR: test_short_jet_veto_xsec (tests.parallel_tests.test_cmd_amcatnlo.MECmdShell)
tests the jet-veto cross section at NNLL+NLO

The associated log file points an issue with the poles:
In this case ML returns 0 which likely points to another issue related to not updating the virtual coupling correctly.

 ---- POLES CANCELLED ----
  COEFFICIENT DOUBLE POLE:
        MadFKS: -0.43764352847099569 OLP: 0.0000000000000000
  COEFFICIENT SINGLE POLE:
        MadFKS: 1.1667772655032984 OLP: 0.0000000000000000
 FIX VETOXSEC in FKS_EW
Time in seconds: 1

[Eleni Vryonidou (evryonidou)]

Hi Olivier, we noticed this (or a very similar) problem with some EFT tests in 3.0.x. Gauthier reported this by email to Valentin. I can forward the emails to you.

The problem occurred with COLLIER, we had this kind of incoherent output produced:

 ---- POLES CANCELLED ----
  COEFFICIENT DOUBLE POLE:
        MadFKS: -8.8973918441133455E-008 OLP: 0.0000000000000000
  COEFFICIENT SINGLE POLE:
        MadFKS: -4.7574099859731733E-008 OLP: -5.1368048587786267E-008

[Rikkert Frederix (frederix)]

Hi Olivier,

The NNLL+NLO jet veto cross sections are not supposed to work yet in the 3.0.2 branch. Hence, it shouldn't be a surprise that this test fails.
However, the factor that the coefficients of the poles are not correct is strange and shouldn't happen.

Best,
Rikkert

[lp:~maddevelopers/mg5amcnlo/3.0.2]

983. By olivier-mattelaer on 2020-05-13

remove polyfit by default

984. By olivier-mattelaer on 2020-05-15

update IOttest

[Olivier Mattelaer (olivier-mattelaer)]

Indeed with ninja we have the correct pole:
 ---- POLES CANCELLED ----
  COEFFICIENT DOUBLE POLE:
 MadFKS: -0.43764352847099569 OLP: -0.43764352847158994
  COEFFICIENT SINGLE POLE:
 MadFKS: 1.1667772655032984 OLP: 1.1667772655048829

Valentin what is your advise on this?
Should we veto the use of Collier for the moment?

Cheers,

Olivier

[marco zaro (marco-zaro)]

Can we set ninja as default (as it should be) and use collier only for
recovery?
I would not disable it completely...
cheers,

Marco Zaro

On Wed, May 20, 2020 at 5:01 PM Olivier Mattelaer <
<email address hidden>> wrote:

> Indeed with ninja we have the correct pole:
> ---- POLES CANCELLED ----
> COEFFICIENT DOUBLE POLE:
> MadFKS: -0.43764352847099569 OLP:
> -0.43764352847158994
> COEFFICIENT SINGLE POLE:
> MadFKS: 1.1667772655032984 OLP:
> 1.1667772655048829
>
> Valentin what is your advise on this?
> Should we veto the use of Collier for the moment?
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>

[Olivier Mattelaer (olivier-mattelaer)]

I face the issue because ninja crashed at the installation actually.
Then collier became the default...

1) if we allow collier they are no guarantee that this is not the default.
2) Having a backup returning wrong value is not a good idea either...

I would personally veto collier here and check in 2.7.3 if we need to veto collier as well or not.
This as long as Valentin does not clarify the situation obviously. (and I do not want to wait for a fix to be ready to release this version)

Cheers,

Olivier

> On 20 May 2020, at 18:00, marco zaro <email address hidden> wrote:
>
> Can we set ninja as default (as it should be) and use collier only for
> recovery?
> I would not disable it completely...
> cheers,
>
> Marco Zaro
>
>
> On Wed, May 20, 2020 at 5:01 PM Olivier Mattelaer <
> <email address hidden>> wrote:
>
>> Indeed with ninja we have the correct pole:
>> ---- POLES CANCELLED ----
>> COEFFICIENT DOUBLE POLE:
>> MadFKS: -0.43764352847099569 OLP:
>> -0.43764352847158994
>> COEFFICIENT SINGLE POLE:
>> MadFKS: 1.1667772655032984 OLP:
>> 1.1667772655048829
>>
>> Valentin what is your advise on this?
>> Should we veto the use of Collier for the moment?
>>
>> Cheers,
>>
>> Olivier
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
>> Your team MadDevelopers is subscribed to branch
>> lp:~maddevelopers/mg5amcnlo/3.0.2.
>>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.

[Valentin Hirschi (valentin-hirschi)]

I think Ninja crashed because Tiziano moved the installation tarball as
part of his migration from bitbucket to github.

We must keep COLLIER as it is crucial for loop-induced efficiency.
I could not reproduce seg-faults in standalone (and neither did Richard I
believe).
So it is either due to an interplay with the MadFKS code or because of some
particular setting of MadLoop when run within MadFKS as opposed to run in
standalone.

On Wed, May 20, 2020 at 6:12 PM Olivier Mattelaer <
<email address hidden>> wrote:

> I face the issue because ninja crashed at the installation actually.
> Then collier became the default...
>
> 1) if we allow collier they are no guarantee that this is not the default.
> 2) Having a backup returning wrong value is not a good idea either...
>
>
> I would personally veto collier here and check in 2.7.3 if we need to veto
> collier as well or not.
> This as long as Valentin does not clarify the situation obviously. (and I
> do not want to wait for a fix to be ready to release this version)
>
> Cheers,
>
> Olivier
>
>
> > On 20 May 2020, at 18:00, marco zaro <email address hidden> wrote:
> >
> > Can we set ninja as default (as it should be) and use collier only for
> > recovery?
> > I would not disable it completely...
> > cheers,
> >
> > Marco Zaro
> >
> >
> > On Wed, May 20, 2020 at 5:01 PM Olivier Mattelaer <
> > <email address hidden>> wrote:
> >
> >> Indeed with ninja we have the correct pole:
> >> ---- POLES CANCELLED ----
> >> COEFFICIENT DOUBLE POLE:
> >> MadFKS: -0.43764352847099569 OLP:
> >> -0.43764352847158994
> >> COEFFICIENT SINGLE POLE:
> >> MadFKS: 1.1667772655032984 OLP:
> >> 1.1667772655048829
> >>
> >> Valentin what is your advise on this?
> >> Should we veto the use of Collier for the moment?
> >>
> >> Cheers,
> >>
> >> Olivier
> >> --
> >> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> >> Your team MadDevelopers is subscribed to branch
> >> lp:~maddevelopers/mg5amcnlo/3.0.2.
> >>
> >
> > --
> > https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> > Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>

--
Valentin

[Stefano Frixione (stefano-frixione)]

Sorry if I didn't follow, but if ninja crashes, we must fix ninja,
not promote collier to being the default.
I think we must not release till ninja is fixed.
Cheers, Stefano.

On Wed, 20 May 2020, Olivier Mattelaer wrote:

> I face the issue because ninja crashed at the installation actually.
> Then collier became the default...
>
> 1) if we allow collier they are no guarantee that this is not the default.
> 2) Having a backup returning wrong value is not a good idea either...
>
>
> I would personally veto collier here and check in 2.7.3 if we need to veto collier as well or not.
> This as long as Valentin does not clarify the situation obviously. (and I do not want to wait for a fix to be ready to release this version)
>
> Cheers,
>
> Olivier
>
>
>> On 20 May 2020, at 18:00, marco zaro <email address hidden> wrote:
>>
>> Can we set ninja as default (as it should be) and use collier only for
>> recovery?
>> I would not disable it completely...
>> cheers,
>>
>> Marco Zaro
>>
>>
>> On Wed, May 20, 2020 at 5:01 PM Olivier Mattelaer <
>> <email address hidden>> wrote:
>>
>>> Indeed with ninja we have the correct pole:
>>> ---- POLES CANCELLED ----
>>> COEFFICIENT DOUBLE POLE:
>>> MadFKS: -0.43764352847099569 OLP:
>>> -0.43764352847158994
>>> COEFFICIENT SINGLE POLE:
>>> MadFKS: 1.1667772655032984 OLP:
>>> 1.1667772655048829
>>>
>>> Valentin what is your advise on this?
>>> Should we veto the use of Collier for the moment?
>>>
>>> Cheers,
>>>
>>> Olivier
>>> --
>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
>>> Your team MadDevelopers is subscribed to branch
>>> lp:~maddevelopers/mg5amcnlo/3.0.2.
>>>
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
>> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.
>

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

Let's try to avoid any confusion here:
1) we do not have any issue with ninja.
   - we can still install it
   - when install it works

2) For some unknow reason, my version of 3.0.2 version was configured as if ninja was not installed in my machine, in that case collier is the next OLP provider and is used as default.

3) We have a bug related to collier that does not return the correct pole (I'm therefore worried about the other output like the finite part of the loop). This bug has been spotted by Eleni and Gauthier and independently by me (thanks to point 2 above)

>We must keep COLLIER as it is crucial for loop-induced efficiency.
I could not reproduce seg-faults in standalone (and neither did Richard I
believe).

Being fast is useless if we are wrong. Like this, this sounds like a bad bug.
If you are sure that this is a MadFKS related bug, then we can forbid the usage of COLLIER only in that case. But it would be nice that you take a real look at this. I prefer to be slow and correct than fast and wrong.

Cheers,

Olivier

[Eleni Vryonidou (evryonidou)]

Hi Olivier,

when I was looking at this a few weeks ago I realised that in BinothLHA.f there is this line:

 if ((dabs(avgPoleRes(1))+dabs(avgPoleRes(2))).ne.0d0) then
               cpol = .not.((((PoleDiff(1)+PoleDiff(2))/
     $ (dabs(avgPoleRes(1))+dabs(avgPoleRes(2)))) .lt.
     $ tolerance*10d0).or.(mod(ret_code,10).eq.7))

and
               cpol = .not.((PoleDiff(1)+PoleDiff(2).lt.tolerance*10d0)
     $ .or.(mod(ret_code,10).eq.7))

I believe Collier corresponds to .eq.7 and this is why the test seems to pass, even if the numbers are not agreeing. E.g. I see things like:

 ---- POLES CANCELLED ----
  COEFFICIENT DOUBLE POLE:
        MadFKS: -5.1118724822382218E-008 OLP: 0.0000000000000000
  COEFFICIENT SINGLE POLE:
        MadFKS: -5.4796083693610551E-008 OLP: -9.2617064392846240E-008

Maybe someone can confirm why that line was added there. It’s not the source of the zero numbers, but it’s weird that this was set up like that to allow the test to pass.

Cheers,

Eleni

> On 20 May 2020, at 20:38, Olivier Mattelaer <email address hidden> wrote:
>
> Hi,
>
> Let's try to avoid any confusion here:
> 1) we do not have any issue with ninja.
> - we can still install it
> - when install it works
>
> 2) For some unknow reason, my version of 3.0.2 version was configured as if ninja was not installed in my machine, in that case collier is the next OLP provider and is used as default.
>
> 3) We have a bug related to collier that does not return the correct pole (I'm therefore worried about the other output like the finite part of the loop). This bug has been spotted by Eleni and Gauthier and independently by me (thanks to point 2 above)
>
>
>> We must keep COLLIER as it is crucial for loop-induced efficiency.
> I could not reproduce seg-faults in standalone (and neither did Richard I
> believe).
>
> Being fast is useless if we are wrong. Like this, this sounds like a bad bug.
> If you are sure that this is a MadFKS related bug, then we can forbid the usage of COLLIER only in that case. But it would be nice that you take a real look at this. I prefer to be slow and correct than fast and wrong.
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.

[Olivier Mattelaer (olivier-mattelaer)]

Those lines were introduced by rikkert here:
https://bazaar.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/revision/841#Template/NLO/SubProcesses/BinothLHA.f

But if you check the previous version [840] you face the exact same issue.

Actually I reproduce this issue up to version 835.
Version before that are crashing in another way preventing me to get more insight on this issue.
This means that both 3.0.0 and 3.0.1 have the exact same issue.

Cheers,

Olivier

[Rikkert Frederix (frederix)]

Hi guys,

These lines were NOT introduced by me. If you go back in through the history you find them all the way back to some merge with a 'COLLIER' branch, (which was merged into 2.5.1 (hence, this is also part of public 2.7.X code) and that branch was later merged to the flattened_FKS_EW_clean branch, which was merged into the current 3.0.X branch).

The COLLIER branch is a branch that Valentin worked on some three years ago:

https://bazaar.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/changes/78.355.27

They were introduced in revision 78.355.19. I couldn't find a commit message, so I don't know why this was introduced here. Valentin, do you remember?

Best,
Rikkert

[Stefano Frixione (stefano-frixione)]

Hi Olivier,
thanks for the clarification.

> If you are sure that this is a MadFKS related bug, then we can forbid
> the usage of COLLIER only in that case. But it would be nice that you
A bug is a bug, but I find it hard to believe it's in MadFKS, given
that in this case the only information returned by the latter is the
residues of the poles -- these are of course independent of the OLP,
so one might verify if there is a problem in that respect by running
the same thing on a machine where ninja works.

Cheers, Stefano.

> take a real look at this. I prefer to be slow and correct than fast and
> wrong.
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.
>

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

Sorry indeed I miss the fact that Rik change the line in that version but only in a cosmetic way.
Consequently to this I have force the printout of the poles in 2.7.3 and found the same printout:

 ---- POLES CANCELLED ----
 POLES ANCELLATION, DIFFERENCE < 1.0000000000000003E-005
  COEFFICIENT DOUBLE POLE:
        MadFKS: -0.48049656339581920 OLP: 0.0000000000000000
  COEFFICIENT SINGLE POLE:
        MadFKS: 1.2727761694678807 OLP: 0.0000000000000000
  FINITE:
           OLP: -0.69598571535491494

  MOMENTA (Exyzm):
           1 39.800372893656203 0.0000000000000000 0.0000000000000000 39.800372893656203 \
  0.0000000000000000
           2 39.800372893656203 -0.0000000000000000 -0.0000000000000000 -39.800372893656203 \
  0.0000000000000000
           3 39.800372893656203 -18.912294290529488 -32.109690658036719 -13.977216210235818 \
  0.0000000000000000
           4 39.800372893656203 18.912294290529488 32.109690658036719 13.977216210235818 \
  0.0000000000000000

Now running in standalone mode on the same PS point we have (with collier):
|| Total(*) Born contribution (GeV^0):
| Born = 6.2823500530737046e+00
|| Total(*) virtual contribution normalized with born*alpha_S/(2*pi):
| Finite = 1.3071069893544380e+00
| Single pole = -4.7247991865699230e+00
| Double pole = -2.6666666666666692e+00
| (*) The results above sum all starred contributions below

Since Valentin does not want to take a look at this and that I'm completely loss on this I have:
1) revert the bypass on the check of the poles for COLLIER
2) set the default for madloop to ninja then cuttools (collier is only on request)
I have done that both in 2.7.3 and 3.0.2
As soon as Valentin found the time to fix this, we can revert this obviously but I have spent too much time on this and want ( and need to focus on something else)

So with this, we are ready to release 3.0.2

Cheers,

Olivier

review: Approve

[lp:~maddevelopers/mg5amcnlo/3.0.2]

985. By olivier-mattelaer on 2020-05-21

deactivating Collier as long as Valentin did not take a look

[Valentin Hirschi (valentin-hirschi)]

Sorry for my slow reaction, I will have a serious look at this.
I'll check with Olivier what is the exact setup to reproduce the problem.
The fact that the problem cannot be reproduced in standalone (unless I
misunderstood the situation, in which case apologies) made me kind of put
it aside in so far as it did not seem like a genuine MadLoop problem but
more its use within MadFKS.
Hopefully I'll be able to easily spot the offending parameter by simply
comparing MadLoop's run for one particular point within MadFKS and in
standalone (which is what I was hoping would straight-forwardly reveal the
issue).
But it seems that the problem is more involved than I thought so I'll
investigate and let you knw.

On Thu, May 21, 2020 at 3:57 PM Olivier Mattelaer <
<email address hidden>> wrote:

> Review: Approve
>
> Hi,
>
> Sorry indeed I miss the fact that Rik change the line in that version but
> only in a cosmetic way.
> Consequently to this I have force the printout of the poles in 2.7.3 and
> found the same printout:
>
> ---- POLES CANCELLED ----
> POLES ANCELLATION, DIFFERENCE < 1.0000000000000003E-005
> COEFFICIENT DOUBLE POLE:
> MadFKS: -0.48049656339581920 OLP:
> 0.0000000000000000
> COEFFICIENT SINGLE POLE:
> MadFKS: 1.2727761694678807 OLP:
> 0.0000000000000000
> FINITE:
> OLP: -0.69598571535491494
>
> MOMENTA (Exyzm):
> 1 39.800372893656203 0.0000000000000000
> 0.0000000000000000 39.800372893656203 \
> 0.0000000000000000
> 2 39.800372893656203 -0.0000000000000000
> -0.0000000000000000 -39.800372893656203 \
> 0.0000000000000000
> 3 39.800372893656203 -18.912294290529488
> -32.109690658036719 -13.977216210235818 \
> 0.0000000000000000
> 4 39.800372893656203 18.912294290529488
> 32.109690658036719 13.977216210235818 \
> 0.0000000000000000
>
> Now running in standalone mode on the same PS point we have (with collier):
> || Total(*) Born contribution (GeV^0):
> | Born = 6.2823500530737046e+00
> || Total(*) virtual contribution normalized with born*alpha_S/(2*pi):
> | Finite = 1.3071069893544380e+00
> | Single pole = -4.7247991865699230e+00
> | Double pole = -2.6666666666666692e+00
> | (*) The results above sum all starred contributions below
>
>
> Since Valentin does not want to take a look at this and that I'm
> completely loss on this I have:
> 1) revert the bypass on the check of the poles for COLLIER
> 2) set the default for madloop to ninja then cuttools (collier is only on
> request)
> I have done that both in 2.7.3 and 3.0.2
> As soon as Valentin found the time to fix this, we can revert this
> obviously but I have spent too much time on this and want ( and need to
> focus on something else)
>
> So with this, we are ready to release 3.0.2
>
> Cheers,
>
> Olivier
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/3.0.2.
>

--
Valentin

[Valentin Hirschi (valentin-hirschi)]

So the problem seem to simply be that these few lines in `Binoth_LHA.f`:

         CALL COLLIER_COMPUTE_UV_POLES(.FALSE.)
         CALL COLLIER_COMPUTE_IR_POLES(.FALSE.)

are called *after* the similar ones in `check_poles.f`:

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.)

which of course means you'd get zero poles when running with COLLIER only.

The introduction of a common block "flag" of the sort of:

logical IS_DOING_POLE_CHECK
common pole_check/IS_DOING_POLE_CHECK
data IS_DOING_POLE_CHECK/.FALSE./

would allow to solve the problem by using it to make sure that
`Binoth_LHA.f` does not disable pole computations with COLLIER when doing
pole checks.

On Thu, May 21, 2020 at 9:42 PM Valentin Hirschi <email address hidden>
wrote:

> Sorry for my slow reaction, I will have a serious look at this.
> I'll check with Olivier what is the exact setup to reproduce the problem.
> The fact that the problem cannot be reproduced in standalone (unless I
> misunderstood the situation, in which case apologies) made me kind of put
> it aside in so far as it did not seem like a genuine MadLoop problem but
> more its use within MadFKS.
> Hopefully I'll be able to easily spot the offending parameter by simply
> comparing MadLoop's run for one particular point within MadFKS and in
> standalone (which is what I was hoping would straight-forwardly reveal the
> issue).
> But it seems that the problem is more involved than I thought so I'll
> investigate and let you knw.
>
> On Thu, May 21, 2020 at 3:57 PM Olivier Mattelaer <
> <email address hidden>> wrote:
>
> > Review: Approve
> >
> > Hi,
> >
> > Sorry indeed I miss the fact that Rik change the line in that version but
> > only in a cosmetic way.
> > Consequently to this I have force the printout of the poles in 2.7.3 and
> > found the same printout:
> >
> > ---- POLES CANCELLED ----
> > POLES ANCELLATION, DIFFERENCE < 1.0000000000000003E-005
> > COEFFICIENT DOUBLE POLE:
> > MadFKS: -0.48049656339581920 OLP:
> > 0.0000000000000000
> > COEFFICIENT SINGLE POLE:
> > MadFKS: 1.2727761694678807 OLP:
> > 0.0000000000000000
> > FINITE:
> > OLP: -0.69598571535491494
> >
> > MOMENTA (Exyzm):
> > 1 39.800372893656203 0.0000000000000000
> > 0.0000000000000000 39.800372893656203 \
> > 0.0000000000000000
> > 2 39.800372893656203 -0.0000000000000000
> > -0.0000000000000000 -39.800372893656203 \
> > 0.0000000000000000
> > 3 39.800372893656203 -18.912294290529488
> > -32.109690658036719 -13.977216210235818 \
> > 0.0000000000000000
> > 4 39.800372893656203 18.912294290529488
> > 32.109690658036719 13.977216210235818 \
> > 0.0000000000000000
> >
> > Now running in standalone mode on the same PS point we have (with
> collier):
> > || Total(*) Born contribution (GeV^0):
> > | Born = 6.2823500530737046e+00
> > || Total(*) virtual contrib...

[marco zaro (marco-zaro)]

Hi all,
just noticed a possible bug:

inside SubProcesses/cuts.f, line 142 reads

                     if (ipdg(i).eq.-ipdg(j)) then

however, it should be ipdg_reco in both cases and not ipdg, right? (that
is, the pdg's of the particles after the recombination with photons...)
If anybody could confirm, I will push the fix.

Best wishes,

Marco Zaro

On Thu, May 21, 2020 at 11:42 PM Valentin Hirschi <
<email address hidden>> wrote:

> So the problem seem to simply be that these few lines in `Binoth_LHA.f`:
>
> CALL COLLIER_COMPUTE_UV_POLES(.FALSE.)
> CALL COLLIER_COMPUTE_IR_POLES(.FALSE.)
>
> are called *after* the similar ones in `check_poles.f`:
>
> 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.)
>
> which of course means you'd get zero poles when running with COLLIER only.
>
> The introduction of a common block "flag" of the sort of:
>
> logical IS_DOING_POLE_CHECK
> common pole_check/IS_DOING_POLE_CHECK
> data IS_DOING_POLE_CHECK/.FALSE./
>
> would allow to solve the problem by using it to make sure that
> `Binoth_LHA.f` does not disable pole computations with COLLIER when doing
> pole checks.
>
> On Thu, May 21, 2020 at 9:42 PM Valentin Hirschi <
> <email address hidden>>
> wrote:
>
> > Sorry for my slow reaction, I will have a serious look at this.
> > I'll check with Olivier what is the exact setup to reproduce the problem.
> > The fact that the problem cannot be reproduced in standalone (unless I
> > misunderstood the situation, in which case apologies) made me kind of put
> > it aside in so far as it did not seem like a genuine MadLoop problem but
> > more its use within MadFKS.
> > Hopefully I'll be able to easily spot the offending parameter by simply
> > comparing MadLoop's run for one particular point within MadFKS and in
> > standalone (which is what I was hoping would straight-forwardly reveal
> the
> > issue).
> > But it seems that the problem is more involved than I thought so I'll
> > investigate and let you knw.
> >
> > On Thu, May 21, 2020 at 3:57 PM Olivier Mattelaer <
> > <email address hidden>> wrote:
> >
> > > Review: Approve
> > >
> > > Hi,
> > >
> > > Sorry indeed I miss the fact that Rik change the line in that version
> but
> > > only in a cosmetic way.
> > > Consequently to this I have force the printout of the poles in 2.7.3
> and
> > > found the same printout:
> > >
> > > ---- POLES CANCELLED ----
> > > POLES ANCELLATION, DIFFERENCE < 1.0000000000000003E-005
> > > COEFFICIENT DOUBLE POLE:
> > > MadFKS: -0.48049656339581920 OLP:
> > > 0.0000000000000000
> > > COEFFICIENT SINGLE POLE:
> > > MadFKS: 1.2727761694678807 OLP:
> > > 0.0000000000000000
> > > FINITE:
> > > OLP: -0.69598571535491494
> > >
> > > MOMENTA (Exyzm):
> > > 1 39.800372893656203 0.0000000000000000
> > > 0.0000000000000000 39.800372893656203 \
> > > 0.0000000000000000
> > > 2 39.800372893656203 -0.0000000000000000
>...

[lp:~maddevelopers/mg5amcnlo/3.0.2]

986. By olivier-mattelaer on 2020-05-23

fixing issuue with the pole check with COLLIER

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

Thanks I have pushed such changes, a variable was actually already existing so I used that one.

@Marco, by reading the code, I would agree with you that this is a bug.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/3.0.2]

987. By olivier-mattelaer on 2020-05-23

fix cuts.f --issue found by marco--

[marco zaro (marco-zaro)]

Thanks for the fix.
Cheers,

Marco

> On 23 May 2020, at 21:27, Olivier Mattelaer <email address hidden> wrote:
>
> Hi Valentin,
>
> Thanks I have pushed such changes, a variable was actually already existing so I used that one.
>
> @Marco, by reading the code, I would agree with you that this is a bug.
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.2/+merge/381317
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/3.0.2.

=== modified file 'HELAS/v2xxxx.f'
--- HELAS/v2xxxx.f	2010-10-30 03:26:37 +0000
+++ HELAS/v2xxxx.f	2020-05-23 19:29:41 +0000
@@ -30,11 +30,11 @@
             vc(2) = p(1)/p(0)
             vc(3) = p(2)/p(0)
             vc(4) = p(3)/p(0)
-         else
+      else
             vc(1) = p(0)/vmass
             vc(2) = p(1)/vmass
             vc(3) = p(2)/vmass
             vc(4) = p(3)/vmass
-         endif
+      endif
       return
       end

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2018-09-21 13:46:51 +0000
+++ MadSpin/decay.py	2020-05-23 19:29:41 +0000
@@ -284,12 +284,13 @@
         line_type = 'none' # support type: init / event / rwgt
         self.diese = ''
         for line in self.inputfile:
+            origline = line
             line = line.lower()
             if line=="":
                 continue 
             # Find special tag in the line
             if line[0]=="#":
-                self.diese+=line
+                self.diese+=origline
                 continue
             if '<event' in line:
                 #start new_event
@@ -319,7 +320,7 @@
             if line_type == 'none':
                 continue
             elif line_type == 'other_block':
-                self.diese += line
+                self.diese += origline
             # read the line and assign the date accordingly                
             elif line_type == 'init':
                 line_type = 'event'
@@ -2244,7 +2245,7 @@
             nb_mc_masses=len(indices_for_mc_masses)
 
             p, p_str=self.curr_event.give_momenta(event_map)
-            stdin_text=' %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, decay_me['max_weight'])
+            stdin_text=' %s %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, decay_me['max_weight'], self.options['frame_id'])
             stdin_text+=p_str
             # here I also need to specify the Monte Carlo Masses
             stdin_text+=" %s \n" % nb_mc_masses
@@ -2372,7 +2373,7 @@
         #no decays for this production mode, run in passthrough mode, only adding the helicities to the events
         nb_mc_masses=0
         p, p_str=self.curr_event.give_momenta(event_map)
-        stdin_text=' %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, 1.0)
+        stdin_text=' %s %s %s %s \n' % ('2', self.options['BW_cut'], self.Ecollider, 1.0, self.options['frame_id'])
         stdin_text+=p_str
         # here I also need to specify the Monte Carlo Masses
         stdin_text+=" %s \n" % nb_mc_masses
@@ -3277,8 +3278,7 @@
         """return the max. weight associated with me decay['path']"""
 
         p, p_str=self.curr_event.give_momenta(event_map)
-
-        std_in=" %s  %s %s %s  \n" % ("1",BWcut, self.Ecollider, nbpoints)
+        std_in=" %s  %s %s %s %s \n" % ("1",BWcut, self.Ecollider, nbpoints, self.options['frame_id'])
         std_in+=p_str
         max_weight = self.loadfortran('maxweight',
                                path, std_in)
@@ -3347,7 +3347,7 @@
                 if nb < cut:
                     if key[0]=='full':
                         path=key[1]
-                        end_signal="5 0 0 0 \n"  # before closing, write down the seed 
+                        end_signal="5 0 0 0 0\n"  # before closing, write down the seed 
                         external.stdin.write(end_signal)
                         ranmar_state=external.stdout.readline()
                         ranmar_file=pjoin(path,'ranmar_state.dat')
@@ -3921,7 +3921,8 @@
                 
         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())):
+        if 'init' in self.banner and (eff!=1 or not any(v==-1 for v in self.br_per_id.values())) \
+                                 and not self.options['onlyhelicity']:
             new_init =''
             curr_proc = 0
             has_missing=False
@@ -3961,7 +3962,7 @@
                     external.terminate()
                     del external
                 elif mode=='full':
-                    stdin_text="5 0 0 0 \n"  # before closing, write down the seed 
+                    stdin_text="5 0 0 0 0\n"  # before closing, write down the seed 
                     external = self.calculator[('full',path)]
                     try:
                         external.stdin.write(stdin_text)
@@ -3988,7 +3989,7 @@
             except Exception:
                 pass
             else:
-                stdin_text="5 0 0 0"
+                stdin_text="5 0 0 0 0"
                 external.stdin.write(stdin_text)
                 external.stdin.close()
                 external.stdout.close()

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2018-10-31 10:07:53 +0000
+++ MadSpin/interface_madspin.py	2020-05-23 19:29:41 +0000
@@ -75,7 +75,7 @@
         self.add_param('cross_section', {'__type__':0.}, comment="forcing normalization of cross-section after MS (for none/onshell)" )
         self.add_param('new_wgt', 'cross-section' ,allowed=['cross-section', 'BR'], comment="if not consistent number of particles, choose what to do for the weight. (BR: means local according to number of part, cross use the force cross-section")
         self.add_param('input_format', 'auto', allowed=['auto','lhe', 'hepmc', 'lhe_no_banner'])
-        
+        self.add_param('frame_id', 6)
         
     ############################################################################
     ##  Special post-processing of the options                                ## 
@@ -211,7 +211,7 @@
         # Read the banner of the inputfile
         self.events_file = open(os.path.realpath(inputfile))
         self.banner = banner.Banner(self.events_file)
-        
+
 
         # Check the validity of the banner:
         if 'slha' not in self.banner:
@@ -226,14 +226,21 @@
             raise self.InvalidCmd('This event file was already decayed by MS. This is not possible to add to it a second decay')
         
         if 'mgruncard' in self.banner:
+            run_card = self.banner.charge_card('run_card')
             if not self.options['Nevents_for_max_weigth']:
-                nevents = int(self.banner.get_detail('run_card', 'nevents'))
+                nevents = run_card['nevents']
                 N_weight = max([75, int(3*nevents**(1/3))])
                 self.options['Nevents_for_max_weigth'] = N_weight
                 N_sigma = max(4.5, math.log(nevents,7.7))
                 self.options['nb_sigma'] = N_sigma
             if self.options['BW_cut'] == -1:
                 self.options['BW_cut'] = float(self.banner.get_detail('run_card', 'bwcutoff'))
+            
+            if isinstance(run_card, banner.RunCardLO):
+                run_card.update_system_parameter_for_include()
+                self.options['frame_id'] = run_card['frame_id']
+            else:
+                self.options['frame_id'] = 6
         else:
             if not self.options['Nevents_for_max_weigth']:
                 self.options['Nevents_for_max_weigth'] = 75
@@ -313,7 +320,6 @@
                 if model_name in line:
                     final_model = True
                     
-            
                 
     def import_model(self, args):
         """syntax: import model NAME CARD_PATH
@@ -400,6 +406,12 @@
         #if self.model and not self.model['case_sensitive']:
         #    decaybranch = decaybranch.lower()
 
+        if self.options['spinmode'] != 'full' and '{' in decaybranch:
+            if self.options['spinmode'] == 'none':
+                logger.warning("polarization option used with spinmode=none. The polarization definition will be done according to the rest-frame of the decaying particles (which is likely not what you expect).")
+            else:
+                logger.warning("polarization option used with spinmode=onshell. This combination is not validated and is by construction using sub-optimal method which can likely lead to bias in some situation. Use at your own risk.")
+
         decay_process, init_part = self.decay.reorder_branch(decaybranch)
         if not self.list_branches.has_key(init_part):
             self.list_branches[init_part] = []
@@ -585,14 +597,16 @@
         self.check_launch(args)
         for part in self.list_branches.keys():
             if part in self.mg5cmd._multiparticles:
+            
                 if any(pid in self.final_state for pid in self.mg5cmd._multiparticles[part]):
                     break
-            pid = self.mg5cmd._curr_model.get('name2pdg')[part]
-            if pid in self.final_state:
-                break
-#        else:
-#            logger.info("Nothing to decay ...")
-#            return
+            else:
+                pid = self.mg5cmd._curr_model.get('name2pdg')[part]
+                if pid in self.final_state:
+                    break
+        else:
+            logger.info("Nothing to decay ...")
+            return
         
 
         model_line = self.banner.get('proc_card', 'full_model_line')
@@ -801,6 +815,7 @@
         elif self.options['input_format'] in ['hepmc']:
             import madgraph.various.hepmc_parser as hepmc_parser
             orig_lhe = hepmc_parser.HEPMC_EventFile(filename)
+            orig_lhe.allow_empty_event = True
             logger.info("Parsing input event to know how many decay to generate. This can takes few minuts.")
         else:
             raise Exception
@@ -1067,8 +1082,13 @@
                     if len(hepmc_output) == 0:
                         hepmc_output.append(lhe_parser.Particle(event=hepmc_output))
                         hepmc_output[0].color2 = 0
+                        hepmc_output[0].status = -1
+                        hepmc_output.nexternal+=1
                     decayed_particle = lhe_parser.Particle(particle, hepmc_output)
+                    decayed_particle.mother1 = hepmc_output[0]
+                    decayed_particle.mother2 = hepmc_output[0]
                     hepmc_output.append(decayed_particle)
+                    hepmc_output.nexternal+=1
                     decayed_particle.add_decay(decay)
             # change the weight associate to the event
             if self.options['new_wgt'] == 'cross-section':
@@ -1086,7 +1106,6 @@
             else:
                 hepmc_output.wgt = event.wgt
                 hepmc_output.nexternal = len(hepmc_output) # the append does not update nexternal
-                hepmc_output.assign_mother()
                 output_lhe.write(str(hepmc_output))
         else:
             if counter==0:
@@ -1226,7 +1245,11 @@
                         pass                 
                     self.me_int[decay_dir] = me5_cmd
                 if self.options["run_card"]:
-                    run_card = self.run_card
+                    if hasattr(self, 'run_card'):
+                        run_card = self.run_card
+                    else:
+                        self.run_card = banner.RunCard(self.options["run_card"])
+                        run_card = self.run_card 
                 else:
                     run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))
                 run_card["nevents"] = int(1.2*nb_event)
@@ -1400,6 +1423,7 @@
                     if name not in self.list_branches or len(self.list_branches[name]) == 0:
                         continue
                     raise self.InvalidCmd("The onshell 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)
@@ -1417,7 +1441,6 @@
         #4. determine the maxwgt
         maxwgt = self.get_maxwgt_for_onshell(orig_lhe, evt_decayfile)
         
-        
         #5. generate the decay 
         orig_lhe.seek(0)
         output_lhe = lhe_parser.EventFile(orig_lhe.name.replace('.lhe', '_decayed.lhe'), 'w')
@@ -1435,8 +1458,10 @@
         for curr_event,production in enumerate(orig_lhe):
             if self.options['fixed_order']:
                 production, counterevt= production[0], production[1:]
-            if curr_event and curr_event % 1000 == 0 and float(str(curr_event)[1:]) ==0:
-                print "decaying event number %s. Efficiency: %s [%s s]" % (curr_event, 1/self.efficiency, time.time()-start)
+            if curr_event and self.efficiency and curr_event % 10 == 0 and float(str(curr_event)[1:]) ==0:
+                logger.info("decaying event number %s. Efficiency: %s [%s s]" % (curr_event, 1/self.efficiency, time.time()-start))
+            else:
+                logger.info("next event [%s]", time.time()-start)
             while 1:
                 nb_try +=1
                 decays = self.get_decay_from_file(production, evt_decayfile, nb_event-curr_event)
@@ -1531,6 +1556,9 @@
         # event_decay is a dict pdg -> list of event file (contain the decay)
         
         
+        if self.options['ms_dir'] and os.path.exists(pjoin(self.options['ms_dir'], 'max_wgt')):
+            return float(open(pjoin(self.options['ms_dir'], 'max_wgt'),'r').read())
+        
         nevents = self.options['Nevents_for_max_weigth']
         if nevents == 0 :
             nevents = 75
@@ -1540,7 +1568,6 @@
         logger.info("*****************************")
         logger.info("Probing the first %s events with %s phase space points" % (nevents, self.options['max_weight_ps_point']))
 
-
         self.efficiency = 1. / self.options['max_weight_ps_point']
         start = time.time()
         for i in range(nevents):
@@ -1553,7 +1580,7 @@
                 base_event = base_event[0]
             for j in range(self.options['max_weight_ps_point']):
                 decays = self.get_decay_from_file(base_event, evt_decayfile, nevents-i)   
-                #carefull base_event is modified by the following function                  
+                #carefull base_event is modified by the following function 
                 _, wgt = self.get_onshell_evt_and_wgt(base_event, decays)
                 maxwgt = max(wgt, maxwgt)
             all_maxwgt.append(maxwgt)
@@ -1569,10 +1596,13 @@
             if len(all_maxwgt) < i:
                 break
             ave_weight, std_weight = decay_tools.get_mean_sd(all_maxwgt[:i])
+            #misc.sprint(ave_weight, std_weight)
             base_max_weight = max(base_max_weight, 1.05 * (ave_weight+self.options['nb_sigma']*std_weight))
                 
             if all_maxwgt[1] > base_max_weight:
                 base_max_weight = 1.05 * all_maxwgt[1]
+        if self.options['ms_dir']:
+            open(pjoin(self.options['ms_dir'], 'max_wgt'),'w').write(str(base_max_weight))
         return base_max_weight
             
             
@@ -1594,13 +1624,11 @@
             raise
         import copy
         
-        
         if hasattr(production, 'me_wgt'):
             production_me = production.me_wgt
         else:
             production_me = self.calculate_matrix_element(production)
             production.me_wgt = production_me
-            
         decay_me = 1.0
         for pdg in decays:
             for dec in decays[pdg]:
@@ -1628,7 +1656,6 @@
             final = tuple(-i for i in final)
             tag = (init, final)
             orig_order = self.all_me[tag]['order']
-            
         pdir = self.all_me[tag]['pdir']
         if pdir in self.all_f2py:
             p = event.get_momenta(orig_order)
@@ -1643,11 +1670,15 @@
             
             mymod = __import__("%s.matrix2py" % (pdir))
             reload(mymod)
-            mymod = getattr(mymod, 'matrix2py')  
+            mymod = getattr(mymod, 'matrix2py')
             with misc.chdir(pjoin(self.path_me, 'madspin_me', 'SubProcesses', pdir)):
                 with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                    mymod.initialisemodel(pjoin(self.path_me, 'Cards','param_card.dat'))
-            self.all_f2py[pdir] = mymod.get_value  
+                    if not os.path.exists(pjoin(self.path_me, 'Cards','param_card.dat')) and \
+                            os.path.exists(pjoin(self.path_me,'param_card.dat')):
+                        mymod.initialisemodel(pjoin(self.path_me,'param_card.dat'))
+                    else:
+                        mymod.initialisemodel(pjoin(self.path_me, 'Cards','param_card.dat'))
+            self.all_f2py[pdir] = mymod.get_value 
             return self.calculate_matrix_element(event)
         
         

=== modified file 'MadSpin/src/driver.f'
--- MadSpin/src/driver.f	2018-09-19 09:41:33 +0000
+++ MadSpin/src/driver.f	2020-05-23 19:29:41 +0000
@@ -18,7 +18,7 @@
       INTEGER I,J,K
       INTEGER HELSET(NEXTERNAL)
       REAL*8 P(0:3,NEXTERNAL_PROD) 
-      REAL*8 PFULL(0:3,NEXTERNAL), Ptrial(0:3,NEXTERNAL) 
+      REAL*8 PFULL(0:3,NEXTERNAL), Ptrial(0:3,NEXTERNAL), P2(0:3,NEXTERNAL) 
       double precision x(36), Ecollider
       CHARACTER*120 BUFF(NEXTERNAL_PROD)
       integer iforest(2,-nexternal:-1,N_MAX_CG)
@@ -39,6 +39,7 @@
       double precision mean, variance, maxweight,weight,std
       double precision temp
       double precision Pprod(0:3,nexternal_prod)
+
       integer nb_mc_masses, indices_mc_masses(nexternal)
       double precision values_mc_masses(nexternal)
 
@@ -54,6 +55,9 @@
        double precision BWcut, maxBW
        common /to_BWcut/BWcut, maxBW
 
+       integer frame_id
+       common /to_me_frame/frame_id
+
 c Conflicting BW stuff
       integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
       double precision cBW_mass(-nexternal:-1,-1:1),
@@ -94,7 +98,6 @@
        read(56,*) P_seed
        close(56)
        iseed = iseed + P_seed
-
 cccccccccccccccccccccccccccccccccccccccccccccccccccc
 c   I. read momenta for the production events
 c
@@ -112,7 +115,7 @@
  
 1     continue
       maxBW=0d0
-      read(*,*) mode,  BWcut, Ecollider, temp
+      read(*,*) mode,  BWcut, Ecollider, temp, frame_id
  
 
       if (mode.eq.1) then    ! calculate the maximum weight
@@ -255,13 +258,14 @@
 
            cycle
            endif
-           !do j=1,nexternal
-           !   write (*,*) (pfull(k,j), k=0,3)  
-           !enddo
-           call SMATRIX(pfull,M_full)
-           call SMATRIX_PROD(pprod,M_prod)
-c           write(*,*) 'M_full ', M_full
-c           write(*,*) 'jac',jac
+          
+           call  boost_to_frame(pfull, frame_id, P2)
+           call SMATRIX(P2,M_full)
+
+           call  boost_to_frame_prod(pprod, frame_id,nexternal_prod, P2)
+           call SMATRIX_PROD(P2,M_prod)
+
+
 
            weight=M_full*jac/M_prod
            if (weight.gt.maxweight) then
@@ -335,7 +339,6 @@
            do i = 1, 3*(nexternal-nexternal_prod)+1
               call  ntuple(x(i),0d0,1d0,i,1)
            enddo
-
            call generate_momenta_conf(jac,x,itree,qmass,qwidth,pfull,pprod,map_external2res) 
            if (jac.lt.0d0) then
              counter2=counter2+1 
@@ -366,8 +369,14 @@
 
              cycle
            endif
-           call SMATRIX(pfull,M_full)
-           call SMATRIX_PROD(pprod,M_prod)
+
+           call  boost_to_frame(pfull, frame_id, P2)
+           call SMATRIX(P2,M_full)
+
+
+           call  boost_to_frame_prod(pprod, frame_id,nexternal_prod, P2)
+           call SMATRIX_PROD(P2,M_prod)
+
 
            weight=M_full*jac/M_prod
 
@@ -1892,3 +1901,162 @@
       END
 
 
+
+
+c======================================================================
+C     Subroutine to return momenta in a dedicated frame
+C     frame_id is the tag of the particle to put at rest
+C     frame_id follow the convention of cluster.f (sum 2**(N-1))
+C     -----------------------------------------
+
+      subroutine boost_to_frame(P1, frame_id, P2)
+
+      implicit none
+
+      include 'nexternal.inc'
+
+      DOUBLE PRECISION P1(0:3,NEXTERNAL)
+      DOUBLE PRECISION P2(0:3,NEXTERNAL)
+      DOUBLE PRECISION PBOOST(0:3)
+      integer frame_id
+
+      integer ids(nexternal)
+      integer i,j
+      Pboost(0)= 0d0
+      Pboost(1)= 0d0
+      Pboost(2)= 0d0
+      Pboost(3)= 0d0
+
+c     uncompress
+      call mapid(frame_id, ids)
+c     find the boost momenta --sum of particles--
+      do i=1,nexternal
+       if (ids(i).eq.1)then
+c             write (999,*) (Pboost(j), j=0,3)
+            do j=0,3
+                   Pboost(j) = Pboost(j) + P1(j,i)
+            enddo
+         endif
+      enddo
+      do j=1,3
+          Pboost(j) = -1 * Pboost(j)
+      enddo
+       
+c              write (999,*) (Pboost(j), j=0,3)
+      do i=1, nexternal
+         call boostx(p1(0,i), pboost, p2(0,i))
+      enddo
+      return
+      end
+
+
+
+
+c========================================================================
+
+
+      subroutine mapid(id,ids)
+c**************************************************************************
+c     input:
+c            id        compressed particle id
+c            ids       array of particle ids
+c**************************************************************************
+      implicit none
+      include 'nexternal.inc'
+      integer i, icd, id, ids(nexternal)
+
+      icd=id
+      do i=1,nexternal
+         ids(i)=0
+         if (btest(id,i)) then
+            ids(i)=1
+         endif
+c         write(*,*) 'cluster.f: uncompressed code ',i,' is ',ids(i)
+      enddo
+
+      return
+      end
+
+
+
+
+C     -----------------------------------------
+
+      subroutine boost_to_frame_prod(P1, frame_id, nexternal_prod, P2)
+
+      implicit none
+
+      include 'nexternal.inc'
+
+      DOUBLE PRECISION P1(0:3,*)
+      DOUBLE PRECISION P2(0:3,*)
+      DOUBLE PRECISION PBOOST(0:3)
+      integer frame_id, nexternal_prod
+
+      integer ids(nexternal_prod)
+      integer i,j
+      Pboost(0)= 0d0
+      Pboost(1)= 0d0
+      Pboost(2)= 0d0
+      Pboost(3)= 0d0
+
+c     uncompress
+      call mapid_prod(frame_id, ids,nexternal_prod)
+c     find the boost momenta --sum of particles--
+      do i=1,nexternal_prod
+       if (ids(i).eq.1)then
+c             write (999,*) (Pboost(j), j=0,3)
+            do j=0,3
+                   Pboost(j) = Pboost(j) + P1(j,i)
+            enddo
+         endif
+      enddo
+      do j=1,3
+          Pboost(j) = -1 * Pboost(j)
+      enddo
+
+c              write (999,*) (Pboost(j), j=0,3)
+      do i=1, nexternal_prod
+         call boostx(p1(0,i), pboost, p2(0,i))
+      enddo
+      return
+      end
+
+
+
+
+
+
+
+
+      subroutine mapid_prod(id,ids,nexternal_prod)
+c**************************************************************************
+c     input:
+c            id        compressed particle id
+c            ids       array of particle ids
+c**************************************************************************
+      implicit none
+      include 'nexternal.inc'
+      integer i, icd, id, nexternal_prod, ids(nexternal_prod)
+
+      icd=id
+      do i=1,nexternal_prod
+         ids(i)=0
+         if (btest(id,i)) then
+            ids(i)=1
+         endif
+c         write(*,*) 'cluster.f: uncompressed code ',i,' is ',ids(i)
+      enddo
+
+      return
+      end
+
+
+
+
+
+
+
+
+
+

=== modified file 'Template/Common/Cards/delphes_card_ATLAS.dat'
--- Template/Common/Cards/delphes_card_ATLAS.dat	2017-05-19 20:52:51 +0000
+++ Template/Common/Cards/delphes_card_ATLAS.dat	2020-05-23 19:29:41 +0000
@@ -405,11 +405,10 @@
 # add InputArray InputArray
   add InputArray ECal/ecalTowers
   add InputArray HCal/hcalTowers
+  add InputArray MuonMomentumSmearing/muons
   set OutputArray towers
 }
 
-
-
 ####################
 # Energy flow merger
 ####################

=== modified file 'Template/Common/Cards/delphes_card_CMS.dat'
--- Template/Common/Cards/delphes_card_CMS.dat	2017-05-19 20:52:51 +0000
+++ Template/Common/Cards/delphes_card_CMS.dat	2020-05-23 19:29:41 +0000
@@ -168,7 +168,7 @@
   # set ResolutionFormula {resolution formula as a function of eta and energy}
 
   # resolution formula for electrons
-  # based on arXiv:1405.6569
+  # based on arXiv:1502.02701
   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}

=== modified file 'Template/Common/Cards/delphes_card_default.dat'
--- Template/Common/Cards/delphes_card_default.dat	2017-05-19 20:52:51 +0000
+++ Template/Common/Cards/delphes_card_default.dat	2020-05-23 19:29:41 +0000
@@ -168,7 +168,7 @@
   # set ResolutionFormula {resolution formula as a function of eta and energy}
 
   # resolution formula for electrons
-  # based on arXiv:1405.6569
+  # based on arXiv:1502.02701
   set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
                          (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2018-05-29 13:11:53 +0000
+++ Template/LO/Cards/run_card.dat	2020-05-23 19:29:41 +0000
@@ -14,10 +14,6 @@
 #   To display more options, you can type the command:               *
 #      update full_run_card                                          *
 #*********************************************************************
-#
-#*******************                                                 
-# Running parameters
-#*******************                                                 
 #                                                                    
 #*********************************************************************
 # Tag name for the run (one word)                                    *
@@ -26,7 +22,6 @@
 #*********************************************************************
 # 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.    *
 #*********************************************************************
   %(nevents)s = nevents ! Number of unweighted events requested 
   %(iseed)s   = iseed   ! rnd seed (0=assigned automatically=default))
@@ -40,6 +35,7 @@
      %(ebeam1)s     = ebeam1  ! beam 1 total energy in GeV
      %(ebeam2)s     = ebeam2  ! beam 2 total energy in GeV
 $beam_pol
+$frame
 #*********************************************************************
 # PDF CHOICE: this automatically fixes also alpha_s and its evol.    *
 #*********************************************************************
@@ -61,20 +57,9 @@
 #*********************************************************************
   %(gridpack)s     = gridpack  !True = setting up the grid pack
   %(time_of_flight)s = time_of_flight ! threshold (in mm) below which the invariant livetime is not written (-1 means not written)
-  %(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
- %(asrwgtflavor)s = asrwgtflavor     ! highest quark flavor for a_s reweight
- %(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
+$MLM
+$CKKW
 #*********************************************************************
 #
 #*********************************************************************
@@ -99,10 +84,10 @@
 # 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)
-#*********************************************************************
+#IF(d)# #*********************************************************************
+#IF(d)# # Apply pt/E/eta/dr/mij/kt_durham cuts on decay products or not
+#IF(d)# # (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
+#IF(d)# #*********************************************************************
    %(cut_decays)s  = cut_decays    ! Cut decay products 
 #*********************************************************************
 # Standard Cuts                                                      *
@@ -121,24 +106,12 @@
  %(missetmax)s  = missetmax ! maximum missing Et (sum of neutrino's momenta)
  %(pt_min_pdg)s = pt_min_pdg ! pt cut for other particles (use pdg code). Applied on particle and anti-particle
  %(pt_max_pdg)s	= pt_max_pdg ! pt cut for other particles (syntax e.g. {6: 100, 25: 50}) 
-#*********************************************************************
-# Minimum and maximum E's (in the center of mass frame)              *
-#*********************************************************************
-  %(ej)s  = ej     ! minimum E for the jets 
-  %(eb)s  = eb     ! minimum E for the b 
-  %(ea)s  = ea     ! minimum E for the photons 
-  %(el)s  = el     ! minimum E for the charged leptons 
-  %(ejmax)s   = ejmax ! maximum E for the jets
- %(ebmax)s   = ebmax ! maximum E for the b
- %(eamax)s   = eamax ! maximum E for the photons
- %(elmax)s   = elmax ! maximum E for the charged leptons
- %(e_min_pdg)s = e_min_pdg ! E cut for other particles (use pdg code). Applied on particle and anti-particle
- %(e_max_pdg)s = e_max_pdg ! E cut for other particles (syntax e.g. {6: 100, 25: 50}) 
+$ecut
 #*********************************************************************
 # Maximum and minimum absolute rapidity (for max, -1 means no cut)   *
 #*********************************************************************
-  %(etaj)s = etaj    ! max rap for the jets 
-  %(etab)s  = etab    ! max rap for the b
+ %(etaj)s = etaj    ! max rap for the jets 
+ %(etab)s  = etab    ! max rap for the b
  %(etaa)s  = etaa    ! max rap for the photons 
  %(etal)s  = etal    ! max rap for the charged leptons 
  %(etajmin)s  = etajmin ! min rap for the jets
@@ -172,8 +145,8 @@
  %(dralmax)s  = dralmax ! maxdistance between gamma and lepton
 #*********************************************************************
 # Minimum and maximum invariant mass for pairs                       *
-# WARNING: for four lepton final state mmll cut require to have      *
-#          different lepton masses for each flavor!                  *           
+#IF(llll)# # WARNING: for four lepton final state mmll cut require to have      *
+#IF(llll)# #          different lepton masses for each flavor!                  *           
 #*********************************************************************
  %(mmjj)s   = mmjj    ! min invariant mass of a jet pair 
  %(mmbb)s   = mmbb    ! min invariant mass of a b pair 
@@ -186,14 +159,14 @@
  %(mxx_min_pdg)s = mxx_min_pdg ! min invariant mass of a pair of particles X/X~ (e.g. {6:250})
  %(mxx_only_part_antipart)s = mxx_only_part_antipart ! if True the invariant mass is applied only 
                        ! to pairs of particle/antiparticle and not to pairs of the same pdg codes.  
-#*********************************************************************
-# Minimum and maximum invariant mass for all letpons                 *
-#*********************************************************************
+#IF(LL)# #*********************************************************************
+#IF(LL)# # Minimum and maximum invariant mass for all letpons                 *
+#IF(LL)# #*********************************************************************
  %(mmnl)s   = mmnl    ! min invariant mass for all letpons (l+- and vl) 
  %(mmnlmax)s  = mmnlmax ! max invariant mass for all letpons (l+- and vl) 
-#*********************************************************************
-# Minimum and maximum pt for 4-momenta sum of leptons                *
-#*********************************************************************
+#IF(LL)# #*********************************************************************
+#IF(LL)# # Minimum and maximum pt for 4-momenta sum of leptons                *
+#IF(LL)# #*********************************************************************
  %(ptllmin)s   = ptllmin  ! Minimum pt for 4-momenta sum of leptons(l and vl)
  %(ptllmax)s  = ptllmax  ! Maximum pt for 4-momenta sum of leptons(l and vl)
 #*********************************************************************
@@ -204,9 +177,9 @@
  %(xptb)s  = xptb ! minimum pt for at least one b 
  %(xpta)s  = xpta ! minimum pt for at least one photon 
  %(xptl)s  = xptl ! minimum pt for at least one charged lepton 
-#*********************************************************************
-# Control the pt's of the jets sorted by pt                          *
-#*********************************************************************
+#IF(jj)# #*********************************************************************
+#IF(jj)# # Control the pt's of the jets sorted by pt                          *
+#IF(jj)# #*********************************************************************
  %(ptj1min)s   = ptj1min ! minimum pt for the leading jet in pt
  %(ptj2min)s   = ptj2min ! minimum pt for the second jet in pt
  %(ptj3min)s   = ptj3min ! minimum pt for the third jet in pt
@@ -216,9 +189,9 @@
  %(ptj3max)s  = ptj3max ! maximum pt for the third jet in pt
  %(ptj4max)s  = ptj4max ! maximum pt for the fourth jet in pt
  %(cutuse)s   = cutuse  ! reject event if fails any (0) / all (1) jet pt cuts
-#*********************************************************************
-# Control the pt's of leptons sorted by pt                           *
-#*********************************************************************
+#IF(ll)# #*********************************************************************
+#IF(ll)# # Control the pt's of leptons sorted by pt                           *
+#IF(ll)# #*********************************************************************
  %(ptl1min)s   = ptl1min ! minimum pt for the leading lepton in pt
  %(ptl2min)s   = ptl2min ! minimum pt for the second lepton in pt
  %(ptl3min)s   = ptl3min ! minimum pt for the third lepton in pt
@@ -227,9 +200,9 @@
  %(ptl2max)s  = ptl2max ! maximum pt for the second lepton in pt
  %(ptl3max)s  = ptl3max ! maximum pt for the third lepton in pt
  %(ptl4max)s  = ptl4max ! maximum pt for the fourth lepton in pt
-#*********************************************************************
-# Control the Ht(k)=Sum of k leading jets                            *
-#*********************************************************************
+#IF(JJ)# #*********************************************************************
+#IF(JJ)# # Control the Ht(k)=Sum of k leading jets                            *
+#IF(JJ)# #*********************************************************************
  %(htjmin)s   = htjmin ! minimum jet HT=Sum(jet pt)
  %(htjmax)s  = htjmax ! maximum jet HT=Sum(jet pt)
  %(ihtmin)s   = ihtmin  !inclusive Ht for all partons (including b)
@@ -240,29 +213,21 @@
  %(ht2max)s  = ht2max ! maximum Ht for the two leading jets
  %(ht3max)s  = ht3max ! maximum Ht for the three leading jets
  %(ht4max)s  = ht4max ! maximum Ht for the four leading jets
-#***********************************************************************
-# Photon-isolation cuts, according to hep-ph/9801442                   *
-# When ptgmin=0, all the other parameters are ignored                  *
-# When ptgmin>0, pta and draj are not going to be used                 *
-#***********************************************************************
- %(ptgmin)s = ptgmin ! Min photon transverse momentum
- %(r0gamma)s = R0gamma ! Radius of isolation code
- %(xn)s = xn ! n parameter of eq.(3.4) in hep-ph/9801442
- %(epsgamma)s = epsgamma ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
- %(isoem)s = isoEM ! isolate photons from EM energy (photons and leptons)
-#*********************************************************************
-# WBF cuts                                                           *
-#*********************************************************************
+#IF(aj)# #***********************************************************************
+#IF(aj)# # Photon-isolation cuts, according to hep-ph/9801442                   *
+#IF(aj)# # When ptgmin=0, all the other parameters are ignored                  *
+#IF(aj)# # When ptgmin>0, pta and draj are not going to be used                 *
+#IF(aj)# #***********************************************************************
+#IF(aj)#  %(ptgmin)s = ptgmin ! Min photon transverse momentum
+#IF(aj)#  %(r0gamma)s = R0gamma ! Radius of isolation code
+#IF(aj)#  %(xn)s = xn ! n parameter of eq.(3.4) in hep-ph/9801442
+#IF(aj)#  %(epsgamma)s = epsgamma ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
+#IF(aj)#  %(isoem)s = isoEM ! isolate photons from EM energy (photons and leptons)
+#IF(jj)# #*********************************************************************
+#IF(jj)# # WBF cuts                                                           *
+#IF(jj)# #*********************************************************************
  %(xetamin)s   = xetamin ! minimum rapidity for two jets in the WBF case  
  %(deltaeta)s   = deltaeta ! minimum rapidity for two jets in the WBF case 
-#***********************************************************************
-# 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)                                     *

=== modified file 'Template/LO/Source/.make_opts'
--- Template/LO/Source/.make_opts	2017-08-14 08:20:30 +0000
+++ Template/LO/Source/.make_opts	2020-05-23 19:29:41 +0000
@@ -49,11 +49,15 @@
 ifeq ($(FC), ftn)
 FFLAGS+= -extend-source # for ifort type of compiler
 else
-FFLAGS+= -ffixed-line-length-132
+    VERS="$(shell $(FC) --version | grep ifort -i)"
+    ifeq ($(VERS), "")
+        FFLAGS+= -ffixed-line-length-132
+    else	     
+        FFLAGS+= -extend-source # for ifort type of compiler
+    endif
 endif
 
 
-
 UNAME := $(shell uname -s)
 ifeq ($(origin LDFLAGS), undefined)
 LDFLAGS=$(STDLIB) $(MACFLAG)

=== modified file 'Template/LO/Source/BIAS/ptj_bias/makefile'
--- Template/LO/Source/BIAS/ptj_bias/makefile	2016-08-13 05:32:47 +0000
+++ Template/LO/Source/BIAS/ptj_bias/makefile	2020-05-23 19:29:41 +0000
@@ -19,9 +19,5 @@
 # 'VALID' is the keyword that *must* be returned if everything is in order.
 #
 requirements:	
-ifeq ($(shell $(call CHECK_MG5AMC_VERSION,2.4.2)),True)
 	@echo "VALID"
-else
-	@echo "Error:: MG5aMC is not recent enough (found "$(MG5AMC_VERSION)")"
-	@echo "FAIL"
-endif
+

=== modified file 'Template/LO/Source/PDF/PhotonFlux.f'
--- Template/LO/Source/PDF/PhotonFlux.f	2018-04-20 08:02:12 +0000
+++ Template/LO/Source/PDF/PhotonFlux.f	2020-05-23 19:29:41 +0000
@@ -38,19 +38,31 @@
 
       end
 
-      real*8 function epa_proton(x,q2max)
+      real*8 function epa_proton(x,q2max,beamid)
       integer i
+      integer beamid
       real*8 x,phi_f
       real*8 xin
       real*8 alpha,qz
       real*8 f, qmi,qma, q2max
       real*8 PI
+    
+      integer nb_proton(2), nb_neutron(2)
+      common/to_heavyion_pdg/ nb_proton, nb_neutron
+      double precision mass_ion(2)
+      common/to_heavyion_mass/mass_ion
+
       data PI/3.14159265358979323846/
 
       data xin/0.938/ ! proton mass in GeV
 
       alpha = .0072992701
       qz = 0.71
+
+      if (nb_proton(beamid).ne.1.or.nb_neutron(beamid).ne.0)then
+         xin = mass_ion(beamid)
+         alpha = alpha * nb_proton(beamid)
+      endif
     
 C     // x = omega/E = (E-E')/E
       if (x.lt.1) then

=== modified file 'Template/LO/Source/PDF/pdg2pdf.f'
--- Template/LO/Source/PDF/pdg2pdf.f	2018-03-07 23:11:46 +0000
+++ Template/LO/Source/PDF/pdg2pdf.f	2020-05-23 19:29:41 +0000
@@ -58,12 +58,14 @@
       endif
 
       ipart=ipdg
-      if(iabs(ipart).eq.21) ipart=0
-      if(iabs(ipart).eq.22) ipart=7
-      iporg=ipart
-
+      if(iabs(ipart).eq.21) then
+         ipart=0
+      else if(iabs(ipart).eq.22) then
+         ipart=7
+      else if(iabs(ipart).eq.7) then
+         ipart=7
 c     This will be called for any PDG code, but we only support up to 7
-      if(iabs(ipart).gt.7)then
+      else if(iabs(ipart).gt.7)then
          write(*,*) 'PDF not supported for pdg ',ipdg
          write(*,*) 'For lepton colliders, please set the lpp* '//
      $    'variables to 0 in the run_card'  
@@ -72,6 +74,7 @@
          stop 1
       endif
 
+      iporg=ipart
       ireuse = 0
       do i=1,2
 c     Check if result can be reused since any of last two calls
@@ -131,7 +134,7 @@
          if(ih.eq.3) then       !from the electron
             pdg2pdf=epa_electron(x,q2max)
          elseif(ih .eq. 2) then !from a proton without breaking
-            pdg2pdf=epa_proton(x,q2max)
+            pdg2pdf=epa_proton(x,q2max,beamid)
          endif 
          pdflast(iporg,ireuse)=pdg2pdf
          return

=== modified file 'Template/LO/Source/PDF/pdg2pdf_lhapdf.f'
--- Template/LO/Source/PDF/pdg2pdf_lhapdf.f	2018-03-07 23:11:46 +0000
+++ Template/LO/Source/PDF/pdg2pdf_lhapdf.f	2020-05-23 19:29:41 +0000
@@ -57,6 +57,7 @@
       ipart=ipdg
       if(ipart.eq.21) ipart=0
       if(iabs(ipart).eq.22) ipart=7
+      if(iabs(ipart).eq.7) ipart=7
       iporg=ipart
 
 c     This will be called for any PDG code, but we only support up to 7

=== modified file 'Template/LO/Source/PDF/pdg2pdf_lhapdf6.f'
--- Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2018-03-07 23:11:46 +0000
+++ Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2020-05-23 19:29:41 +0000
@@ -24,7 +24,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(-7:7,20)
-      double precision epa_proton
+      double precision epa_proton, epa_electron
       save ihlast,xlast,xmulast,pdflast,imemlast,ipartlast
       data ihlast/20*-99/
       data ipartlast/20*-99/
@@ -54,12 +54,14 @@
       endif
 
       ipart=ipdg
-      if(iabs(ipart).eq.21) ipart=0
-      if(iabs(ipart).eq.22) ipart=7
-      iporg=ipart
-
+      if(iabs(ipart).eq.21) then
+         ipart=0
+      else if(iabs(ipart).eq.22) then 
+         ipart=7
+      else if (iabs(ipart).eq.7) then
+         ipart=7
 c     This will be called for any PDG code, but we only support up to 7
-      if(iabs(ipart).gt.7)then
+      else if(iabs(ipart).gt.7)then
          write(*,*) 'PDF not supported for pdg ',ipdg
          write(*,*) 'For lepton colliders, please set the lpp* '//
      $    'variables to 0 in the run_card'  
@@ -68,6 +70,7 @@
          stop 1
       endif
 
+      iporg=ipart
 c     Determine the iset used in lhapdf
       call getnset(iset)
       if (iset.ne.1) then
@@ -136,6 +139,8 @@
             pdg2pdf = get_ion_pdf(pdflast(-7, i_replace), ipart, nb_proton(beamid), nb_neutron(beamid))
          endif
          pdg2pdf=pdg2pdf/x
+      else if(ih.eq.3) then       !from the electron
+            pdg2pdf=epa_electron(x,xmu*xmu)
       else if(ih.eq.2) then ! photon from a proton without breaking
           pdg2pdf = epa_proton(x,xmu*xmu)
       else

=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2017-09-24 21:24:55 +0000
+++ Template/LO/Source/dsample.f	2020-05-23 19:29:41 +0000
@@ -316,6 +316,7 @@
                fx =0d0
                wgt=0d0
             endif
+            
             if (nzoom .le. 0) then
                call sample_put_point(wgt,x(1),iter,ipole,itmin) !Store result
             else
@@ -705,6 +706,7 @@
       tmean = 0d0
       trmean = 0d0
       tsigma = 0d0
+      nb_pass_cuts = 0
       kn = 0
       cur_it = 1
       do j=1,ng
@@ -1029,7 +1031,7 @@
 
       integer i,j
 
-      open(26, file='ftn26',status='unknown')
+      open(26, file=name, status='unknown')
       write(26,fmt='(4f21.17)') ((grid(2,i,j),i=1,ng),j=1,invar)
       write(26,*) twgt, force_max_wgt
 c                  write(26,fmt='(4f21.16)') (alpha(i),i=1,maxconfigs)
@@ -1516,7 +1518,7 @@
 c
       include 'genps.inc'
       integer    max_events
-      parameter (max_events=500000) !Maximum # events before get non_zero
+      parameter (max_events=5000000) !Maximum # events before get non_zero
 c
 c     Arguments
 c
@@ -1756,6 +1758,11 @@
          if (kn .ge. max_events .and. non_zero .le. 5) then
             call none_pass(max_events)
          endif
+         if (iteration.eq.1) then
+           if (nb_pass_cuts.ge.1000 .and. non_zero.eq.0) then
+              call none_pass(1000)
+           endif
+         endif
          if (non_zero .ge. events .or. (kn .gt. 200*events .and.
      $        non_zero .gt. 5)) then
 
@@ -2267,8 +2274,20 @@
 c----
 c  Begin Code
 c----
-      write(*,*) 'No points passed cuts!'
-      write(*,*) 'Loosen cuts or increase max_events',max_events
+      if (1000.eq.max_events) then
+         write(*,*) nb_pass_cuts, 
+     &    ' points passed the cut but all returned zero'
+         write(*,*) 'therefore considering this contribution as zero'
+      else if (nb_pass_cuts.gt.0.and.nb_pass_cuts.lt.1000)then
+         write(*,*) 'only', nb_pass_cuts, 
+     &    ' points passed the cut and they all returned zero'
+         write(*,*) 'therefore considering this contribution as zero'
+         write(*,*) 'Loosen cuts or increase max_events if you believe this is not zero'
+      else
+         write(*,*) 'No points passed cuts!'
+         write(*,*) 'Loosen cuts or increase max_events',max_events
+      endif
+
 c      open(unit=22,file=result_file,status='old',access='append',
 c     &           err=23)
 c      write(22,222) 'Iteration',0,'Mean: ',0d0,
@@ -2540,6 +2559,7 @@
 
       write(*,*) "RESET CUMULATIVE VARIABLE"
       non_zero = 0
+      nb_pass_cuts = 0
       do j=1,maxinvar
          do i=1,ng -1
             inon_zero = 0

=== modified file 'Template/LO/Source/genps.inc'
--- Template/LO/Source/genps.inc	2016-04-02 14:23:21 +0000
+++ Template/LO/Source/genps.inc	2020-05-23 19:29:41 +0000
@@ -38,3 +38,8 @@
       double precision hel_jacobian
 	  INTEGER HEL_PICKED
 	  COMMON/HEL_PICKED/HEL_PICKED,hel_jacobian
+c**********************************************************************************
+c       parameters to control is zero-cross section is due to PS or matrix-element
+c**********************************************************************************
+      integer nb_pass_cuts
+      common/cut_efficiency/nb_pass_cuts
\ No newline at end of file

=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2018-02-01 22:30:56 +0000
+++ Template/LO/Source/run.inc	2020-05-23 19:29:41 +0000
@@ -16,6 +16,9 @@
       integer                                        lpp(2)
       double precision    ebeam(2), xbk(2),q2fact(2)
       common/to_collider/ ebeam   , xbk   ,q2fact,   lpp
+
+      integer frame_id
+      common/to_frame_me/frame_id
 c
 c     Number of events
 c
@@ -91,4 +94,7 @@
       double precision mxxmin4pdg(0:25) 
       logical mxxpart_antipart(1:25)
       common/TO_PDG_SPECIFIC_CUT/pdg_cut, ptmin4pdg,ptmax4pdg, Emin4pdg, Emax4pdg, etamin4pdg,
-     &etamax4pdg, mxxmin4pdg,mxxpart_antipart
\ No newline at end of file
+     &etamax4pdg, mxxmin4pdg,mxxpart_antipart
+
+       double precision small_width_treatment
+       common/narrow_width/small_width_treatment
\ No newline at end of file

=== modified file 'Template/LO/Source/transpole.f'
--- Template/LO/Source/transpole.f	2016-06-15 22:46:42 +0000
+++ Template/LO/Source/transpole.f	2020-05-23 19:29:41 +0000
@@ -28,13 +28,24 @@
       double precision z,zmin,zmax,xmin,xmax,ez
       double precision pole1,width1,x,xc
       double precision a,b
+c
+c     small width treatment
+c
+      double precision small_width_treatment
+      common/narrow_width/small_width_treatment
 c-----
 c  Begin Code
 c-----
       pole=pole1
       width=width1
+
       x = x1
       if (pole .gt. 0d0) then
+         if (width.lt.pole*small_width_treatment)then
+            width = pole * small_width_treatment
+            jac = jac * width/width1
+         endif
+
          zmin = atan((-pole)/width)/width
          zmax = atan((1d0-pole)/width)/width
          if (x .gt. del .and. x .lt. 1d0-del) then
@@ -175,7 +186,11 @@
 c
       double precision pole1,width1,y1,jac
       real*8 x
-
+c
+c     small width treatment
+c
+      double precision small_width_treatment
+      common/narrow_width/small_width_treatment
 c
 c     Local
 c
@@ -191,6 +206,10 @@
       width=width1
       y = y1
       if (pole .gt. 0d0) then                   !BW 
+         if (width.lt.pole*small_width_treatment)then
+            width = pole * small_width_treatment
+            jac = jac * width/width1
+         endif
          zmin = atan((-pole)/width)/width
          zmax = atan((1d0-pole)/width)/width
          z = atan((y-pole)/width)/width

=== modified file 'Template/LO/SubProcesses/cuts.f'
--- Template/LO/SubProcesses/cuts.f	2018-01-11 16:16:25 +0000
+++ Template/LO/SubProcesses/cuts.f	2020-05-23 19:29:41 +0000
@@ -1592,6 +1592,7 @@
       m0d = 0.0
       m0c = 1.5
       m0s = 0.0
+      m0b = 4.7 
       m0t = 172.5
       m0w = 80.4
       m0z = 91.188

=== modified file 'Template/LO/SubProcesses/dummy_fct.f'
--- Template/LO/SubProcesses/dummy_fct.f	2018-01-11 16:16:25 +0000
+++ Template/LO/SubProcesses/dummy_fct.f	2020-05-23 19:29:41 +0000
@@ -86,3 +86,13 @@
       shat = x(1)*x(2)*ebeam(1)*ebeam(2)
       return 
       end
+
+
+      logical  function dummy_boostframe()
+      implicit none
+c
+c      
+      dummy_boostframe = .false.
+      return
+      end
+      

=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2018-06-27 09:03:22 +0000
+++ Template/LO/SubProcesses/genps.f	2020-05-23 19:29:41 +0000
@@ -185,7 +185,8 @@
 c     External function
       double precision SumDot
       external SumDot
-
+      logical dummy_boostframe
+      external dummy_boostframe
 c
 c     data
 c
@@ -304,8 +305,12 @@
 c     First Generate Momentum for initial state particles
 c
       if (lpp(1).eq.9.or.lpp(2).eq.9)then
-         p(:,1) = pi1(:)
-         p(:,2) = pi2(:)
+         if (dummy_boostframe())then
+            call mom2cx(m(-nbranch),m(1),m(2),1d0,0d0,p(0,1),p(0,2))
+         else
+            p(:,1) = pi1(:)
+            p(:,2) = pi2(:)
+         endif
       else if(nincoming.eq.2) then
         call mom2cx(m(-nbranch),m(1),m(2),1d0,0d0,p(0,1),p(0,2))
       else
@@ -840,7 +845,7 @@
          m12 = m(itree(2,ibranch))**2
          mn2 = m(ibranch-1)**2
 c         write(*,*) 'Enertering yminmax',sqrt(s1),sqrt(m12),sqrt(mn2)
-         call yminmax(s1,t,m12,ma2,mb2,mn2,tmin,tmax)
+         call yminmax(s1,0d0,m12,ma2,mb2,mn2,tmin,tmax)
 c
 c     Call for 0<x<1
 c
@@ -1114,6 +1119,7 @@
 C     This is the G function from Particle Kinematics by
 C     E. Byckling and K. Kajantie, Chapter 4 p. 91 eqs 5.28
 C     It is used to determine physical limits for Y based on inputs
+C     Y is not used in this formula (called with dummy value)
 C**************************************************************************
       implicit none
 c
@@ -1258,5 +1264,47 @@
       END
       
 
+C     -----------------------------------------
+C     Subroutine to return momenta in a dedicated frame
+C     frame_id is the tag of the particle to put at rest
+C     frame_id follow the convention of cluster.f (sum 2**(N-1))
+C     -----------------------------------------
+
+      subroutine boost_to_frame(P1, frame_id, P2)
+
+      implicit none
+
+      include 'nexternal.inc'
+
+      DOUBLE PRECISION P1(0:3,NEXTERNAL)
+      DOUBLE PRECISION P2(0:3,NEXTERNAL)
+      DOUBLE PRECISION PBOOST(0:3)
+      integer frame_id
+
+      integer ids(nexternal)
+      integer i,j
+
+c     uncompress
+      call mapid(frame_id, ids)
+      pboost(:) = 0d0
+      p2(:,:) = 0d0
+c     find the boost momenta --sum of particles--
+      do i=1,nexternal
+       if (ids(i).eq.1)then
+            do j=0,3
+	           Pboost(j) = Pboost(j) + P1(j,i)
+            enddo
+         endif
+      enddo
+      do j=1,3	
+          Pboost(j) = -1 * Pboost(j)
+      enddo	    
+      do i=1, nexternal
+         call boostx(p1(0,i), pboost, p2(0,i))
+      enddo   
+      return
+      end
+
+
 
 

=== modified file 'Template/LO/SubProcesses/myamp.f'
--- Template/LO/SubProcesses/myamp.f	2016-12-17 16:26:56 +0000
+++ Template/LO/SubProcesses/myamp.f	2020-05-23 19:29:41 +0000
@@ -1,91 +1,3 @@
-      double precision function testamp(p)
-c*****************************************************************************
-c     Approximates matrix element by propagators
-c*****************************************************************************
-      implicit none
-c
-c     Constants
-c     
-      include 'genps.inc'
-      include 'maxconfigs.inc'
-      include 'nexternal.inc'
-      double precision   zero
-      parameter (zero = 0d0)
-c
-c     Arguments
-c
-      double precision p(0:3,nexternal)
-c      integer iconfig
-c
-c     Local
-c
-      double precision xp(0:3,-nexternal:nexternal)
-      double precision mpole(-nexternal:0),shat,tsgn
-      integer i,j,iconfig
-
-      double precision prmass(-nexternal:0,lmaxconfigs)
-      double precision prwidth(-nexternal:0,lmaxconfigs)
-      integer pow(-nexternal:0,lmaxconfigs)
-      logical first_time
-c
-c     Global
-c
-      integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
-      integer            mapconfig(0:lmaxconfigs), this_config
-      common/to_mconfigs/mapconfig, this_config
-      
-      include 'coupl.inc'
-c
-c     External
-c
-      double precision dot
-
-      save prmass,prwidth,pow
-      data first_time /.true./
-c-----
-c  Begin Code
-c-----      
-      iconfig = this_config
-      if (first_time) then
-c         include 'props.inc'
-         first_time=.false.
-      endif
-
-      do i=1,nexternal
-         mpole(-i)=0d0
-         do j=0,3
-            xp(j,i)=p(j,i)
-         enddo
-      enddo
-c      mpole(-3) = 174**2
-c      shat = dot(p(0,1),p(0,2))/(1800)**2
-      shat = dot(p(0,1),p(0,2))/(10)**2
-c      shat = 1d0
-      testamp = 1d0
-      tsgn    = +1d0
-      do i=-1,-(nexternal-3),-1              !Find all the propagotors
-         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
-         do j=0,3
-            xp(j,i) = xp(j,iforest(1,i,iconfig))
-     $           +tsgn*xp(j,iforest(2,i,iconfig))
-         enddo
-         if (prwidth(i,iconfig) .ne. 0d0 .and. .false.) then
-            testamp=testamp/((dot(xp(0,i),xp(0,i))
-     $                        -prmass(i,iconfig)**2)**2
-     $         -(prmass(i,iconfig)*prwidth(i,iconfig))**2)
-         else
-            testamp = testamp/((dot(xp(0,i),xp(0,i)) -
-     $                          prmass(i,iconfig)**2)
-     $                          **(pow(i,iconfig)))
-         endif
-        testamp=testamp*shat**(pow(i,iconfig))
-c        write(*,*) i,iconfig,pow(i,iconfig),prmass(i,iconfig)
-      enddo
-c      testamp = 1d0/dot(xp(0,-1),xp(0,-1))
-      testamp=abs(testamp)
-c      testamp = 1d0
-      end
 
       logical function cut_bw(p)
 c*****************************************************************************
@@ -114,6 +26,7 @@
 
       double precision prmass(-nexternal:0,lmaxconfigs)
       double precision prwidth(-nexternal:0,lmaxconfigs)
+      double precision prwidth_tmp(-nexternal:0,lmaxconfigs)
       integer pow(-nexternal:0,lmaxconfigs)
       logical first_time, onshell
       double precision xmass
@@ -152,13 +65,14 @@
 c
       double precision dot
 
-      save prmass,prwidth,pow
+      save prmass,prwidth,pow,prwidth_tmp
       data first_time /.true./
 c-----
 c  Begin Code
 c-----      
       cut_bw = .false.    !Default is we passed the cut
       iconfig = this_config
+
       if (first_time) then
          include 'props.inc'
          nbw = 0
@@ -209,10 +123,12 @@
 c            write(*,*) 'xmass',xmass,prmass(i,iconfig)
 c
 c           Here we set if the BW is "on-shell" for LesHouches
-c
+c            
+            prwidth_tmp(i,iconfig) = max(prwidth(i,iconfig), prmass(i,iconfig)*small_width_treatment)
+ 
             onshell = (abs(xmass - prmass(i,iconfig)) .lt.
-     $           bwcutoff*prwidth(i,iconfig).and.
-     $           (prwidth(i,iconfig)/prmass(i,iconfig).lt.0.1d0.or.
+     $           bwcutoff*prwidth_tmp(i,iconfig).and.
+     $           (prwidth_tmp(i,iconfig)/prmass(i,iconfig).lt.0.1d0.or.
      $            gForceBW(i,iconfig).eq.1))
             if(onshell)then
 c     Remove on-shell forbidden s-channels (gForceBW=2) (JA 2/10/11)
@@ -264,10 +180,10 @@
 c     For decay-chain syntax use BWcutoff here too (22/12/14)
             if (gForceBW(i, iconfig).eq.1) then
                onshell = (abs(xmass - prmass(i,iconfig)) .lt.
-     $           bwcutoff*prwidth(i,iconfig))
+     $           bwcutoff*prwidth_tmp(i,iconfig))
             else
                onshell = (abs(xmass - prmass(i,iconfig)) .lt.
-     $           5d0*prwidth(i,iconfig))
+     $           5d0*prwidth_tmp(i,iconfig))
             endif
 
             if (onshell .and. (lbw(nbw).eq. 2) .or.
@@ -313,6 +229,7 @@
 
       double precision prmass(-nexternal:0,lmaxconfigs)
       double precision prwidth(-nexternal:0,lmaxconfigs)
+      double precision prwidth_tmp(-nexternal:0,lmaxconfigs)
       integer pow(-nexternal:0,lmaxconfigs)
 
       integer idup(nexternal,maxproc,maxsproc)
@@ -382,11 +299,21 @@
 
 c-----
 c  Begin Code
-c-----      
+c-----     
+      iconfig = this_config
+c     needs to be initialise to avoid segfault
+      do i = -nexternal,-1
+         prwidth(i,iconfig) = 0
+         prmass(i,iconfig) =0
+      enddo 
       include 'props.inc'
 c      etmin = 10
       nt = 0
-      iconfig = this_config
+      do i = -nexternal,-1
+         prwidth_tmp(i,iconfig) = max(prwidth(i,iconfig), prmass(i,iconfig)*small_width_treatment)
+      enddo
+
+
       mtot = 0d0
       etot = 0d0   !Total energy needed
       spmass = 0d0 !Keep track of BW masses for shat
@@ -423,7 +350,7 @@
 c     Look for identical particles to map radiation processes
       call idenparts(iden_part, iforest(1,-max_branch,iconfig),
      $     sprop(1,-max_branch,iconfig), gForceBW(-max_branch,iconfig),
-     $     prwidth(-nexternal,iconfig))
+     $     prwidth_tmp(-nexternal,iconfig))
 
 c     Start loop over propagators
       do i=-1,-(nexternal-3),-1
@@ -446,15 +373,15 @@
               xm(i)=max(xm(i),max(xqcutij(l1,l2),0d0))
             endif
 c            write(*,*) 'iconfig,i',iconfig,i
-c            write(*,*) prwidth(i,iconfig),prmass(i,iconfig)
-            if (prwidth(i,iconfig) .gt. 0 ) then
+c            write(*,*) prwidth_tmp(i,iconfig),prmass(i,iconfig)
+            if (prwidth_tmp(i,iconfig) .gt. 0 ) then
                nbw=nbw+1
 c              JA 6/8/2011 Set xe(i) for resonances
                if (gforcebw(i,iconfig).eq.1) then
-                  xm(i) = max(xm(i), prmass(i,iconfig)-bwcutoff*prwidth(i,iconfig))
+                  xm(i) = max(xm(i), prmass(i,iconfig)-bwcutoff*prwidth_tmp(i,iconfig))
                   bwcut_for_PS(i) = bwcutoff
                else if (lbw(nbw).eq.1) then
-                  xm(i) = max(xm(i), prmass(i,iconfig)-5d0*prwidth(i,iconfig))
+                  xm(i) = max(xm(i), prmass(i,iconfig)-5d0*prwidth_tmp(i,iconfig))
                   bwcut_for_PS(i) = 5d0
                else
                   bwcut_for_PS(i) = 5d0
@@ -465,19 +392,19 @@
 c     Either: required onshell and daughter masses too large
 c     Or: forced and daughter masses too large
 c     Or: required offshell and forced
-            if(prwidth(i,iconfig) .gt. 0.and.
+            if(prwidth_tmp(i,iconfig) .gt. 0.and.
      $         (lbw(nbw).eq.1.and.
-     $          (prmass(i,iconfig)+bwcut_for_PS(i)*prwidth(i,iconfig).lt.xm(i)
-     $           .or.prmass(i,iconfig)-bwcut_for_PS(i)*prwidth(i,iconfig).gt.dsqrt(stot))
+     $          (prmass(i,iconfig)+bwcut_for_PS(i)*prwidth_tmp(i,iconfig).lt.xm(i)
+     $           .or.prmass(i,iconfig)-bwcut_for_PS(i)*prwidth_tmp(i,iconfig).gt.dsqrt(stot))
      $          .or.gforcebw(i,iconfig).eq.1.and.
-     $              prmass(i,iconfig)+bwcutoff*prwidth(i,iconfig).lt.xm(i)
+     $              prmass(i,iconfig)+bwcutoff*prwidth_tmp(i,iconfig).lt.xm(i)
      $          .or.lbw(nbw).eq.2.and.gforcebw(i,iconfig).eq.1))
      $        then
 c     Write results.dat and quit
                call write_null_results()
                stop
             endif
-            if (prwidth(i,iconfig) .gt. 0 .and. lbw(nbw) .le. 1) then         !B.W.
+            if (prwidth_tmp(i,iconfig) .gt. 0 .and. lbw(nbw) .le. 1) then         !B.W.
                if (i .eq. -(nexternal-(nincoming+1))) then  !This is s-hat
                   j = 3*(nexternal-2)-4+1    !set i to ndim+1
 c-----
@@ -489,27 +416,27 @@
      $                 .or. lbw(nbw).eq.1) then
                      write(*,*) 'Setting PDF BW',j,nbw,prmass(i,iconfig)
                      spole(j)=prmass(i,iconfig)*prmass(i,iconfig)/stot
-                     swidth(j) = prwidth(i,iconfig)*prmass(i,iconfig)/stot
+                     swidth(j) = prwidth(i,iconfig)*prmass(i,iconfig)/stot ! keep the real width here (important for the jacobian)
                   endif
-               else if((prmass(i,iconfig)+bwcut_for_PS(i)*prwidth(i,iconfig)).ge.xm(i)
+               else if((prmass(i,iconfig)+bwcut_for_PS(i)*prwidth_tmp(i,iconfig)).ge.xm(i)
      $                  .and. iden_part(i).eq.0 .or. lbw(nbw).eq.1) then
 c              JA 02/13 Only allow BW if xm below M+5*Gamma
                   write(*,*) 'Setting BW',i,nbw,prmass(i,iconfig)
                   spole(-i)=prmass(i,iconfig)*prmass(i,iconfig)/stot
-                  swidth(-i) = prwidth(i,iconfig)*prmass(i,iconfig)/stot
+                  swidth(-i) = prwidth(i,iconfig)*prmass(i,iconfig)/stot ! keep the real width here (important for the jacobian)
                endif
 c     JA 4/1/2011 Set grid in case there is no BW (radiation process)
                if (swidth(-i) .eq. 0d0 .and.
      $              i.ne.-(nexternal-(nincoming+1)))then
                   a=prmass(i,iconfig)**2/stot
                   xo = min(xm(i)**2/stot, 1-1d-8)
-                  if (xo.eq.0d0) xo=1d0/stot
+                  if (xo.eq.0d0) xo=MIN(10d0/stot, stot/50d0, 0.5)
                   call setgrid(-i,xo,a,1)
                endif
 c     Set spmass for BWs
                if (swidth(-i) .ne. 0d0)
      $              spmass=spmass-xm(i) +
-     $              max(xm(i),prmass(i,iconfig)-bwcut_for_PS(i)*prwidth(i,iconfig))
+     $              max(xm(i),prmass(i,iconfig)-bwcut_for_PS(i)*prwidth_tmp(i,iconfig))
             else                                  !1/x^pow
               a=prmass(i,iconfig)**2/stot
 c     JA 4/1/2011 always set grid
@@ -529,8 +456,8 @@
                     xo = (MMJJ * 0.8)**2/stot
                  endif
               endif
-              if (xo.eq.0d0) xo=1d0/stot
-c              if (prwidth(i, iconfig) .eq. 0d0.or.iden_part(i).gt.0) then 
+              if (xo.eq.0d0) xo=MIN(10d0/stot, stot/50d0, 0.5)
+c              if (prwidth_tmp(i, iconfig) .eq. 0d0.or.iden_part(i).gt.0) then 
               call setgrid(-i,xo,a,1)
 c              else 
 c                 write(*,*) 'Using flat grid for BW',i,nbw,
@@ -574,7 +501,7 @@
 c           Use 1/10000 of sqrt(s) as minimum, to always get integration
             xo = xo*xo/stot
             if (xo.eq.0d0)then
-               xo=1d0/stot
+               xo=1/10000d0
                write(*,*) 'Warning: No cutoff for shat integral found'
                write(*,*) '         Minimum set to ', xo
             endif
@@ -587,7 +514,7 @@
          endif
       enddo
 c     Perform setting for shat (PDF BW or 1/s)
-      if (abs(lpp(1)) .eq. 1 .or. abs(lpp(2)) .eq. 1) then
+      if (abs(lpp(1)) .ge. 1 .or. abs(lpp(2)) .ge. 1) then
 c     Set minimum based on: 1) required energy 2) resonances 3) 1/10000 of sqrt(s)
          i = max(1,3*(nexternal-2) - 4 + 1)
          xo = max(min(etot**2/stot, 1d0-1d-8),1d0/stot)
@@ -610,7 +537,7 @@
             spole(i)= -2.0d0    ! 1/s pole
             write(*,*) "Transforming s_hat 1/s ",i,xo, smin, stot
          else
-            write(*,*) "Transforming s_hat BW ",spole(i),swidth(i)
+            write(*,*) "Transforming s_hat BW ",spole(i), max(swidth(i), spole(i)*small_width_treatment)
          endif
       endif
 

=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2018-07-22 11:40:04 +0000
+++ Template/LO/SubProcesses/reweight.f	2020-05-23 19:29:41 +0000
@@ -316,7 +316,7 @@
           ipart(1,imo)=ipart(1,ida2)
           ipart(2,imo)=ipart(2,ida2)
         endif
-      else if(idmo.eq.21 .and. abs(idda1).le.6 .and.
+      else if(idmo.eq.21.and. abs(idda1).le.6 .and.
      $        abs(idda2).le.6) then
 c     gluon -> quark anti-quark: use both, but take hardest as 1
         if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
@@ -327,6 +327,17 @@
           ipart(1,imo)=ipart(1,ida2)
           ipart(2,imo)=ipart(1,ida1)
         endif
+      else if (get_color(idmo).eq.8.and.iabs(get_color(idda1)).eq.3.and.
+     $         iabs(get_color(idda2)).eq.3) then
+c     gluon' -> quark' anti-quark': use both, but take hardest as 1
+        if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
+     $     p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
+          ipart(1,imo)=ipart(1,ida1)
+          ipart(2,imo)=ipart(1,ida2)
+        else
+          ipart(1,imo)=ipart(1,ida2)
+          ipart(2,imo)=ipart(1,ida1)
+        endif
       else if(idmo.eq.21.and.(idda1.eq.21.or.idda2.eq.21))then
          if(idda1.eq.21) then
             iddgluon = idda1
@@ -380,10 +391,32 @@
 c       exotic q > Scalar q'
         ipart(1,imo)=ipart(1,ida2)
         ipart(2,imo)=0
-      else if (get_color(idmo).eq.1) then
+      else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda2)).eq.8.and.iabs(get_color(idda1)).eq.3) then
+c       exotic q >  q' gluon
+        ipart(1,imo)=ipart(1,ida2)
+        ipart(2,imo)=0
+      else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda2)).eq.3.and.get_color(idda1).eq.8) then
+c       exotic q > gluon q'
+        ipart(1,imo)=ipart(1,ida1)
+        ipart(2,imo)=0
+      else if (get_color(idmo).eq.1.or.get_color(idmo).eq.2.or.get_color(idda1).eq.2.or.get_color(idda2).eq.2) then
 c     Color singlet
          ipart(1,imo)=ipart(1,ida1)
          ipart(2,imo)=ipart(1,ida2)
+      else if (abs(get_color(idmo)).eq.3.and.abs(get_color(idda1)).eq.3.and.abs(get_color(idda2)).eq.3)then
+c     EPSILON^IJK color structure -- handle it as photon emission --
+         ipart(1,imo)=ipart(1,ida1)
+         ipart(2,imo)=0
+      else if (abs(get_color(idmo)).eq.6.and.abs(get_color(idda1)).eq.3.and.abs(get_color(idda2)).eq.3)then
+c     sextet -> (anti-)quark (anti-)quark': use both, but take hardest as 1
+        if(p(1,ipart(1,ida1))**2+p(2,ipart(1,ida1))**2.gt.
+     $     p(1,ipart(1,ida2))**2+p(2,ipart(1,ida2))**2) then
+          ipart(1,imo)=ipart(1,ida1)
+          ipart(2,imo)=ipart(1,ida2)
+        else
+          ipart(1,imo)=ipart(1,ida2)
+          ipart(2,imo)=ipart(1,ida1)
+        endif
       else
          write(*,*) idmo,'>', idda1, idda2, 'color', get_color(idmo),'>', get_color(idda1), get_color(idda2)
          write(*,*) "failed for ipartupdate. Please retry without MLM/default dynamical scale"
@@ -641,12 +674,12 @@
 C   If we have fixed factorization scale, for ickkw>0 means central
 C   scale, i.e. last two scales (ren. scale for these vertices are
 C   anyway already set by "scale" above)
-      if(ickkw.gt.0) then
-         if(fixed_fac_scale.and.first)then
+      if (first) then
             q2bck(1)=q2fact(1)
             q2bck(2)=q2fact(2)
             first=.false.
-         else if(fixed_fac_scale) then
+      else if(ickkw.gt.0) then
+         if(fixed_fac_scale) then
             q2fact(1)=q2bck(1)
             q2fact(2)=q2bck(2)
          endif
@@ -775,15 +808,17 @@
 c          Check QCD jet, take care so not a decay
            if(.not.isjetvx(imocl(n),idacl(n,1),idacl(n,2),
      $        ipdgcl(1,igraphs(1),iproc),ipart,n.eq.nexternal-2)) then
+              pdgm   = ipdgcl(imocl(n),igraphs(1),iproc)
+              pdgid1 = ipdgcl(idacl(n,1),igraphs(1),iproc)
+              pdgid2 = ipdgcl(idacl(n,2),igraphs(1),iproc)
+              if (isqcd(pdgm).and.isqcd(pdgid1).and.isqcd(pdgid2))then
+                  continue
 c          Remove non-gluon jets that lead up to non-jet vertices
-           if(ipart(1,imocl(n)).gt.2)then ! ipart(1) set and not IS line
+           elseif(ipart(1,imocl(n)).gt.2)then ! ipart(1) set and not IS line
 c          The ishft gives the FS particle corresponding to imocl
               if(.not.is_octet(ipdgcl(ishft(1,ipart(1,imocl(n))-1),igraphs(1),iproc)))then
                  ! split case for q a > q and for g > g h (with the gluon splitting into quark)
                  ! also check for case of three scalar interaction (then do nothing)
-                 pdgm   = ipdgcl(imocl(n),igraphs(1),iproc)
-                 pdgid1 = ipdgcl(idacl(n,1),igraphs(1),iproc)
-                 pdgid2 = ipdgcl(idacl(n,2),igraphs(1),iproc)
                  if (.not.isqcd(pdgm).and..not.isqcd(pdgid1).and..not.isqcd(pdgid2)) then
                     ! this is to avoid to do weird stuff for w+ w- z (or h h h) 
                     ! this fix an issue for qq_zttxqq G1594.08
@@ -955,7 +990,7 @@
      $              '. Have jets (>0)',(iqjets(i),i=1,nexternal),
      $              ', should be ',
      $              (iqjetstore(i,iconfig),i=1,njetstore(iconfig))
-               stop
+               stop 4
             endif
             if (btest(mlevel,3))
      $           write(*,*) 'Bad clustering, jets fail. Reclustering ',
@@ -1697,6 +1732,10 @@
          if (btest(mlevel,3))
      $        write(*,*)' set fact scales for PS to ',
      $        sqrt(q2fact(1)),sqrt(q2fact(2))
+      else if (abs(lpp(1)).eq.2.or.abs(lpp(1)).eq.3) then
+         q2fact(1)=q2bck(1)
+      else if (abs(lpp(2)).eq.2.or.abs(lpp(2)).eq.3) then
+         q2fact(2)=q2bck(2)
       endif
 
       if (btest(mlevel,3)) then

=== modified file 'Template/LO/SubProcesses/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2018-03-11 12:10:19 +0000
+++ Template/LO/SubProcesses/setcuts.f	2020-05-23 19:29:41 +0000
@@ -664,14 +664,21 @@
          endif
       enddo
 
-c     ensure symmetry of s_min(i,j)
+c     ensure symmetry of s_min(i,j) + mass effect
       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))
+               s_min(i,j) = max(s_min(j,i), s_min(i,j), (pmass(i)+pmass(j))**2)
             endif
         enddo
       enddo
+
+c     ensure that smin is bigger than mass threshold
+      smin_p = 0d0
+      do i=nincoming+1, nexternal
+         smin_p = smin_p + pmass(i)
+      enddo
+      smin = max(smin, smin_p**2)
       write(*,*) "Define smin to", smin
 
 c Check that results are consistent among all the grouped subprocesses

=== modified file 'Template/LO/bin/internal/Gridpack/compile'
--- Template/LO/bin/internal/Gridpack/compile	2015-09-28 20:34:44 +0000
+++ Template/LO/bin/internal/Gridpack/compile	2020-05-23 19:29:41 +0000
@@ -63,7 +63,7 @@
 
 if [[ -d Source ]]; then
     cd Source
-    for i in ../bin/internal/gen_ximprove all ../bin/internal/combine_events; do
+    for i in ../bin/internal/gen_ximprove all; do
     make $i -j 1 > /dev/null
     if [[ $? -ne 0 ]];then
         # Make didn't exit successfully

=== modified file 'Template/NLO/Cards/FKS_params.dat'
--- Template/NLO/Cards/FKS_params.dat	2017-08-04 19:52:15 +0000
+++ Template/NLO/Cards/FKS_params.dat	2020-05-23 19:29:41 +0000
@@ -153,6 +153,16 @@
 0.005d0
 ! Default :: 0.005d0
 !
+! Use a polynomial fit for the virtual corrections. Gets a better accuracy,
+! but uses more memmory and disk space.
+#UsePolyVirtual
+.False.
+! Default :: .True.
+!
+! ==========================================================================
+! Combination of matrix elements
+! ==========================================================================
+!
 ! For fixed order calculations only (ignored for (N)LO+PS runs): This
 ! parameter determines if parton flavour configurations with idential
 ! matrix elements should be separated at the level of the analysis or
@@ -170,11 +180,5 @@
 ! Default :: .false.
 !
 ! ==========================================================================
-! Arbitrary numerical parameters used in the FKS formalism
-! ==========================================================================
-!
-! To be implemented by the FKS authors
-!
-! ==========================================================================
 ! End of FKS_params.dat file
 ! ==========================================================================

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIGPP.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIGPP.Script	2016-10-26 13:53:02 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIGPP.Script	2020-05-23 19:29:41 +0000
@@ -149,12 +149,15 @@
 if [ -f $HWPPPATH/bin/Herwig ]; then HWEXE=Herwig ; fi
 HWVER=`$HWPPPATH/bin/$HWEXE --version`
 HWVER=($(echo "$HWVER" | sed 's/ / /g'))
-HWVER=($(echo ${HWVER[1]} | tr "." "\n"))
-HWVER=${HWVER[0]}
-if [ "$HWVER" -ge 7 ]; then HWNEW=1 ; HWEXE=Herwig   ; fi
-if [ "$HWVER" -le 2 ]; then HWNEW=0 ; HWEXE=Herwig++ ; fi
+HWVER1=($(echo ${HWVER[1]} | tr "." "\n"))
+HWVER2=($(echo ${HWVER1[1]} | tr "." "\n"))
+HWVER=$(echo "$HWVER1$HWVER2" | bc)
+if [ "$HWVER" -ge 72 ]; then HWNEW=3 ; HWEXE=Herwig   ; fi
+if [ "$HWVER" -eq 71 ]; then HWNEW=2 ; HWEXE=Herwig   ; fi
+if [ "$HWVER" -eq 70 ]; then HWNEW=1 ; HWEXE=Herwig   ; fi
+if [ "$HWVER" -le 69 ]; then HWNEW=0 ; HWEXE=Herwig++ ; fi
 
-if [[ "$PDFCODE" -ne 0 && "$HWNEW" -eq 0 ]] || [ "$HWNEW" -eq 1 ]
+if [[ "$PDFCODE" -ne 0 && "$HWNEW" -eq 0 ]] || [ "$HWNEW" -ge 1 ]
 then
   HERPDF='EXTPDF'
   PDFLIBRARY='LHAPDF'
@@ -166,7 +169,7 @@
   if [ "$PDFCODE" -gt 1 ]; then PDFSET=$PDFCODE; fi
   if [ "$PDFCODE" -lt 0 ]; then LHAOFL=EXTRAPOLATE; fi
   if [ "$PDFCODE" -lt -1 ]; then PDFSET=$((PDFCODE*-1)); fi
-  if [[ "$PDFCODE" -eq 0 && "$HWNEW" -eq 1 ]]; then PDFSET=25000 ; LHAOFL=FREEZE ; fi
+  if [[ "$PDFCODE" -eq 0 && "$HWNEW" -ge 1 ]]; then PDFSET=25000 ; LHAOFL=FREEZE ; fi
   linklhapdf
 # convert PDF number to name reading PDFsets.index
   ISLHA5=0
@@ -184,7 +187,7 @@
       echo 'Unable to understand LHAPDF version!'
       exit 1
   fi
-  if [[ "$ISLHA5" -eq 1 && "$HWNEW" -eq 1 ]]
+  if [[ "$ISLHA5" -eq 1 && "$HWNEW" -ge 1 ]]
   then
       echo " Please, use LHAPDF6 with HW7 "
       exit 1
@@ -282,7 +285,7 @@
   HWPPUTIMAKE="HWPPUTI="$HWPPUTI
   FORFILEMAKE="FORFILE="$HWPPUTI
   HEPMCF=HepMCFortran7.h
-  if [ "$HWVER" -le 2 ]; then HEPMCF=HepMCFortran2.h; fi
+  if [ "$HWVER" -le 69 ]; then HEPMCF=HepMCFortran2.h; fi
   cp $thisdir/HWPPAnalyzer/$HEPMCF $thisdir/HWPPAnalyzer/HepMCFortran.h
   HWPPANALYZER=HepMCFortran
   if ! make -C $thisdir/HWPPAnalyzer $HWPPUTI $HWPPANALYZER.so $HWPPMAKE "$LIBSMAKE" "$LIBSMAKEP" "$INCLMAKE" "$FORFILEMAKE"
@@ -385,6 +388,16 @@
 mv $EXEC_NAMEFF $exefile
 
 # create input file
+if [ "$HWNEW" -le 1 ]
+then
+GENERATOR="LHCGenerator"
+EXTRACTOR="QCDExtractor"
+elif [ "$HWNEW" -ge 2 ]
+then
+GENERATOR="EventGenerator"
+EXTRACTOR="PPExtractor"
+fi
+
 cat <<EOF > ./$ifile
 #
 #  Create the Les Houches handler for MC@NLO
@@ -394,35 +407,47 @@
 
 # --------------------------------------------------------------------- #
 # 1.) NECESSARY SETTINGS FOR RUNNING WITH MC@NLO EVENTS (DO NOT MODIFY)
+EOF
 
+if [ "$HWNEW" -le 1 ]
+then
+cat <<EOF >> ./$ifile
 set /Herwig/Shower/Evolver:HardVetoMode 1
 set /Herwig/Shower/Evolver:HardVetoScaleSource 1
 set /Herwig/Shower/Evolver:MECorrMode 0
-
-
+EOF
+fi
+
+cat <<EOF >> ./$ifile
 # create the Handler & Reader
 create ThePEG::LesHouchesFileReader /Herwig/EventHandlers/LHEReader
 create ThePEG::LesHouchesEventHandler /Herwig/EventHandlers/LHEHandler 
 insert /Herwig/EventHandlers/LHEHandler:LesHouchesReaders 0 /Herwig/EventHandlers/LHEReader
-
+EOF
+
+
+cat <<EOF >> ./$ifile
 set /Herwig/EventHandlers/LHEReader:AllowedToReOpen 0
 set /Herwig/EventHandlers/LHEReader:MomentumTreatment RescaleEnergy
 set /Herwig/EventHandlers/LHEReader:WeightWarnings 0
-
+EOF
+
+
+cat <<EOF >> ./$ifile
 set /Herwig/EventHandlers/LHEHandler:WeightOption VarNegWeight
-set /Herwig/EventHandlers/LHEHandler:PartonExtractor /Herwig/Partons/QCDExtractor
+set /Herwig/EventHandlers/LHEHandler:PartonExtractor /Herwig/Partons/$EXTRACTOR
 set /Herwig/EventHandlers/LHEHandler:CascadeHandler /Herwig/Shower/ShowerHandler 
 set /Herwig/EventHandlers/LHEHandler:HadronizationHandler /Herwig/Hadronization/ClusterHadHandler
 set /Herwig/EventHandlers/LHEHandler:DecayHandler /Herwig/Decays/DecayHandler
-set /Herwig/Generators/LHCGenerator:EventHandler /Herwig/EventHandlers/LHEHandler
-set /Herwig/Generators/LHCGenerator:NumberOfEvents $NEVENTS
-set /Herwig/Generators/LHCGenerator:RandomNumberGenerator:Seed $RNDEVSEED_HWPP
-set /Herwig/Generators/LHCGenerator:PrintEvent $MAXPR_HWPP
-set /Herwig/Generators/LHCGenerator:MaxErrors $ERR_NUM_HWPP 
-set /Herwig/Generators/LHCGenerator:DebugLevel 1
+set /Herwig/Generators/$GENERATOR:EventHandler /Herwig/EventHandlers/LHEHandler
+set /Herwig/Generators/$GENERATOR:NumberOfEvents $NEVENTS
+set /Herwig/Generators/$GENERATOR:RandomNumberGenerator:Seed $RNDEVSEED_HWPP
+set /Herwig/Generators/$GENERATOR:PrintEvent $MAXPR_HWPP
+set /Herwig/Generators/$GENERATOR:MaxErrors $ERR_NUM_HWPP 
+set /Herwig/Generators/$GENERATOR:DebugLevel 1
 EOF
 
-if [ $HWNEW = 0 ]
+if [ "$HWNEW" -eq 0 ]
 then
 cat <<EOF >> ./$ifile
 insert /Herwig/EventHandlers/LHEHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler
@@ -430,7 +455,7 @@
 fi
 
 whichherpdf $HERPDF
-if [[ "$pdftype" -eq 1 && "$HWNEW" -eq 0 ]] || [ "$HWNEW" -eq 1 ]
+if [[ "$pdftype" -eq 1 && "$HWNEW" -eq 0 ]] || [ "$HWNEW" -ge 1 ]
 then
 cat <<EOF >> ./$ifile
 # Define PDF from MCatNLO.inputs PDFSET number
@@ -438,6 +463,16 @@
 create ThePEG::LHAPDF /LHAPDF/MCNLOPDF ThePEGLHAPDF.so
 set /LHAPDF/MCNLOPDF:PDFName $PDFSETN
 set /LHAPDF/MCNLOPDF:RemnantHandler /Herwig/Partons/HadronRemnants
+EOF
+
+if [ "$HWNEW" -ge 1 ]
+then
+cat <<EOF >> ./$ifile
+set /Herwig/Partons/RemnantDecayer:AllowTop Yes
+EOF
+fi
+
+cat <<EOF >> ./$ifile
 set /Herwig/EventHandlers/LHEReader:PDFA /LHAPDF/MCNLOPDF
 set /Herwig/EventHandlers/LHEReader:PDFB /LHAPDF/MCNLOPDF
 set /Herwig/Particles/p+:PDF /LHAPDF/MCNLOPDF
@@ -494,9 +529,9 @@
 # 3.) ADDITIONAL SETTINGS
 
 # SHOWER SETTINGS
-set /Herwig/Shower/Evolver:IntrinsicPtGaussian 2.2*GeV
 EOF
 
+
 if ((`bc <<< "$LAMBDAHERW >= 0.0"`))
 then
 LAMBDAHERW=$(echo "$LAMBDAHERW * 1000" | bc)
@@ -512,23 +547,43 @@
 EOF
 fi
 
+if [ "$HWNEW" -le 1 ]
+then
+cat <<EOF >> ./$ifile
+set /Herwig/Shower/Evolver:IntrinsicPtGaussian 2.2*GeV
+EOF
+elif [ "$HWNEW" -ge 2 ]
+then
+cat <<EOF >> ./$ifile
+set /Herwig/Shower/ShowerHandler:IntrinsicPtGaussian 1.3*GeV
+EOF
+fi
 
-if [ $HWNEW = 1 ]
+if [ "$HWNEW" -ge 1 ]
 then
 cat <<EOF >> ./$ifile
 set /Herwig/Shower/PartnerFinder:PartnerMethod Random
 set /Herwig/Shower/PartnerFinder:ScaleChoice Partner
+EOF
+fi
+
+if [ "$HWNEW" -le 1 ]
+then
+cat <<EOF >> ./$ifile
 set /Herwig/Shower/GtoQQbarSplitFn:AngularOrdered Yes
-
-#set /Herwig/Shower/Evolver:Interactions QCDOnly
-#set /Herwig/Shower/Evolver:SoftCorrelations No
-#set /Herwig/Shower/Evolver:SpinCorrelations Yes
-
-set /Herwig/Shower/Evolver:SpinCorrelations No
-set /Herwig/Shower/KinematicsReconstructor:FinalStateReconOption Default
-set /Herwig/Shower/KinematicsReconstructor:InitialStateReconOption Rapidity
+set /Herwig/Shower/GammatoQQbarSplitFn:AngularOrdered Yes
+#set /Herwig/Shower/ShowerHandler:Interactions QCDOnly
+#set /Herwig/Shower/ShowerHandler:SoftCorrelations No
+#set /Herwig/Shower/ShowerHandler:SpinCorrelations Yes
+EOF
+fi
+
+if [ "$HWNEW" -ge 2 ]
+then
+cat <<EOF >> ./$ifile
+set /Herwig/Shower/ShowerHandler:SpinCorrelations No
 set /Herwig/Shower/ShowerHandler:MaxPtIsMuF Yes
-set /Herwig/Shower/ShowerHandler:RestrictPhasespace On
+set /Herwig/Shower/ShowerHandler:RestrictPhasespace Yes
 
 EOF
 fi
@@ -807,36 +862,53 @@
 # 4.) HW++ SETTINGS THAT ARE NOT ALLOWED TO BE TOUCHED BY THE USER
 
 #Boost and reconstruction stuff
+EOF
+
+if [ "$HWNEW" -eq 1 ] || [ "$HWNEW" -eq 2 ]
+then
+cat <<EOF >> ./$ifile
+set /Herwig/Shower/KinematicsReconstructor:FinalStateReconOption Default
+set /Herwig/Shower/KinematicsReconstructor:InitialStateReconOption Rapidity
+EOF
+fi
+
+if [ "$HWNEW" -eq 3 ]
+then
+cat <<EOF >> ./$ifile
+set /Herwig/Shower/KinematicsReconstructor:InitialStateReconOption Rapidity
+EOF
+fi
+
+cat <<EOF >> ./$ifile
 set /Herwig/Shower/KinematicsReconstructor:ReconstructionOption General
 set /Herwig/Shower/KinematicsReconstructor:InitialInitialBoostOption LongTransBoost
 
 # --------------------------------------------------------------------- #
 # 5.) CREATE THE ANALYZER AND SAVE THE RUN (DO NOT MODIFY)
-
 EOF
 
 if [ "$HWPPUTI" != "" ]
 then
 cat <<EOF >> ./$ifile
 create MCatNLO::$HWPPANALYZER $HWPPANALYZER $HWPPANALYZER.so
-insert /Herwig/Generators/LHCGenerator:AnalysisHandlers 0 $HWPPANALYZER
+insert /Herwig/Generators/$GENERATOR:AnalysisHandlers 0 $HWPPANALYZER
 
 EOF
 else
 cat <<EOF >> ./$ifile
 #create MCatNLO::hepfortr hepfortr hepfortr.so
-#insert /Herwig/Generators/LHCGenerator:AnalysisHandlers 0 hepfortr
+#insert /Herwig/Generators/$GENERATOR:AnalysisHandlers 0 hepfortr
 
 #  Useful analysis handlers for HepMC related output
 ##################################################
 # Schematic overview of an event (requires --with-hepmc to be set at configure time
 # and the graphviz program 'dot' to produce a plot)
 #create MCatNLO::$HWPPANALYZER $HWPPANALYZER $HWPPANALYZER.so
-#insert /Herwig/Generators/LHCGenerator:AnalysisHandlers 0 $HWPPANALYZER
-# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/Plot
+#insert /Herwig/Generators/$GENERATOR:AnalysisHandlers 0 $HWPPANALYZER
+# insert $GENERATOR:AnalysisHandlers 0 /Herwig/Analysis/Plot
 # A HepMC dump file (requires --with-hepmc to be set at configure time)
-# insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
-insert /Herwig/Generators/LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
+# insert $GENERATOR:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
+insert /Herwig/Generators/$GENERATOR:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
 set /Herwig/Analysis/HepMCFile:PrintEvent $NEVENTS
 set /Herwig/Analysis/HepMCFile:Format GenEvent
 set /Herwig/Analysis/HepMCFile:Units GeV_mm
@@ -845,9 +917,10 @@
 fi
 
 cat <<EOF >> ./$ifile
-saverun $runfile /Herwig/Generators/LHCGenerator
+saverun $runfile /Herwig/Generators/$GENERATOR
 EOF
 
+
 ( cd . ;\
      $exefile )
 

=== modified file 'Template/NLO/MCatNLO/shower_template.sh'
--- Template/NLO/MCatNLO/shower_template.sh	2017-01-16 14:11:24 +0000
+++ Template/NLO/MCatNLO/shower_template.sh	2020-05-23 19:29:41 +0000
@@ -9,8 +9,11 @@
 
 export %(ld_library_path)s=$%(ld_library_path)s:%(extralibs)s
 
-# this is for py8 runs
-export PYTHIA8DATA=`pwd`/xmldoc
+if [ $SHOWER == "HERWIGPP" ] ; then
+    export PYTHIA8DATA=""
+else
+    export PYTHIA8DATA=`pwd`/xmldoc
+fi
 
 # if one is splitting file cd to a new dir and link all files here
 if [[ "$NFILE" != "" ]]; then

=== modified file 'Template/NLO/Source/PDF/pdg2pdf.f'
--- Template/NLO/Source/PDF/pdg2pdf.f	2015-11-06 16:54:54 +0000
+++ Template/NLO/Source/PDF/pdg2pdf.f	2020-05-23 19:29:41 +0000
@@ -66,12 +66,14 @@
       endif
 
       ipart=ipdg
-      if(iabs(ipart).eq.21) ipart=0
-      if(iabs(ipart).eq.22) ipart=7
-      iporg=ipart
-
+      if(iabs(ipart).eq.21) then
+         ipart=0
+      else if(iabs(ipart).eq.22) then
+         ipart=7
+      else if(iabs(ipart).eq.7) then
+         ipart=7
+      else if(iabs(ipart).gt.7)then
 c     This will be called for any PDG code, but we only support up to 7
-      if(iabs(ipart).gt.7)then
 C         write(*,*) 'PDF not supported for pdg ',ipdg
 C         write(*,*) 'For lepton colliders, please set the lpp* '//
 C     $    'variables to 0 in the run_card'  
@@ -82,6 +84,7 @@
          return
       endif
 
+      iporg=ipart
       ireuse = 0
       ii=i_replace
       do i=1,20

=== modified file 'Template/NLO/Source/PDF/pdg2pdf_lhapdf.f'
--- Template/NLO/Source/PDF/pdg2pdf_lhapdf.f	2016-05-24 12:48:05 +0000
+++ Template/NLO/Source/PDF/pdg2pdf_lhapdf.f	2020-05-23 19:29:41 +0000
@@ -64,12 +64,14 @@
       endif
 
       ipart=ipdg
-      if(iabs(ipart).eq.21) ipart=0
-      if(iabs(ipart).eq.22) ipart=7
-      iporg=ipart
-
+      if(iabs(ipart).eq.21) then
+         ipart=0
+      else if(iabs(ipart).eq.22) then
+         ipart=7
+      else if(iabs(ipart).eq.7) then
+         ipart=7
+      else if(iabs(ipart).gt.7)then
 c     This will be called for any PDG code, but we only support up to 7
-      if(iabs(ipart).gt.7)then
 C         write(*,*) 'PDF not supported for pdg ',ipdg
 C         write(*,*) 'For lepton colliders, please set the lpp* '//
 C     $    'variables to 0 in the run_card'  
@@ -86,6 +88,7 @@
 c     Determine the member of the set (function of lhapdf)
       call getnmem(iset,imem)
 
+      iporg=ipart
       ireuse = 0
       ii=i_replace
       do i=1,20

=== modified file 'Template/NLO/Source/PDF/pdg2pdf_lhapdf6.f'
--- Template/NLO/Source/PDF/pdg2pdf_lhapdf6.f	2017-01-10 15:57:32 +0000
+++ Template/NLO/Source/PDF/pdg2pdf_lhapdf6.f	2020-05-23 19:29:41 +0000
@@ -65,12 +65,14 @@
       endif
 
       ipart=ipdg
-      if(iabs(ipart).eq.21) ipart=0
-      if(iabs(ipart).eq.22) ipart=7
-      iporg=ipart
-
+      if(iabs(ipart).eq.21) then
+         ipart=0
+      else if(iabs(ipart).eq.22) then
+         ipart=7
+      else if(iabs(ipart).eq.7) then
+         ipart=7
+      else if(iabs(ipart).gt.7)then
 c     This will be called for any PDG code, but we only support up to 7
-      if(iabs(ipart).gt.7)then
 C         write(*,*) 'PDF not supported for pdg ',ipdg
 C         write(*,*) 'For lepton colliders, please set the lpp* '//
 C     $    'variables to 0 in the run_card'  
@@ -86,6 +88,7 @@
 c     Determine the member of the set (function of lhapdf)
       call getnmem(iset,imem)
 
+      iporg=ipart
       ireuse = 0
       ii=i_replace
       do i=1,20

=== modified file 'Template/NLO/Source/kin_functions.f'
--- Template/NLO/Source/kin_functions.f	2017-07-14 08:23:18 +0000
+++ Template/NLO/Source/kin_functions.f	2020-05-23 19:29:41 +0000
@@ -236,7 +236,7 @@
 c-----
 c  Begin Code
 c-----
-      if ((lpp(1).eq.0).and.(lpp(2).eq.0)) then
+      if ((lpp(1).ne.1).and.(lpp(2).ne.1)) then
          p1a = dsqrt(p1(1)**2+p1(2)**2+p1(3)**2)
          p2a = dsqrt(p2(1)**2+p2(2)**2+p2(3)**2)
          if (p1a*p2a .ne. 0d0) then

=== modified file 'Template/NLO/Source/makefile'
--- Template/NLO/Source/makefile	2017-07-31 08:48:20 +0000
+++ Template/NLO/Source/makefile	2020-05-23 19:29:41 +0000
@@ -13,6 +13,7 @@
 	  zerox64_cernlib.o
 
 .f.o: ; $(FC) $(FFLAGS) -c $*.f
+.f90.o: ; $(FC) $(FFLAGS) -c $*.f90
 
 all: $(LIBDIR)libdhelas.a $(LIBDIR)libgeneric.a $(LIBDIR)libpdf.a	\
     $(LIBDIR)libmodel.a $(LIBDIR)libcernlib.a param_card.inc

=== modified file 'Template/NLO/Source/ranmar.f'
--- Template/NLO/Source/ranmar.f	2017-01-30 15:18:32 +0000
+++ Template/NLO/Source/ranmar.f	2020-05-23 19:29:41 +0000
@@ -1,18 +1,3 @@
-      function ran2()
-c Wrapper for the random numbers; needed for the NLO stuff
-      implicit none
-      include '../SubProcesses/mint.inc' ! includes iconfig common
-      double precision ran2,x,a,b
-      integer ii,jconfig
-      a=0d0           ! min allowed value for x
-      b=1d0           ! max allowed value for x
-      ii=0            ! dummy argument of ntuple
-      jconfig=iconfig ! integration channel (for off-set)
-      call ntuple(x,a,b,ii,jconfig)
-      ran2=x
-      return
-      end
-
       subroutine ntuple(x,a,b,ii,jconfig)
 c-------------------------------------------------------
 c     Front to ranmar which allows user to easily

=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
--- Template/NLO/SubProcesses/BinothLHA.f	2017-06-23 14:14:03 +0000
+++ Template/NLO/SubProcesses/BinothLHA.f	2020-05-23 19:29:41 +0000
@@ -4,12 +4,11 @@
 c that calls the OLP and returns the virtual weights. For convenience
 c also the born_wgt is passed to this subroutine.
 c
+      use FKSParams
       implicit none
       include "nexternal.inc"
       include "coupl.inc"
       include 'born_nhel.inc'
-c general MadFKS parameters
-      include 'FKSParams.inc'
       double precision pi, zero,mone
       parameter (pi=3.1415926535897932385d0)
       parameter (zero=0d0)
@@ -72,6 +71,9 @@
       integer ret_code_common
       common /to_ret_code/ret_code_common
       double precision born_hel_from_virt
+
+      logical updateloop
+      common /to_updateloop/updateloop
 c masses
       include 'pmass.inc'
       data nbad / 0 /
@@ -80,7 +82,11 @@
 c update the ren_scale for MadLoop and the couplings (should be the
 c Ellis-Sexton scale)
       mu_r = sqrt(QES2)
+      ! force to update also loop-related parameters
+      updateloop=.true.
       call update_as_param()
+      updateloop=.false.
+
       alpha_S=g**2/(4d0*PI)
       ao2pi= alpha_S/(2d0*PI)
       virt_wgt= 0d0
@@ -111,8 +117,13 @@
 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.)
+         IF (.not. force_polecheck) THEN ! still have the pole for the pole check
+            CALL COLLIER_COMPUTE_UV_POLES(.FALSE.)
+            CALL COLLIER_COMPUTE_IR_POLES(.FALSE.)
+         else
+            CALL COLLIER_COMPUTE_UV_POLES(.TRUE.)
+            CALL COLLIER_COMPUTE_IR_POLES(.TRUE.)
+         endif
          firsttime_run = .false.
       endif
       firsttime=firsttime.or.force_polecheck
@@ -290,7 +301,7 @@
             if ((dabs(avgPoleRes(1))+dabs(avgPoleRes(2))).ne.0d0) then
                cpol = .not.((((PoleDiff(1)+PoleDiff(2))/
      $              (dabs(avgPoleRes(1))+dabs(avgPoleRes(2)))) .lt.
-     $              tolerance*10d0).or.(mod(ret_code,10).eq.7))
+     $              tolerance*10d0).or.(mod(ret_code,10).eq.7.and..not.force_polecheck))
             else
                cpol = .not.((PoleDiff(1)+PoleDiff(2).lt.tolerance*10d0)
      $              .or.(mod(ret_code,10).eq.7))

=== modified file 'Template/NLO/SubProcesses/BinothLHA_OLP.f'
--- Template/NLO/SubProcesses/BinothLHA_OLP.f	2015-08-13 09:26:37 +0000
+++ Template/NLO/SubProcesses/BinothLHA_OLP.f	2020-05-23 19:29:41 +0000
@@ -6,11 +6,11 @@
 c
 C************************************************************************
 c
+      use FKSParams
       implicit none
       include "nexternal.inc"
       include "coupl.inc"
       include "Binoth_proc.inc"
-      include "FKSParams.inc"
       double precision pi
       parameter (pi=3.1415926535897932385d0)
       double precision pin(0:3,nexternal-1),p(0:4,nexternal-1)

=== removed file 'Template/NLO/SubProcesses/FKSParamReader.f'
=== added file 'Template/NLO/SubProcesses/FKSParams.f90'
--- Template/NLO/SubProcesses/FKSParams.f90	1970-01-01 00:00:00 +0000
+++ Template/NLO/SubProcesses/FKSParams.f90	2020-05-23 19:29:41 +0000
@@ -0,0 +1,248 @@
+!====================================================================
+!
+!  Define common block with all general parameters used by MadFKS 
+!  See their definitions in the file FKS_params.dat.
+!
+!====================================================================
+module FKSParams
+  character(len=64), parameter ::  paramFileName='FKS_params.dat'
+  integer,parameter :: maxContribsSelected=100, &
+                       maxCouplingsSelected=100, &
+                       maxContribType=15, &
+                       maxCouplingTypes=20
+  real*8 :: IRPoleCheckThreshold,Virt_fraction, PrecisionVirtualAtRunTime,Min_virt_fraction
+  integer  :: NHelForMCoverHels,VetoedContributionTypes(0:maxContribsSelected), &
+              SelectedContributionTypes(0:maxContribsSelected),QED_squared_selected, &
+              SelectedCouplingOrders(maxCouplingTypes,0:maxCouplingsSelected), &
+              QCD_squared_selected
+  logical :: separate_flavour_configs,IncludeBornContributions,use_poly_virtual
+
+contains
+
+  subroutine FKSParamReader(filename, printParam, force)
+    ! Reads the file 'filename' and sets the parameters found in that file.
+    implicit none
+    logical, save :: HasReadOnce=.False.,paramPrinted=.false.
+    logical :: force,couldRead,printParam
+    character(*) :: filename
+    CHARACTER(len=64) :: buff, buff2, mode
+    include "orders.inc"
+    integer :: i,j
+    couldRead=.False.
+    if (HasReadOnce.and..not.force) then
+       goto 901
+    endif
+! Make sure to have default parameters if not set in the FKSParams.dat card
+! (if it is an old one for instance)
+    call DefaultFKSParam()
+! Overwrite the default parameters from file:
+    open(68, file=fileName, err=676, action='READ')
+    do
+       read(68,*,end=999) buff
+       if(index(buff,'#').eq.1) then
+          if (buff .eq. '#IRPoleCheckThreshold') then
+             read(68,*,end=999) IRPoleCheckThreshold
+             if (IRPoleCheckThreshold .lt. -1.01d0 ) then
+                stop 'IRPoleCheckThreshold must be >= -1.0d0.'
+             endif
+          elseif (buff .eq. '#PrecisionVirtualAtRunTime') then
+             read(68,*,end=999) PrecisionVirtualAtRunTime
+             if (IRPoleCheckThreshold .lt. -1.01d0 ) then
+                stop 'PrecisionVirtualAtRunTime must be >= -1.0d0.'
+             endif
+          else if (buff .eq. '#NHelForMCoverHels') then
+             read(68,*,end=999) NHelForMCoverHels
+             if (NHelForMCoverHels .lt. -1) then
+                stop 'NHelForMCoverHels must be >= -1.'
+             endif
+          else if (buff .eq. '#QCD^2==') then
+             read(68,*,end=999) QCD_squared_selected
+             if (QCD_squared_selected .lt. -1) then
+                stop 'QCD_squared_selected must be >= -1.'
+             endif
+          else if (buff .eq. '#QED^2==') then
+             read(68,*,end=999) QED_squared_selected
+             if (QED_squared_selected .lt. -1) then
+                stop 'QED_squared_selected must be >= -1.'
+             endif
+          else if (buff .eq. '#VirtualFraction') then
+             read(68,*,end=999) Virt_fraction
+             if (Virt_fraction .lt. 0 .or. virt_fraction .gt.1) then
+                stop 'VirtualFraction should be a fraction between 0 and 1'
+             endif
+          else if (buff .eq. '#MinVirtualFraction') then
+             read(68,*,end=999) Min_Virt_fraction
+             if (min_virt_fraction .lt. 0 .or. min_virt_fraction .gt.1) then
+                stop 'VirtualFraction should be a fraction between 0 and 1'
+             endif
+          else if (buff .eq. '#SeparateFlavourConfigurations') then
+             read(68,*,end=999) separate_flavour_configs
+          else if (buff .eq. '#UsePolyVirtual') then
+             read(68,*,end=999) use_poly_virtual
+          else if (buff .eq. '#VetoedContributionTypes') then
+             read(68,*,end=999) VetoedContributionTypes(0)
+             if (VetoedContributionTypes(0) .lt. 0.or. &
+                  VetoedContributionTypes(0) .gt. maxContribsSelected) then
+                write(*,*) 'VetoedContributionTypes length should be >= 0 and <=', &
+                     maxContribsSelected
+                stop 'Format error in FKS_params.dat.'
+             endif
+             read(68,*,end=999) (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
+             do I=1,VetoedContributionTypes(0)
+                if (VetoedContributionTypes(I).lt.1.or. &
+                     VetoedContributionTypes(I).gt.maxContribType) then
+                   write(*,*) 'VetoedContributionTypes must be >=1 and <=',maxContribType
+                   stop 'Format error in FKS_params.dat.'
+                endif
+             enddo
+             do I=VetoedContributionTypes(0)+1,maxContribsSelected
+                VetoedContributionTypes(I)=-1
+             enddo
+
+          else if (buff .eq. '#SelectedContributionTypes') then
+             read(68,*,end=999) SelectedContributionTypes(0)
+             if (SelectedContributionTypes(0) .lt. 0 .or. &
+                  SelectedContributionTypes(0) .gt. maxContribsSelected) then
+                write(*,*) 'SelectedContributionTypes length should be >= 0 and <=', &
+                     maxContribsSelected
+                stop 'Format error in FKS_params.dat.'
+             endif
+             read(68,*,end=999) (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
+             do I=1,SelectedContributionTypes(0)
+                if (SelectedContributionTypes(I).lt.1.or. &
+                     SelectedContributionTypes(I).gt.maxContribType) then
+                   write(*,*) 'SelectedContributionTypes must be >=1 and <=',maxContribType
+                   stop 'Format error in FKS_params.dat.'
+                endif
+             enddo
+             do I=SelectedContributionTypes(0)+1,maxContribsSelected
+                SelectedContributionTypes(I)=-1
+             enddo
+          else if (buff .eq. '#SelectedCouplingOrders') then
+             read(68,*,end=999) SelectedCouplingOrders(1,0)
+             if (SelectedCouplingOrders(1,0) .lt. 0 .or. &
+                  SelectedCouplingOrders(1,0) .gt. maxCouplingsSelected) then
+                write(*,*) 'SelectedCouplingOrders length should be >= 0 and <=', &
+                     maxCouplingsSelected
+                stop 'Format error in FKS_params.dat.'                
+             endif
+             do j = 2, maxCouplingTypes
+                SelectedCouplingOrders(j,0) = SelectedCouplingOrders(1,0)
+             enddo
+             do j = 1, SelectedCouplingOrders(1,0)
+                read(68,*,end=999) (SelectedCouplingOrders(i,j),i=1,nsplitorders)
+                do i=nsplitorders+1,maxCouplingTypes
+                   SelectedCouplingOrders(i,j)=-1
+                enddo
+             enddo
+          else
+             write(*,*) 'The parameter name ',buff(2:),'is not reckognized.'
+             stop 'Format error in FKS_params.dat.'
+          endif
+       endif
+    enddo
+999 continue
+    couldRead=.True.
+    goto 998      
+
+676 continue
+    write(*,*) 'ERROR :: MadFKS parameter file ',fileName, &
+         ' could not be found or is malformed. Please specify it.'
+    stop 
+    !   Below is the code if one desires to let the code continue with
+    !   a non existing or malformed parameter file
+    write(*,*) 'WARNING :: The file ',fileName,' could not be ', &
+         ' open or did not contain the necessary information. The ', &
+         ' default MadFKS parameters will be used.'
+    call DefaultFKSParam()
+    goto 998
+998 continue
+
+    if(printParam.and..not.paramPrinted) then
+       write(*,*) &
+            '==============================================================='
+       if (couldRead) then      
+          write(*,*) 'INFO: MadFKS read these parameters from ',filename
+       else
+          write(*,*) 'INFO: MadFKS used the default parameters.'
+       endif
+       write(*,*) &
+            '==============================================================='
+       write(*,*) ' > IRPoleCheckThreshold      = ',IRPoleCheckThreshold
+       write(*,*) ' > PrecisionVirtualAtRunTime = ',PrecisionVirtualAtRunTime
+       if (SelectedContributionTypes(0).gt.0) then
+          write(*,*) ' > SelectedContributionTypes = ', &
+               (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
+       else
+          write(*,*) ' > SelectedContributionTypes = All'
+       endif
+       if (VetoedContributionTypes(0).gt.0) then
+          write(*,*) ' > VetoedContributionTypes   = ', &
+               (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
+       else
+          write(*,*) ' > VetoedContributionTypes   = None'
+       endif
+       if (QCD_squared_selected.eq.-1) then
+          write(*,*) ' > QCD_squared_selected      = All'
+       else
+          write(*,*) ' > QCD_squared_selected      = ',QCD_squared_selected
+       endif
+       if (QED_squared_selected.eq.-1) then
+          write(*,*) ' > QED_squared_selected      = All'
+       else
+          write(*,*) ' > QED_squared_selected      = ',QED_squared_selected
+       endif
+       if (SelectedCouplingOrders(1,0).gt.0) then
+          do j=1,SelectedCouplingOrders(1,0)
+             write(*,*) ' > SelectedCouplingOrders(',j,') = ', &
+                  (SelectedCouplingOrders(i,j),i=1,nsplitorders)
+          enddo
+       else
+          write(*,*) ' > SelectedCouplingOrders    = All'
+       endif
+       write(*,*) ' > NHelForMCoverHels         = ',NHelForMCoverHels
+       write(*,*) ' > VirtualFraction           = ',Virt_fraction
+       write(*,*) ' > MinVirtualFraction        = ',Min_virt_fraction
+       write(*,*) ' > SeparateFlavourConfigs    = ',separate_flavour_configs
+       write(*,*) ' > UsePolyVirtual            = ',use_poly_virtual
+       write(*,*) &
+            '==============================================================='
+       paramPrinted=.TRUE.
+    endif
+
+    close(68)
+    HasReadOnce=.TRUE.
+901 continue
+  end subroutine FKSParamReader
+
+  subroutine DefaultFKSParam() 
+    ! Sets the default parameters
+    implicit none
+    integer i,j
+    IRPoleCheckThreshold=1.0d-5
+    NHelForMCoverHels=5
+    PrecisionVirtualAtRunTime=1d-3
+    Virt_fraction=1d0
+    QED_squared_selected=-1
+    QCD_squared_selected=-1
+    Min_virt_fraction=0.005d0
+    separate_flavour_configs=.false.
+    use_poly_virtual=.true.
+    IncludeBornContributions=.true.
+    SelectedContributionTypes(0)=0
+    VetoedContributionTypes(0)=0
+    do i=1, maxContribsSelected
+       SelectedContributionTypes(I)=-1
+       VetoedContributionTypes(I)=-1
+    enddo
+    do j=1,maxCouplingTypes
+       SelectedCouplingOrders(j,0) = 0
+    enddo
+    do j=1,maxCouplingsSelected
+       do i=1,maxCouplingTypes
+          SelectedCouplingOrders(i,j) = -1
+       enddo
+    enddo
+  end subroutine DefaultFKSParam
+
+end module FKSParams

=== removed file 'Template/NLO/SubProcesses/FKSParams.inc'
=== modified file 'Template/NLO/SubProcesses/MC_integer.f'
--- Template/NLO/SubProcesses/MC_integer.f	2018-06-04 07:36:30 +0000
+++ Template/NLO/SubProcesses/MC_integer.f	2020-05-23 19:29:41 +0000
@@ -58,6 +58,7 @@
       data firsttime/maxdim*.true./
       data realfirsttime/.true./
       character*1 cdum
+      character*3 action
       integer nintervals(maxdim),maxintervals,niint_thisd
       parameter (maxintervals=200)
       integer ncall(0:maxintervals,maxdim)
@@ -95,7 +96,7 @@
             do i=1,this_dim-1 ! skip the lines not needed for 'this_dim'
                read(52,*,end=999,err=999) cdum
             enddo
-            read(52,*,end=999,err=999)
+            read(52,*,end=999,err=999) action,
      &           (grid(i,this_dim),i=0,nintervals(this_dim)) ! here is what we want
             do i=this_dim+1,maxdim ! make sure that there are enough lines in this file
                read(52,*,end=999,err=999) cdum
@@ -266,7 +267,7 @@
       enddo
       open(unit=52,file='grid.MC_integer',status='unknown',err=999)
       do this_dim=1,maxdim
-         write(52,*) (grid(i,this_dim),i=0,nintervals(this_dim))
+         write(52,*) 'AVE',(grid(i,this_dim),i=0,nintervals(this_dim))
       enddo
       close(52)
 c

=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2018-09-21 13:46:51 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2020-05-23 19:29:41 +0000
@@ -2,6 +2,8 @@
 c**************************************************************************
 c     This is the driver for the whole calulation
 c**************************************************************************
+      use mint_module
+      use FKSParams
       implicit none
 C
 C     CONSTANTS
@@ -9,7 +11,7 @@
       double precision pi, zero
       parameter (pi=3.1415926535897932385d0)
       parameter (zero = 0d0)
-      integer npoints, npointsChecked
+      integer npointsChecked
       integer i, j, k
       integer return_code
       double precision tolerance, tolerance_default
@@ -63,8 +65,6 @@
       common /to_polecheck/force_polecheck, polecheck_passed
       integer ret_code_ml
       common /to_ret_code/ret_code_ml
-      include 'FKSParams.inc'
-      include 'mint.inc'
       
 C-----
 C  BEGIN CODE
@@ -462,16 +462,3 @@
       return 
       end
 
-      subroutine outfun(p, a, b, i)
-c     just a dummy subroutine
-      implicit none
-      include 'nexternal.inc'
-      double precision p(0:3, nexternal), a, b
-      integer i
-      write(*,*) 'THIS FUNCTION SHOULD NEVER BE CALLED'
-      return
-      end
-
-      
-      subroutine initplot
-      end

=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2018-04-16 14:08:47 +0000
+++ Template/NLO/SubProcesses/cuts.f	2020-05-23 19:29:41 +0000
@@ -139,7 +139,7 @@
                            return
                         endif
                      endif
-                     if (ipdg(i).eq.-ipdg(j)) then
+                     if (ipdg_reco(i).eq.-ipdg_reco(j)) then
                         if (drll_sf.gt.0d0) then
                            if (R2_04(p_reco(0,i),p_reco(0,j)).lt.drll_sf**2) then
                               passcuts_user=.false.

=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2018-09-21 13:46:51 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2020-05-23 19:29:41 +0000
@@ -3,6 +3,8 @@
 c     This is the driver for the whole calculation
 c**************************************************************************
       use extra_weights
+      use mint_module
+      use FKSParams
       implicit none
 C
 C     CONSTANTS
@@ -11,13 +13,10 @@
       parameter       (ZERO = 0d0)
       include 'nexternal.inc'
       include 'genps.inc'
-      INTEGER    ITMAX,   NCALL
-
-      common/citmax/itmax,ncall
 C
 C     LOCAL
 C
-      integer i,j,k,l,l1,l2
+      integer i,j,k,l,l1,l2,kchan
       character*130 buf
 c
 c     Global
@@ -27,8 +26,8 @@
       include 'coupl.inc'
       
 c     Vegas stuff
-      integer         ndim
-      common/tosigint/ndim
+      integer         nndim
+      common/tosigint/nndim
 
       real*8 sigint
       external sigint
@@ -43,21 +42,11 @@
       double precision xratmax
       common/ccheckcnt/i_momcmp_count,xratmax
 
+      character*4      abrv
+      common /to_abrv/ abrv
       integer n_mp, n_disc
-c For MINT:
-      include "mint.inc"
-      integer nhits_in_grids(maxchannels)
-      real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
-     $     ,ndimmax,maxchannels),ymax_virt(0:maxchannels),ans(nintegrals
-     $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
-     $     ,0:maxchannels),x(ndimmax),itmax_fl
       integer ixi_i,iphi_i,iy_ij,vn
-      integer ifold(ndimmax) 
-      common /cifold/ifold
-      integer ifold_energy,ifold_phi,ifold_yij
-      common /cifoldnumbers/ifold_energy,ifold_phi,ifold_yij
       logical putonshell
-      integer imode,dummy
       logical unwgt
       double precision evtsgn
       common /c_unwgt/evtsgn,unwgt
@@ -72,20 +61,12 @@
 
       double precision virtual_over_born
       common/c_vob/virtual_over_born
-      double precision average_virtual(0:n_ave_virt,maxchannels)
-     $     ,virtual_fraction(maxchannels)
-      common/c_avg_virt/average_virtual,virtual_fraction
       include 'orders.inc'
-      integer              n_ord_virt
-      common /c_n_ord_virt/n_ord_virt
 
 c timing statistics
       include "timing_variables.inc"
       real*4 tOther, tTot
 
-c general MadFKS parameters
-      include "FKSParams.inc"
-
 c applgrid
       integer iappl
       common /for_applgrid/ iappl
@@ -95,13 +76,17 @@
       integer ntot_granny,derntot,ncase(0:6)
       common /c_granny_counters/ ntot_granny,ncase,derntot,deravg,derstd
      &     ,dermax,xi_i_fks_ev_der_max,y_ij_fks_ev_der_max
-      logical              fixed_order,nlo_ps
-      common /c_fnlo_nlops/fixed_order,nlo_ps
 
+      logical useitmax
+      common/cuseitmax/useitmax
 
 C-----
 C  BEGIN CODE
-C-----  
+C-----
+c Write the process PID in the log.txt files (i.e., to the screen)
+      write (*,*) getpid()
+      
+      useitmax=.false. ! to be overwritten in open_output_files.f if need be
 c
 c     Setup the timing variable
 c
@@ -112,6 +97,7 @@
 c     Read general MadFKS parameters
 c
       call FKSParamReader(paramFileName,.TRUE.,.FALSE.)
+      min_virt_fraction_mint=min_virt_fraction
       do kchan=1,maxchannels
          do i=0,n_ave_virt
             average_virtual(i,kchan)=0d0
@@ -153,7 +139,7 @@
 c     Get user input
 c
       write(*,*) "getting user params"
-      call get_user_params(ncall,itmax,imode)
+      call get_user_params(ncalls0,itmax,imode)
       if(imode.eq.0)then
         flat_grid=.true.
       else
@@ -162,6 +148,7 @@
       ndim = 3*(nexternal-nincoming)-4
       if (abs(lpp(1)) .ge. 1) ndim=ndim+1
       if (abs(lpp(2)) .ge. 1) ndim=ndim+1
+      nndim=ndim
 c Don't proceed if muF1#muF2 (we need to work out the relevant formulae
 c at the NLO)
       if( ( fixed_fac_scale .and.
@@ -172,7 +159,7 @@
         write(*,*)'NLO computations require muF1=muF2'
         stop
       endif
-      write(*,*) "about to integrate ", ndim,ncall,itmax
+      write(*,*) "about to integrate ", ndim,ncalls0,itmax
 c APPLgrid
       if (imode.eq.0) iappl=0 ! overwrite when starting completely fresh
       if(iappl.ne.0) then
@@ -184,6 +171,23 @@
          call find_iproc_map
          write(6,*) "   ... done."
       endif
+      if (abrv(1:4).eq.'virt') then
+         only_virt=.true.
+      else
+         only_virt=.false.
+      endif
+c     Prepare the MINT folding
+      do j=1,ndimmax
+         if (j.le.ndim) then
+            ifold(j)=1
+         else
+            ifold(j)=0
+         endif
+      enddo
+      ifold_energy=ndim-2
+      ifold_yij=ndim-1
+      ifold_phi=ndim
+c      
       i_momcmp_count=0
       xratmax=0.d0
       unwgt=.false.
@@ -194,38 +198,8 @@
             doreweight=.false.
             do_rwgt_scale=.false.
             do_rwgt_pdf=.false.
-            do kchan=1,nchans
-               do i=1,ndimmax
-                  do j=0,nintervals
-                     xgrid(j,i,kchan)=0.d0
-                  enddo
-               enddo
-            enddo
          else
-            doreweight=do_rwgt_scale.or.do_rwgt_pdf
-c to restore grids:
-            open (unit=12, file='mint_grids',status='old')
-            ans(1,0)=0d0
-            unc(1,0)=0d0
-            do kchan=1,nchans
-               do j=0,nintervals
-                  read (12,*) (xgrid(j,i,kchan),i=1,ndim)
-               enddo
-               do j=1,nintervals_virt
-                  do k=0,n_ord_virt
-                     read (12,*) (ave_virt(j,i,k,kchan),i=1,ndim)
-                  enddo
-               enddo
-               read(12,*) ans(1,kchan),unc(1,kchan),dummy,dummy
-     $              ,nhits_in_grids(kchan)
-               read(12,*) virtual_fraction(kchan),average_virtual(0
-     $              ,kchan)
-               ans(1,0)=ans(1,0)+ans(1,kchan)
-               unc(1,0)=unc(1,0)+unc(1,kchan)**2
-            enddo
-            unc(1,0)=sqrt(unc(1,0))
-            close (12)
-            write (*,*) "Update iterations and points to",itmax,ncall
+            doreweight=do_rwgt_scale.or.do_rwgt_pdf.or.store_rwgt_info
          endif
 c
          write (*,*) 'imode is ',imode
@@ -236,44 +210,10 @@
                virtual_fraction(kchan)=1d0
             enddo
          endif
-C check for zero cross-section
-C if restoring grids corresponding to sigma=0, just terminate the run
-         if (imode.ne.0.and.ans(1,0).eq.0d0.and.unc(1,0).eq.0d0) then
-            call initplot()
-            call close_run_zero_res(ncall, itmax, ndim)
-            stop
-         endif
-         call mint(sigint,ndim,ncall,itmax,imode,xgrid,ymax
-     $        ,ymax_virt,ans,unc,chi2,nhits_in_grids)
+         call mint(sigint)
          call topout
          call deallocate_weight_lines
-         write(*,*)'Final result [ABS]:',ans(1,0),' +/-',unc(1,0)
-         write(*,*)'Final result:',ans(2,0),' +/-',unc(2,0)
-         write(*,*)'chi**2 per D.o.F.:',chi2(1,0)
-         open(unit=58,file='results.dat',status='unknown')
-         do kchan=0,nchans
-            write(58,*) ans(1,kchan),unc(2,kchan),0d0,0,0,0,0,0d0,0d0
-     $           ,ans(2,kchan)
-         enddo
-         close(58)
 c
-c to save grids:
-         open (unit=12, file='mint_grids',status='unknown')
-         do kchan=1,nchans
-            do j=0,nintervals
-               write (12,*) (xgrid(j,i,kchan),i=1,ndim)
-            enddo
-            do j=1,nintervals_virt
-               do k=0,n_ord_virt
-                  write (12,*) (ave_virt(j,i,k,kchan),i=1,ndim)
-               enddo
-            enddo
-            write (12,*) ans(1,kchan),unc(1,kchan),ncall,itmax
-     $           ,nhits_in_grids(kchan)
-            write (12,*) virtual_fraction(kchan),average_virtual(0
-     $           ,kchan)
-         enddo
-         close (12)
       else
          write (*,*) 'Unknown imode',imode
          stop
@@ -351,7 +291,7 @@
       open (unit=12, file='res.dat',status='unknown')
       do kchan=0,nchans
          write (12,*)ans(1,kchan),unc(1,kchan),ans(2,kchan),unc(2,kchan)
-     $        ,itmax,ncall,tTot
+     $        ,itmax,ncalls0,tTot
       enddo
       close(12)
 
@@ -392,9 +332,9 @@
       double precision function sigint(xx,vegas_wgt,ifl,f)
       use weight_lines
       use extra_weights
+      use mint_module
       implicit none
       include 'nexternal.inc'
-      include 'mint.inc'
       include 'nFKSconfigs.inc'
       include 'run.inc'
       include 'orders.inc'
@@ -410,8 +350,8 @@
       integer             ini_fin_fks(maxchannels)
       common/fks_channels/ini_fin_fks
       data sum /.false./
-      integer         ndim
-      common/tosigint/ndim
+      integer         nndim
+      common/tosigint/nndim
       logical       nbody
       common/cnbody/nbody
       double precision p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
@@ -419,9 +359,6 @@
       common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
       double precision p_born(0:3,nexternal-1)
       common /pborn/   p_born
-      double precision           virt_wgt_mint(0:n_ave_virt),
-     &                           born_wgt_mint(0:n_ave_virt)
-      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
       double precision virtual_over_born
       common/c_vob/virtual_over_born
       logical                calculatedBorn
@@ -434,6 +371,9 @@
       common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
       integer ini_fin_fks_map(0:2,0:fks_configs)
       save ini_fin_fks_map
+      if (new_point .and. ifl.ne.2) then
+         pass_cuts_check=.false.
+      endif
       if (firsttime) then
          firsttime=.false.
          call setup_ini_fin_fks_map(ini_fin_fks_map)
@@ -479,14 +419,16 @@
       else
          jac=0.5d0
       endif
-      call generate_momenta(ndim,iconfig,jac,x,p)
+      call generate_momenta(nndim,iconfig,jac,x,p)
       if (p_born(0,1).lt.0d0) goto 12
       call compute_prefactors_nbody(vegas_wgt)
       call set_cms_stuff(izero)
       passcuts_nbody=passcuts(p1_cnt(0,1,0),rwgt)
       if (passcuts_nbody) then
+         pass_cuts_check=.true.
          if (ickkw.eq.3) call set_FxFx_scale(1,p1_cnt(0,1,0))
          call set_alphaS(p1_cnt(0,1,0))
+         call include_multichannel_enhance(1)
          if (abrv(1:2).ne.'vi') then
             call compute_born
          endif
@@ -520,7 +462,7 @@
          wgt_me_real=0d0
          jac=MC_int_wgt
          call update_fks_dir(iFKS)
-         call generate_momenta(ndim,iconfig,jac,x,p)
+         call generate_momenta(nndim,iconfig,jac,x,p)
          if (p_born(0,1).lt.0d0) cycle
          call compute_prefactors_n1body(vegas_wgt,jac)
          call set_cms_stuff(izero)
@@ -528,9 +470,11 @@
          call set_cms_stuff(mohdr)
          passcuts_n1body=passcuts(p,rwgt)
          if (passcuts_nbody .and. abrv.ne.'real') then
+            pass_cuts_check=.true.
             call set_cms_stuff(izero)
             if (ickkw.eq.3) call set_FxFx_scale(2,p1_cnt(0,1,0))
             call set_alphaS(p1_cnt(0,1,0))
+            call include_multichannel_enhance(3)
             call compute_soft_counter_term(0d0)
             call set_cms_stuff(ione)
             call compute_collinear_counter_term(0d0)
@@ -538,9 +482,11 @@
             call compute_soft_collinear_counter_term(0d0)
          endif
          if (passcuts_n1body) then
+            pass_cuts_check=.true.
             call set_cms_stuff(mohdr)
             if (ickkw.eq.3) call set_FxFx_scale(3,p)
             call set_alphaS(p)
+            call include_multichannel_enhance(2)
             call compute_real_emission(p,1d0)
          endif
       enddo
@@ -703,26 +649,26 @@
       end
 
       subroutine update_vegas_x(xx,x)
+      use mint_module
       implicit none
-      include 'mint.inc'
       integer i
       double precision xx(ndimmax),x(99),ran2
       external ran2
-      integer ndim
-      common/tosigint/ndim
-      character*4 abrv
+      integer         nndim
+      common/tosigint/nndim
+      character*4      abrv
       common /to_abrv/ abrv
       do i=1,99
          if (abrv.eq.'born'.or.abrv(1:2).eq.'vi') then
-            if(i.le.ndim-3)then
+            if(i.le.nndim-3)then
                x(i)=xx(i)
-            elseif(i.le.ndim) then
+            elseif(i.le.nndim) then
                x(i)=ran2()      ! Choose them flat when not including real-emision
             else
                x(i)=0.d0
             endif
          else
-            if(i.le.ndim)then
+            if(i.le.nndim)then
                x(i)=xx(i)
             else
                x(i)=0.d0
@@ -733,10 +679,11 @@
       end
 
 c
-      subroutine get_user_params(ncall,itmax,irestart)
+      subroutine get_user_params(ncall,nitmax,irestart)
 c**********************************************************************
 c     Routine to get user specified parameters for run
 c**********************************************************************
+      use mint_module
       implicit none
 c
 c     Constants
@@ -746,16 +693,15 @@
       include 'nFKSconfigs.inc'
       include 'fks_info.inc'
       include 'run.inc'
-      include 'mint.inc'
       include 'orders.inc'
 c
 c     Arguments
 c
-      integer ncall,itmax
+      integer ncall,nitmax
 c
 c     Local
 c
-      integer i, j
+      integer i, j, kchan
       double precision dconfig(maxchannels)
 c
 c     Global
@@ -789,8 +735,6 @@
       character * 70 idstring
       logical savegrid
 
-      logical usexinteg,mint
-      common/cusexinteg/usexinteg,mint
       logical unwgt
       double precision evtsgn
       common /c_unwgt/evtsgn,unwgt
@@ -805,7 +749,6 @@
 c-----
 c  Begin Code
 c-----
-      mint=.true.
       unwgt=.false.
       open (unit=83,file='input_app.txt',status='old')
       done=.false.
@@ -817,8 +760,8 @@
             read(buffer,*) ncall
             write (*,*) 'Number of phase-space points per iteration:',ncall
          elseif(buffer(1:11).eq.'NITERATIONS') then
-            read(buffer(14:),*) itmax
-            write (*,*) 'Maximum number of iterations is:',itmax
+            read(buffer(14:),*) nitmax
+            write (*,*) 'Maximum number of iterations is:',nitmax
          elseif(buffer(1:8).eq.'ACCURACY') then
             read(buffer(11:),*) accuracy
             write (*,*) 'Desired accuracy is:',accuracy

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2018-09-21 13:46:51 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2020-05-23 19:29:41 +0000
@@ -3,6 +3,8 @@
 c     This is the driver for the whole calculation
 c**************************************************************************
       use extra_weights
+      use mint_module
+      use FKSParams
       implicit none
 C
 C     CONSTANTS
@@ -11,16 +13,13 @@
       parameter       (ZERO = 0d0)
       include 'nexternal.inc'
       include 'genps.inc'
-      INTEGER    ITMAX,   NCALL
-
-      common/citmax/itmax,ncall
       integer ncall_virt,ncall_novi
       character*4 abrv
       common /to_abrv/ abrv
 C
 C     LOCAL
 C
-      integer i,j,k,l,l1,l2,ndim,nevts
+      integer i,j,k,l,l1,l2,nndim,nevts
 
       integer lunlhe
       parameter (lunlhe=98)
@@ -30,10 +29,9 @@
 cc
       include 'run.inc'
       include 'coupl.inc'
-      include "mint.inc"
       
 c Vegas stuff
-      common/tosigint/ndim
+      common/tosigint/nndim
 
       real*8 sigintF
       external sigintF
@@ -47,34 +45,17 @@
 
       double precision virtual_over_born
       common/c_vob/virtual_over_born
-      double precision average_virtual(0:n_ave_virt,maxchannels)
-     $     ,virtual_fraction(maxchannels)
-      common/c_avg_virt/average_virtual,virtual_fraction
       include 'orders.inc'
-      integer              n_ord_virt
-      common /c_n_ord_virt/n_ord_virt
 
       double precision weight,event_weight,inv_bias
       character*7 event_norm
       common /event_normalisation/event_norm
-c For MINT:
-      real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
-     $     ,ndimmax,maxchannels),ymax_virt(0:maxchannels),ans(nintegrals
-     $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
-     $     ,0:maxchannels),x(ndimmax)
-      integer ixi_i,iphi_i,iy_ij,vn,nhits_in_grids(maxchannels)
-      integer ifold(ndimmax) 
-      common /cifold/ifold
-      integer ifold_energy,ifold_phi,ifold_yij
-      common /cifoldnumbers/ifold_energy,ifold_phi,ifold_yij
+      integer ixi_i,iphi_i,iy_ij,vn
       logical putonshell
-      logical only_virt
-      integer imode
-      common /c_imode/imode,only_virt
       logical unwgt
       double precision evtsgn
       common /c_unwgt/evtsgn,unwgt
-      double precision ran2
+      double precision ran2,x(ndimmax)
       external ran2
       
       integer ifile,ievents
@@ -94,11 +75,6 @@
       include "timing_variables.inc"
       real*4 tOther, tTot
 
-c general MadFKS parameters
-      include "FKSParams.inc"
-      logical              fixed_order,nlo_ps
-      common /c_fnlo_nlops/fixed_order,nlo_ps
-
       double precision deravg,derstd,dermax,xi_i_fks_ev_der_max
      &     ,y_ij_fks_ev_der_max
       integer ntot_granny,derntot,ncase(0:6)
@@ -110,6 +86,9 @@
 C-----
 C  BEGIN CODE
 C-----  
+c Write the process PID in the log.txt files (i.e., to the screen)
+      write (*,*) getpid()
+
       call cpu_time(tBefore)
       fixed_order=.false.
       nlo_ps=.true.
@@ -122,6 +101,7 @@
 c     Read general MadFKS parameters
 c
       call FKSParamReader(paramFileName,.TRUE.,.FALSE.)
+      min_virt_fraction_mint=min_virt_fraction
       do i=0,n_ave_virt
          average_virtual(i,1)=0d0
       enddo
@@ -158,11 +138,11 @@
 c     Get user input
 c
       write(*,*) "getting user params"
-      call get_user_params(ncall,itmax,imode,
+      call get_user_params(ncalls0,itmax,
      &     ixi_i,iphi_i,iy_ij,SHsep)
 c Only do the reweighting when actually generating the events
       if (imode.eq.2) then
-         doreweight=do_rwgt_scale.or.do_rwgt_pdf
+         doreweight=do_rwgt_scale.or.do_rwgt_pdf.or.store_rwgt_info
       else
          doreweight=.false.
          do_rwgt_scale=.false.
@@ -182,6 +162,7 @@
       ndim = 3*(nexternal-nincoming)-4
       if (abs(lpp(1)) .ge. 1) ndim=ndim+1
       if (abs(lpp(2)) .ge. 1) ndim=ndim+1
+      nndim=ndim
 c Don''t proceed if muF1#muF2 (we need to work out the relevant formulae
 c at the NLO)
       if( ( fixed_fac_scale .and.
@@ -192,146 +173,46 @@
         write(*,*)'NLO computations require muF1=muF2'
         stop
       endif
-      write(*,*) "about to integrate ", ndim,ncall,itmax,iconfig
+      write(*,*) "about to integrate ", ndim,ncalls0,itmax,iconfig
       i_momcmp_count=0
       xratmax=0.d0
       unwgt=.false.
+c     Prepare the MINT folding
+      do j=1,ndimmax
+         if (j.le.ndim) then
+            ifold(j)=1
+         else
+            ifold(j)=0
+         endif
+      enddo
+      ifold_energy=ndim-2
+      ifold_yij=ndim-1
+      ifold_phi=ndim
+      ifold(ifold_energy)=ixi_i
+      ifold(ifold_phi)=iphi_i
+      ifold(ifold_yij)=iy_ij
 
 c*************************************************************
 c     setting of the grids
 c*************************************************************
       if (imode.eq.-1.or.imode.eq.0) then
-         if(imode.eq.0)then
-c initialize grids
-            do j=0,nintervals
-               do i=1,ndimmax
-                  xgrid(j,i,1)=0.d0
-               enddo
-            enddo
-         else
-c to restore grids:
-            open (unit=12, file='mint_grids',status='old')
-            do j=0,nintervals
-               read (12,*) (xgrid(j,i,1),i=1,ndim)
-            enddo
-            do j=1,nintervals_virt
-               do k=0,n_ord_virt
-                  read (12,*) (ave_virt(j,i,k,1),i=1,ndim)
-               enddo
-            enddo
-            read (12,*) (ans(i,1),i=1,nintegrals)
-            read (12,*) ifold_energy,ifold_phi,ifold_yij
-            read (12,*) virtual_fraction(1),average_virtual(0,1)
-            close (12)
-         endif
-c
          write (*,*) 'imode is ',imode
-         call mint(sigintF,ndim,ncall,itmax,imode,xgrid,ymax,ymax_virt
-     $        ,ans,unc,chi2,nhits_in_grids)
+         call mint(sigintF)
          call deallocate_weight_lines
-         open(unit=58,file='res_0',status='unknown')
-         write(58,*)'Final result [ABS]:',ans(1,1),' +/-',unc(1,1)
-         write(58,*)'Final result:',ans(2,1),' +/-',unc(2,1)
-         close(58)
-         write(*,*)'Final result [ABS]:',ans(1,1),' +/-',unc(1,1)
-         write(*,*)'Final result:',ans(2,1),' +/-',unc(2,1)
-         write(*,*)'chi**2 per D.o.F.:',chi2(1,1)
          open(unit=58,file='results.dat',status='unknown')
          write(58,*) ans(1,1),unc(2,1),0d0,0,0,0,0,0d0,0d0,ans(2,1)
          close(58)
-c
-c to save grids:
-         open (unit=12, file='mint_grids',status='unknown')
-         do j=0,nintervals
-            write (12,*) (xgrid(j,i,1),i=1,ndim)
-         enddo
-         do j=1,nintervals_virt
-            do k=0,n_ord_virt
-               write (12,*) (ave_virt(j,i,k,1),i=1,ndim)
-            enddo
-         enddo
-         write (12,*) (ans(i,1),i=1,nintegrals)
-         write (12,*) ifold_energy,ifold_phi,ifold_yij
-         write (12,*) virtual_fraction(1),average_virtual(0,1)
-         close (12)
-
 c*************************************************************
 c     computation of upper bounding envelope
 c*************************************************************
       elseif(imode.eq.1) then
-c to restore grids:
-         open (unit=12, file='mint_grids',status='old')
-         do j=0,nintervals
-            read (12,*) (xgrid(j,i,1),i=1,ndim)
-         enddo
-         do j=1,nintervals_virt
-            do k=0,n_ord_virt
-               read (12,*) (ave_virt(j,i,k,1),i=1,ndim)
-            enddo
-         enddo
-         read (12,*) (ans(i,1),i=1,nintegrals)
-         read (12,*) ifold_energy,ifold_phi,ifold_yij
-         read (12,*) virtual_fraction(1),average_virtual(0,1)
-         close (12)
-
-c Prepare the MINT folding
-         do j=1,ndimmax
-            if (j.le.ndim) then
-               ifold(j)=1
-            else
-               ifold(j)=0
-            endif
-         enddo
-         ifold(ifold_energy)=ixi_i
-         ifold(ifold_phi)=iphi_i
-         ifold(ifold_yij)=iy_ij
-         
          write (*,*) 'imode is ',imode
-         call mint(sigintF,ndim,ncall,itmax,imode,xgrid,ymax,ymax_virt
-     $        ,ans,unc,chi2,nhits_in_grids)
+         call mint(sigintF)
          call deallocate_weight_lines
-         
-c If integrating the virtuals alone, we include the virtuals in
-c ans(1). Therefore, no need to have them in ans(5) and we have to set
-c them to zero.
-         if (only_virt) then
-            ans(3,1)=0d0 ! virtual Xsec
-            ans(5,1)=0d0 ! ABS virtual Xsec
-         endif
-
-         open(unit=58,file='res_1',status='unknown')
-         write(58,*)'Final result [ABS]:',ans(1,1)+ans(5,1),' +/-'
-     $        ,sqrt(unc(1,1)**2+unc(5,1)**2)
-         write(58,*)'Final result:',ans(2,1),' +/-',unc(2,1)
-         close(58)
-         write(*,*)'Final result [ABS]:',ans(1,1)+ans(5,1),' +/-'
-     $        ,sqrt(unc(1,1)**2+unc(5,1)**2)
-         write(*,*)'Final result:',ans(2,1),' +/-',unc(2,1)
-         write(*,*)'chi**2 per D.o.F.:',chi2(1,1)
-c write the results.dat file 
          open(unit=58,file='results.dat',status='unknown')
-         write(58,*)ans(1,1)+ans(5,1), unc(2,1), 0d0, 0, 0, 0, 0, 0d0 ,0d0, ans(2,1) 
+         write(58,*) ans(1,1)+ans(5,1),unc(2,1),0d0,0,0,0,0,0d0,0d0
+     $        ,ans(2,1) 
          close(58)
-
-c to save grids:
-         open (unit=12, file='mint_grids',status='unknown')
-         write (12,*) (xgrid(0,i,1),i=1,ndim)
-         do j=1,nintervals
-            write (12,*) (xgrid(j,i,1),i=1,ndim)
-            write (12,*) (ymax(j,i,1),i=1,ndim)
-         enddo
-         do j=1,nintervals_virt
-            do k=0,n_ord_virt
-               write (12,*) (ave_virt(j,i,k,1),i=1,ndim)
-            enddo
-         enddo
-         write (12,*) ymax_virt(1)
-         write (12,*) (ifold(i),i=1,ndim)
-         write (12,*) (ans(i,1),i=1,nintegrals)
-         write (12,*) (unc(i,1),i=1,nintegrals)
-         write (12,*) virtual_fraction(1),average_virtual(0,1)
-         close (12)
-
 c*************************************************************
 c     event generation
 c*************************************************************
@@ -354,30 +235,15 @@
      &           'No events needed for this channel...skipping it'
             stop
          endif
-         ncall=nevts ! Update ncall with the number found in 'nevts'
-
-c to restore grids:
-         open (unit=12, file='mint_grids',status='unknown')
-         read (12,*) (xgrid(0,i,1),i=1,ndim)
-         do j=1,nintervals
-            read (12,*) (xgrid(j,i,1),i=1,ndim)
-            read (12,*) (ymax(j,i,1),i=1,ndim)
-         enddo
-         do j=1,nintervals_virt
-            do k=0,n_ord_virt
-               read (12,*) (ave_virt(j,i,k,1),i=1,ndim)
-            enddo
-         enddo
-         read (12,*) ymax_virt(1)
-         read (12,*) (ifold(i),i=1,ndim)
-         read (12,*) (ans(i,1),i=1,nintegrals)
-         read (12,*) (unc(i,1),i=1,nintegrals)
-         read (12,*) virtual_fraction(1),average_virtual(0,1)
-         close (12)
+         ncalls0=nevts ! Update ncall with the number found in 'nevts'
+
+c     to restore grids:
+
+         call read_grids_from_file
 
 c determine how many events for the virtual and how many for the no-virt
-         ncall_virt=int(ans(5,1)/(ans(1,1)+ans(5,1)) * ncall)
-         ncall_novi=ncall-ncall_virt
+         ncall_virt=int(ans(5,1)/(ans(1,1)+ans(5,1)) * ncalls0)
+         ncall_novi=ncalls0-ncall_virt
 
          write (*,*) "Generating virt :: novi approx.",ncall_virt
      $        ,ncall_novi
@@ -386,13 +252,13 @@
 
 c fill the information for the write_header_init common block
          ifile=lunlhe
-         ievents=ncall
+         ievents=ncalls0
          inter=ans(2,1)
          absint=ans(1,1)+ans(5,1)
          uncer=unc(2,1)
 
          if (event_norm(1:4).ne.'bias') then
-            weight=(ans(1,1)+ans(5,1))/ncall
+            weight=(ans(1,1)+ans(5,1))/ncalls0
          else
             weight=event_weight
          endif
@@ -405,25 +271,25 @@
 
          write (*,*) 'imode is ',imode
          vn=-1
-         call gen(sigintF,ndim,xgrid,ymax,ymax_virt,0,x,vn)
-         do j=1,ncall
+         call gen(sigintF,0,vn,x)
+         do j=1,ncalls0
             if (abrv(1:4).eq.'born') then
                vn=3
-               call gen(sigintF,ndim,xgrid,ymax,ymax_virt,1,x,vn)
+               call gen(sigintF,1,vn,x)
             else
                if (ran2().lt.ans(5,1)/(ans(1,1)+ans(5,1)) .or. only_virt) then
                   abrv='virt'
                   if (only_virt) then
                      vn=2
-                     call gen(sigintF,ndim,xgrid,ymax,ymax_virt,1,x,vn)
+                     call gen(sigintF,1,vn,x)
                   else
                      vn=1
-                     call gen(sigintF,ndim,xgrid,ymax,ymax_virt,1,x,vn)
+                     call gen(sigintF,1,vn,x)
                   endif
                else
                   abrv='novi'
                   vn=2
-                  call gen(sigintF,ndim,xgrid,ymax,ymax_virt,1,x,vn)
+                  call gen(sigintF,1,vn,x)
                endif
             endif
 c Randomly pick the contribution that will be written in the event file
@@ -438,14 +304,7 @@
          enddo
          call deallocate_weight_lines
          vn=-1
-         call gen(sigintF,ndim,xgrid,ymax,ymax_virt,3,x,vn)
-         write (*,*) 'Generation efficiencies:',x(1),x(4)
-c Uncomment the next to lines to print the integral from the PS points
-c trown during event generation. This corresponds only to the cross
-c section if these points are thrown flat, so not using the xmmm() stuff
-c in mint.
-c         write (*,*) 'Integral from novi points computed',x(2),x(3)
-c         write (*,*) 'Integral from virt points computed',x(5),x(6)
+         call gen(sigintF,3,vn,x) ! print counters generation efficiencies
          write (lunlhe,'(a)') "</LesHouchesEvents>"
          close(lunlhe)
       endif
@@ -528,10 +387,10 @@
 
       open (unit=12, file='res.dat',status='unknown')
       if (imode.eq.0) then
-         write (12,*)ans(1,1),unc(1,1),ans(2,1),unc(2,1),itmax,ncall,tTot
+         write (12,*)ans(1,1),unc(1,1),ans(2,1),unc(2,1),itmax,ncalls0,tTot
       else
          write (12,*)ans(1,1)+ans(5,1),sqrt(unc(1,1)**2+unc(5,1)**2),ans(2,1)
-     $        ,unc(2,1),itmax,ncall,tTot
+     $        ,unc(2,1),itmax,ncalls0,tTot
       endif
       close(12)
 
@@ -565,25 +424,25 @@
       end
 
 
-      subroutine get_user_params(ncall,itmax,
-     &     imode,ixi_i,iphi_i,iy_ij,SHsep)
+      subroutine get_user_params(ncall,nitmax,
+     &     ixi_i,iphi_i,iy_ij,SHsep)
 c**********************************************************************
 c     Routine to get user specified parameters for run
 c**********************************************************************
+      use mint_module
       implicit none
 c
 c     Constants
 c
       include 'nexternal.inc'
       include 'genps.inc'
-      include 'mint.inc'
       include 'nFKSconfigs.inc'
       include 'fks_info.inc'
       include 'run.inc'
 c
 c     Arguments
 c
-      integer ncall,itmax,jconfig
+      integer ncall,nitmax
 c
 c     Local
 c
@@ -633,19 +492,14 @@
 c
 c MINT stuff
 c
-      integer imode,ixi_i,iphi_i,iy_ij
-
-      logical usexinteg,mint
-      common/cusexinteg/usexinteg,mint
+      integer ixi_i,iphi_i,iy_ij
 
 c-----
 c  Begin Code
 c-----
-      mint=.true.
-      usexinteg=.false.
       write(*,'(a)') 'Enter number of events and iterations: '
-      read(*,*) ncall,itmax
-      write(*,*) 'Number of events and iterations ',ncall,itmax
+      read(*,*) ncall,nitmax
+      write(*,*) 'Number of events and iterations ',ncall,nitmax
 
       write(*,'(a)') 'Enter desired fractional accuracy: '
       read(*,*) accuracy
@@ -736,7 +590,7 @@
       write (*,*) 'MINT running mode:',imode
       if (imode.eq.2)then
          write (*,*) 'Generating events, doing only one iteration'
-         itmax=1
+         nitmax=1
       endif
 
       write (*,'(a)') 'Set the three folding parameters for MINT'
@@ -790,8 +644,8 @@
 
       function sigintF(xx,vegas_wgt,ifl,f)
       use weight_lines
+      use mint_module
       implicit none
-      include 'mint.inc'
       include 'nexternal.inc'
       include 'nFKSconfigs.inc'
       include 'run.inc'
@@ -817,25 +671,22 @@
       common/ccalculatedBorn/calculatedBorn
       logical              MCcntcalled
       common/c_MCcntcalled/MCcntcalled
-      double precision           virt_wgt_mint(0:n_ave_virt),
-     &                           born_wgt_mint(0:n_ave_virt)
-      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
       double precision virtual_over_born
       common /c_vob/   virtual_over_born
       logical       nbody
       common/cnbody/nbody
-      integer         ndim
-      common/tosigint/ndim
+      integer         nndim
+      common/tosigint/nndim
       character*4      abrv
       common /to_abrv/ abrv
       double precision p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
      $     ,pswgt_cnt(-2:2),jac_cnt(-2:2)
       common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
-      logical               only_virt
-      integer         imode
-      common /c_imode/imode,only_virt
       double precision       wgt_ME_born,wgt_ME_real
       common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
+      if (new_point .and. ifl.ne.2) then
+         pass_cuts_check=.false.
+      endif
       sigintF=0d0
 c Find the nFKSprocess for which we compute the Born-like contributions
       if (firsttime) then
@@ -887,15 +738,17 @@
          else
             jac=0.5d0
          endif
-         call generate_momenta(ndim,iconfig,jac,x,p)
+         call generate_momenta(nndim,iconfig,jac,x,p)
          if (p_born(0,1).lt.0d0) goto 12
          call compute_prefactors_nbody(vegas_wgt)
          call set_cms_stuff(izero)
          call set_shower_scale_noshape(p,nFKS_picked_nbody*2-1)
          passcuts_nbody=passcuts(p1_cnt(0,1,0),rwgt)
          if (passcuts_nbody) then
+            pass_cuts_check=.true.
             if (ickkw.eq.3) call set_FxFx_scale(1,p1_cnt(0,1,0))
             call set_alphaS(p1_cnt(0,1,0))
+            call include_multichannel_enhance(1)
             if (abrv(1:2).ne.'vi') then
                call compute_born
             endif
@@ -926,7 +779,7 @@
             gfactsf=1.d0
             gfactcl=1.d0
             MCcntcalled=.false.
-            call generate_momenta(ndim,iconfig,jac,x,p)
+            call generate_momenta(nndim,iconfig,jac,x,p)
 c Every contribution has to have a viable set of Born momenta (even if
 c counter-event momenta do not exist).
             if (p_born(0,1).lt.0d0) cycle
@@ -970,11 +823,13 @@
                endif
             endif               
             if (passcuts_nbody .and. abrv.ne.'real') then
+               pass_cuts_check=.true.
 c Include the MonteCarlo subtraction terms
                if (ickkw.ne.4) then
                   call set_cms_stuff(mohdr)
                   if (ickkw.eq.3) call set_FxFx_scale(-3,p)
                   call set_alphaS(p)
+                  call include_multichannel_enhance(4)
                   call compute_MC_subt_term(p,gfactsf,gfactcl,probne)
                else
 c For UNLOPS all real-emission contributions need to be added to the
@@ -989,6 +844,7 @@
                call set_cms_stuff(izero)
                if (ickkw.eq.3) call set_FxFx_scale(-2,p1_cnt(0,1,0))
                call set_alphaS(p1_cnt(0,1,0))
+               call include_multichannel_enhance(3)
                replace_MC_subt=(1d0-gfactsf)*probne
                call compute_soft_counter_term(replace_MC_subt)
                call set_cms_stuff(ione)
@@ -1000,9 +856,11 @@
             endif
 c Include the real-emission contribution.
             if (passcuts_n1body) then
+               pass_cuts_check=.true.
                call set_cms_stuff(mohdr)
                if (ickkw.eq.3) call set_FxFx_scale(-3,p)
                call set_alphaS(p)
+               call include_multichannel_enhance(2)
                sudakov_damp=probne
                call compute_real_emission(p,sudakov_damp)
             endif
@@ -1280,26 +1138,26 @@
       end
 
       subroutine update_vegas_x(xx,x)
+      use mint_module
       implicit none
-      include 'mint.inc'
       integer i
       double precision xx(ndimmax),x(99),ran2
       external ran2
-      integer         ndim
-      common/tosigint/ndim
-      character*4 abrv
+      integer         nndim
+      common/tosigint/nndim
+      character*4      abrv
       common /to_abrv/ abrv
       do i=1,99
          if (abrv.eq.'born') then
-            if(i.le.ndim-3)then
+            if(i.le.nndim-3)then
                x(i)=xx(i)
-            elseif(i.le.ndim) then
+            elseif(i.le.nndim) then
                x(i)=ran2()      ! Choose them flat when not including real-emision
             else
                x(i)=0.d0
             endif
          else
-            if(i.le.ndim)then
+            if(i.le.nndim)then
                x(i)=xx(i)
             else
                x(i)=0.d0

=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2018-11-21 15:25:36 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2020-05-23 19:29:41 +0000
@@ -176,13 +176,13 @@
 c This subroutine computes the soft-virtual matrix elements and adds its
 c value to the list of weights using the add_wgt subroutine
       use extra_weights
+      use mint_module
       implicit none
       include 'nexternal.inc'
       include 'coupl.inc'
       include 'run.inc'
       include 'timing_variables.inc'
       include 'orders.inc'
-      include 'mint.inc'
       integer orders(nsplitorders)
       integer iamp, i
       double precision amp_split_virt(amp_split_size),
@@ -203,9 +203,6 @@
       double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
      $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
       common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
-      double precision           virt_wgt_mint(0:n_ave_virt),
-     &                           born_wgt_mint(0:n_ave_virt)
-      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
       double precision   xiimax_cnt(-2:2)
       common /cxiimaxcnt/xiimax_cnt
       double precision  xi_i_hat_ev,xi_i_hat_cnt(-2:2)
@@ -976,23 +973,17 @@
 c i.e. all the nbody contributions. Also initialises the plots and
 c bpower.
       use extra_weights
+      use mint_module
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'genps.inc'
       include 'timing_variables.inc'
-      double precision pi,vegas_wgt,enhance,xnoborn_cnt,xtot
-     $     ,bpower,cpower,tiny
-      data xnoborn_cnt /0d0/
-      integer inoborn_cnt,i
-      data inoborn_cnt /0/
-      double precision wgt_c
+      double precision pi,vegas_wgt
+      integer i
       logical firsttime
       data firsttime /.true./
       parameter (pi=3.1415926535897932385d0)
-      parameter (tiny=1d-6)
-      double precision p_born(0:3,nexternal-1)
-      common/pborn/    p_born
       double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
      $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
       common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
@@ -1009,12 +1000,6 @@
       integer                  ngluons,nquarks(-6:6),nphotons
       common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
      &                  fkssymmetryfactorDeg,ngluons,nquarks,nphotons
-      integer            mapconfig(0:lmaxconfigs), iconfig
-      common/to_mconfigs/mapconfig,                iconfig
-      Double Precision amp2(ngraphs), jamp2(0:ncolor)
-      common/to_amps/  amp2,          jamp2
-      double precision   diagramsymmetryfactor
-      common /dsymfactor/diagramsymmetryfactor
       double precision      f_b,f_nb
       common /factor_nbody/ f_b,f_nb
       integer iappl
@@ -1025,8 +1010,6 @@
       external ran2
       real*8 rndec(10)
       common/crndec/rndec
-      logical              fixed_order,nlo_ps
-      common /c_fnlo_nlops/fixed_order,nlo_ps
       include "appl_common.inc" 
       !!!include "orders.inc"
       call cpu_time(tBefore)
@@ -1043,138 +1026,70 @@
          if (fixed_order) call initplot
          firsttime=.false.
       endif
-c Compute the multi-channel enhancement factor 'enhance'.
-      enhance=1.d0
-      if (p_born(0,1).gt.0d0) then
-         call sborn(p_born,wgt_c)
-      elseif(p_born(0,1).lt.0d0)then
-         enhance=0d0
-      endif
-      if (enhance.eq.0d0)then
-         xnoborn_cnt=xnoborn_cnt+1.d0
-         if(log10(xnoborn_cnt).gt.inoborn_cnt)then
-            write (*,*) 'WARNING: no Born momenta more than 10**',
-     $           inoborn_cnt,'times'
-            inoborn_cnt=inoborn_cnt+1
-         endif
-      else
-         xtot=0d0
-         if (mapconfig(0).eq.0) then
-            write (*,*) 'Fatal error in compute_prefactor_nbody:'/
-     &           /' no Born diagrams ',mapconfig,
-     &           '. Check bornfromreal.inc'
-            write (*,*) 'Is fks_singular compiled correctly?'
-            stop 1
-         endif
-         do i=1, mapconfig(0)
-            xtot=xtot+amp2(mapconfig(i))
-         enddo
-         if (xtot.ne.0d0) then
-            enhance=amp2(mapconfig(iconfig))/xtot
-            enhance=enhance*diagramsymmetryfactor
-         else
-            enhance=0d0
-         endif
-      endif
       call set_cms_stuff(0)
 c f_* multiplication factors for Born and nbody
       f_b=jac_cnt(0)*xinorm_ev/(min(xiimax_ev,xiBSVcut_used)*shat/(16
-     $     *pi**2))*enhance*fkssymmetryfactorBorn*vegas_wgt
+     $     *pi**2))*fkssymmetryfactorBorn*vegas_wgt
       f_nb=f_b
       call cpu_time(tAfter)
       tf_nb=tf_nb+(tAfter-tBefore)
       return
       end
 
-      subroutine compute_prefactors_n1body(vegas_wgt,jac_ev)
-c Compute all relevant prefactors for the real emission and counter
-c terms.
+
+      subroutine include_multichannel_enhance(imode)
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'genps.inc'
-      include 'fks_powers.inc'
-      include 'coupl.inc'
       include 'timing_variables.inc'
-      double precision vegas_wgt,enhance,xnoborn_cnt,xtot
-     &     ,prefact,prefact_cnt_ssc,prefact_deg,prefact_c,prefact_coll
-     &     ,jac_ev,pi,prefact_cnt_ssc_c,prefact_coll_c,prefact_deg_slxi
-     &     ,prefact_deg_sxi,zero,enhance_real
-      parameter (pi=3.1415926535897932385d0, zero=0d0)
+      double precision xnoborn_cnt,xtot,wgt_c,enhance,enhance_real
+     $     ,pas(0:3,nexternal)
       data xnoborn_cnt /0d0/
-      integer inoborn_cnt,i
+      integer inoborn_cnt,i,imode
       data inoborn_cnt /0/
-      double precision wgt_c
       double precision p_born(0:3,nexternal-1)
       common/pborn/    p_born
       double precision p_born_ev(0:3,nexternal-1)
-      common/pborn_ev/p_born_ev
+      common/pborn_ev/ p_born_ev
+      double precision p_ev(0:3,nexternal)
+      common/pev/      p_ev
       double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
      $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
       common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
-      double precision    xi_i_fks_ev,y_ij_fks_ev
-      double precision    p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
-      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
-      integer            i_fks,j_fks
-      common/fks_indices/i_fks,j_fks
-      double precision   xi_i_fks_cnt(-2:2)
-      common /cxiifkscnt/xi_i_fks_cnt
-      double precision  xinorm_ev
-      common /cxinormev/xinorm_ev
-      double precision  xiimax_ev
-      common /cxiimaxev/xiimax_ev
-      double precision   xiimax_cnt(-2:2)
-      common /cxiimaxcnt/xiimax_cnt
-      double precision   xinorm_cnt(-2:2)
-      common /cxinormcnt/xinorm_cnt
-      double precision    delta_used
-      common /cdelta_used/delta_used
-      double precision    xicut_used
-      common /cxicut_used/xicut_used
-      double precision     xiScut_used,xiBSVcut_used
-      common /cxiScut_used/xiScut_used,xiBSVcut_used
-      double precision        ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
-      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
-      double precision         fkssymmetryfactor,fkssymmetryfactorBorn,
-     &                         fkssymmetryfactorDeg
-      integer ngluons,nquarks(-6:6),nphotons
-      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
-     &                  fkssymmetryfactorDeg,ngluons,nquarks,nphotons
       integer            mapconfig(0:lmaxconfigs), iconfig
       common/to_mconfigs/mapconfig,                iconfig
       Double Precision amp2(ngraphs), jamp2(0:ncolor)
       common/to_amps/  amp2,          jamp2
       double precision   diagramsymmetryfactor
       common /dsymfactor/diagramsymmetryfactor
-      logical nocntevents
-      common/cnocntevents/nocntevents
+      double precision      f_b,f_nb
+      common /factor_nbody/ f_b,f_nb
+      double precision     f_r,f_s,f_c,f_dc,f_sc,f_dsc(4)
+      common/factor_n1body/f_r,f_s,f_c,f_dc,f_sc,f_dsc
+      double precision           f_s_MC_S,f_s_MC_H,f_c_MC_S,f_c_MC_H
+     $     ,f_sc_MC_S,f_sc_MC_H,f_MC_S,f_MC_H
+      common/factor_n1body_NLOPS/f_s_MC_S,f_s_MC_H,f_c_MC_S,f_c_MC_H
+     $     ,f_sc_MC_S,f_sc_MC_H,f_MC_S,f_MC_H
+      double precision f_dis_d,f_dis_p,f_dis_l
+      common/factor_dis/f_dis_d,f_dis_p,f_dis_l
       integer igranny,iaunt
       logical granny_chain(-nexternal:nexternal),granny_is_res
      &     ,granny_chain_real_final(-nexternal:nexternal)
       common /c_granny_res/igranny,iaunt,granny_is_res,granny_chain
      &     ,granny_chain_real_final
-      double precision     f_r,f_s,f_c,f_dc,f_sc,f_dsc(4)
-      common/factor_n1body/f_r,f_s,f_c,f_dc,f_sc,f_dsc
-      double precision           f_s_MC_S,f_s_MC_H,f_c_MC_S,f_c_MC_H
-     $     ,f_sc_MC_S,f_sc_MC_H,f_MC_S,f_MC_H
-      common/factor_n1body_NLOPS/f_s_MC_S,f_s_MC_H,f_c_MC_S,f_c_MC_H
-     $     ,f_sc_MC_S,f_sc_MC_H,f_MC_S,f_MC_H
-      ! prefactors for the DIS scheme
-      double precision prefact_dis_d,prefact_dis_p,prefact_dis_l
-      double precision f_dis_d,f_dis_p,f_dis_l
-      common/factor_dis/f_dis_d,f_dis_p,f_dis_l
       logical calculatedBorn
       common/ccalculatedBorn/calculatedBorn
-      double precision pmass(nexternal)
-      include 'pmass.inc'
+
       call cpu_time(tBefore)
+
+c Compute the multi-channel enhancement factor 'enhance'.
       enhance=1.d0
       if (p_born(0,1).gt.0d0) then
          call sborn(p_born,wgt_c)
       elseif(p_born(0,1).lt.0d0)then
          enhance=0d0
       endif
-c Compute the multi-channel enhancement factor 'enhance'.
       if (enhance.eq.0d0)then
          xnoborn_cnt=xnoborn_cnt+1.d0
          if(log10(xnoborn_cnt).gt.inoborn_cnt)then
@@ -1185,9 +1100,9 @@
       else
          xtot=0d0
          if (mapconfig(0).eq.0) then
-            write (*,*) 'Fatal error in compute_prefactor_n1body,'/
-     &           /' no Born diagrams ',mapconfig
-     &           ,'. Check bornfromreal.inc'
+            write (*,*) 'Fatal error in compute_prefactor_nbody:'/
+     &           /' no Born diagrams ',mapconfig,
+     &           '. Check bornfromreal.inc'
             write (*,*) 'Is fks_singular compiled correctly?'
             stop 1
          endif
@@ -1202,11 +1117,17 @@
          endif
       endif
 
+c In the case there is the special phase-space mapping for resonances,
+c use the Born computed with those as the mapping.
       enhance_real=1.d0
-      if (granny_is_res) then
+      if (granny_is_res .and. imode.eq.2) then
          if (p_born_ev(0,1).gt.0d0) then
             calculatedBorn=.false.
+            pas(0:3,nexternal)=0d0
+            pas(0:3,1:nexternal-1)=p_born_ev(0:3,1:nexternal-1)
+            call set_alphas(pas)
             call sborn(p_born_ev,wgt_c)
+            call set_alphas(p_ev)
             calculatedBorn=.false.
          elseif(p_born_ev(0,1).lt.0d0)then
             if (enhance.ne.0d0) then 
@@ -1246,18 +1167,114 @@
          enhance_real=enhance
       endif
 
+      if (imode.eq.1) then
+         f_b=      f_b      *enhance
+         f_nb=     f_nb     *enhance
+      elseif(imode.eq.2) then
+         f_r=      f_r      *enhance_real
+      elseif(imode.eq.4) then
+         f_MC_S=   f_MC_S   *enhance
+         f_MC_H=   f_MC_H   *enhance
+      elseif(imode.eq.3) then
+         f_s=      f_s      *enhance
+         f_s_MC_S= f_s_MC_S *enhance
+         f_S_MC_H= f_S_MC_H *enhance
+         f_c=      f_c      *enhance
+         f_c_MC_S= f_c_MC_S *enhance
+         f_c_MC_H= f_c_MC_H *enhance
+         f_dc=     f_dc     *enhance
+         f_sc=     f_sc     *enhance
+         f_sc_MC_S=f_sc_MC_S*enhance
+         f_sc_MC_H=f_sc_MC_H*enhance
+         f_dsc(1)= f_dsc(1) *enhance
+         f_dsc(2)= f_dsc(2) *enhance
+         f_dsc(3)= f_dsc(3) *enhance
+         f_dsc(4)= f_dsc(4) *enhance
+         f_dis_d=  f_dis_d  *enhance
+         f_dis_p=  f_dis_p  *enhance
+         f_dis_l=  f_dis_l  *enhance
+      endif
+      call cpu_time(tAfter)
+      tf_nb=tf_nb+(tAfter-tBefore)
+
+      return
+      end
+      
+
+      subroutine compute_prefactors_n1body(vegas_wgt,jac_ev)
+c Compute all relevant prefactors for the real emission and counter
+c terms.
+      implicit none
+      include 'nexternal.inc'
+      include 'run.inc'
+      include 'genps.inc'
+      include 'fks_powers.inc'
+      include 'coupl.inc'
+      include 'timing_variables.inc'
+      double precision vegas_wgt,prefact,prefact_cnt_ssc,prefact_deg
+     $     ,prefact_c,prefact_coll,jac_ev,pi,prefact_cnt_ssc_c
+     $     ,prefact_coll_c,prefact_deg_slxi,prefact_deg_sxi,zero
+      integer i
+      parameter (pi=3.1415926535897932385d0, ZERO=0d0)
+      double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
+     $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
+      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      double precision    xi_i_fks_ev,y_ij_fks_ev
+      double precision    p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
+      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
+      integer            i_fks,j_fks
+      common/fks_indices/i_fks,j_fks
+      double precision   xi_i_fks_cnt(-2:2)
+      common /cxiifkscnt/xi_i_fks_cnt
+      double precision  xinorm_ev
+      common /cxinormev/xinorm_ev
+      double precision  xiimax_ev
+      common /cxiimaxev/xiimax_ev
+      double precision   xiimax_cnt(-2:2)
+      common /cxiimaxcnt/xiimax_cnt
+      double precision   xinorm_cnt(-2:2)
+      common /cxinormcnt/xinorm_cnt
+      double precision    delta_used
+      common /cdelta_used/delta_used
+      double precision    xicut_used
+      common /cxicut_used/xicut_used
+      double precision     xiScut_used,xiBSVcut_used
+      common /cxiScut_used/xiScut_used,xiBSVcut_used
+      double precision        ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
+      double precision         fkssymmetryfactor,fkssymmetryfactorBorn,
+     &                         fkssymmetryfactorDeg
+      integer ngluons,nquarks(-6:6),nphotons
+      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
+     &                  fkssymmetryfactorDeg,ngluons,nquarks,nphotons
+      logical nocntevents
+      common/cnocntevents/nocntevents
+      double precision     f_r,f_s,f_c,f_dc,f_sc,f_dsc(4)
+      common/factor_n1body/f_r,f_s,f_c,f_dc,f_sc,f_dsc
+      double precision           f_s_MC_S,f_s_MC_H,f_c_MC_S,f_c_MC_H
+     $     ,f_sc_MC_S,f_sc_MC_H,f_MC_S,f_MC_H
+      common/factor_n1body_NLOPS/f_s_MC_S,f_s_MC_H,f_c_MC_S,f_c_MC_H
+     $     ,f_sc_MC_S,f_sc_MC_H,f_MC_S,f_MC_H
+      ! prefactors for the DIS scheme
+      double precision prefact_dis_d,prefact_dis_p,prefact_dis_l
+      double precision f_dis_d,f_dis_p,f_dis_l
+      common/factor_dis/f_dis_d,f_dis_p,f_dis_l
+      double precision pmass(nexternal)
+      include 'pmass.inc'
+      call cpu_time(tBefore)
+
 c f_* multiplication factors for real-emission, soft counter, ... etc.       
       prefact=xinorm_ev/xi_i_fks_ev/(1-y_ij_fks_ev)
-      f_r=prefact*jac_ev*enhance_real*fkssymmetryfactor*vegas_wgt
+      f_r=prefact*jac_ev*fkssymmetryfactor*vegas_wgt
       f_MC_S=f_r
       f_MC_H=f_r
       if (.not.nocntevents) then
          prefact_cnt_ssc=xinorm_ev/min(xiimax_ev,xiScut_used)*
      $        log(xicut_used/min(xiimax_ev,xiScut_used))/(1
      $        -y_ij_fks_ev)
-         f_s=(prefact+prefact_cnt_ssc)*jac_cnt(0)*enhance
+         f_s=(prefact+prefact_cnt_ssc)*jac_cnt(0)
      $        *fkssymmetryfactor*vegas_wgt
-         f_s_MC_S=prefact*jac_cnt(0)*enhance
+         f_s_MC_S=prefact*jac_cnt(0)
      $        *fkssymmetryfactor*vegas_wgt
          f_s_MC_H=f_s_MC_S
 
@@ -1268,9 +1285,9 @@
             prefact_coll=xinorm_cnt(1)/xi_i_fks_cnt(1)*log(delta_used
      $           /deltaS)/deltaS
             f_c=(prefact_c+prefact_coll)*jac_cnt(1)
-     $           *enhance*fkssymmetryfactor*vegas_wgt
+     $           *fkssymmetryfactor*vegas_wgt
             f_c_MC_S=prefact_c*jac_cnt(1)
-     $           *enhance*fkssymmetryfactor*vegas_wgt
+     $           *fkssymmetryfactor*vegas_wgt
             f_c_MC_H=f_c_MC_S
 
             call set_cms_stuff(1)
@@ -1281,13 +1298,13 @@
             prefact_coll_c=xinorm_cnt(1)/min(xiimax_cnt(1),xiScut_used)
      $           *log(xicut_used/min(xiimax_cnt(1),xiScut_used))
      $           *log(delta_used/deltaS)/deltaS
-            f_dc=jac_cnt(1)*prefact_deg/(shat/(32*pi**2))*enhance
+            f_dc=jac_cnt(1)*prefact_deg/(shat/(32*pi**2))
      $           *fkssymmetryfactorDeg*vegas_wgt
             f_sc=(prefact_c+prefact_coll+prefact_cnt_ssc_c
-     &           +prefact_coll_c)*jac_cnt(2)*enhance
+     &           +prefact_coll_c)*jac_cnt(2)
      &           *fkssymmetryfactorDeg*vegas_wgt
             f_sc_MC_S=prefact_c*jac_cnt(2)
-     $           *enhance*fkssymmetryfactor*vegas_wgt
+     $           *fkssymmetryfactor*vegas_wgt
             f_sc_MC_H=f_sc_MC_S
 
             call set_cms_stuff(2)
@@ -1298,24 +1315,24 @@
      &           ,xiScut_used)*( log(xicut_used)**2
      &           -log(min(xiimax_cnt(1),xiScut_used))**2 )*1/(2.d0
      &           *deltaS)
-            f_dsc(1)=prefact_deg*jac_cnt(2)/(shat/(32*pi**2))*enhance
+            f_dsc(1)=prefact_deg*jac_cnt(2)/(shat/(32*pi**2))
      &           *fkssymmetryfactorDeg*vegas_wgt
             f_dsc(2)=prefact_deg_sxi*jac_cnt(2)/(shat/(32*pi**2))
-     &           *enhance*fkssymmetryfactorDeg*vegas_wgt
+     &           *fkssymmetryfactorDeg*vegas_wgt
             f_dsc(3)=prefact_deg_slxi*jac_cnt(2)/(shat/(32*pi**2))
-     &           *enhance*fkssymmetryfactorDeg*vegas_wgt
+     &           *fkssymmetryfactorDeg*vegas_wgt
             f_dsc(4)=( prefact_deg+prefact_deg_sxi )*jac_cnt(2)/(shat
-     &           /(32*pi**2))*enhance*fkssymmetryfactorDeg
+     &           /(32*pi**2))*fkssymmetryfactorDeg
      &           *vegas_wgt
             ! prefactor for the DIS scheme
             prefact_dis_d=xinorm_cnt(1)/xiScut_used/deltaS
-            f_dis_d=prefact_dis_d*jac_cnt(2)/(shat/(32*pi**2))*enhance
+            f_dis_d=prefact_dis_d*jac_cnt(2)/(shat/(32*pi**2))
      &           *fkssymmetryfactorDeg*vegas_wgt
             prefact_dis_p=xinorm_cnt(1)*dlog(xiScut_used)/xiScut_used/deltaS
-            f_dis_p=prefact_dis_p*jac_cnt(2)/(shat/(32*pi**2))*enhance
+            f_dis_p=prefact_dis_p*jac_cnt(2)/(shat/(32*pi**2))
      &           *fkssymmetryfactorDeg*vegas_wgt
             prefact_dis_l=xinorm_cnt(1)*dlog(xiScut_used)**2/2d0/xiScut_used/deltaS
-            f_dis_l=prefact_dis_l*jac_cnt(2)/(shat/(32*pi**2))*enhance
+            f_dis_l=prefact_dis_l*jac_cnt(2)/(shat/(32*pi**2))
      &           *fkssymmetryfactorDeg*vegas_wgt
          else
             f_c=0d0
@@ -1439,6 +1456,7 @@
 c        contribution
       use weight_lines
       use extra_weights
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
@@ -1446,7 +1464,6 @@
       include 'coupl.inc'
       include 'fks_info.inc'
       include 'q_es.inc'
-      include 'FKSParams.inc'
       include 'orders.inc'
       integer type,i,j
       logical foundIt,foundOrders
@@ -1633,6 +1650,8 @@
 c or to fill histograms.
       use weight_lines
       use extra_weights
+      use mint_module
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
@@ -1640,10 +1659,8 @@
       include 'timing_variables.inc'
       include 'genps.inc'
       include 'orders.inc'
-      include 'mint.inc'
-      include 'FKSParams.inc'
       integer orders(nsplitorders)
-      integer i,j,k,iamp
+      integer i,j,k,iamp,icontr_orig
       logical virt_found
       double precision xlum,dlum,pi,mu2_r,mu2_f,mu2_q,rwgt_muR_dep_fac
      $     ,wgt_wo_pdf,conv
@@ -1652,9 +1669,6 @@
       external dlum
       integer              nFKSprocess
       common/c_nFKSprocess/nFKSprocess
-      double precision           virt_wgt_mint(0:n_ave_virt),
-     &                           born_wgt_mint(0:n_ave_virt)
-      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
       INTEGER              IPROC
       DOUBLE PRECISION PD(0:MAXPROC)
       COMMON /SUBPROC/ PD, IPROC
@@ -1662,6 +1676,9 @@
       call cpu_time(tBefore)
       if (icontr.eq.0) return
       virt_found=.false.
+c number of contributions before they are (possibly) increased through a
+c call to separate_flavour_config().
+      icontr_orig=icontr
       i=0
       do while (i.lt.icontr)
          i=i+1
@@ -1678,8 +1695,12 @@
 c iwgt=1 is the central value (i.e. no scale/PDF reweighting).
          iwgt=1
          call weight_lines_allocated(nexternal,max_contr,iwgt,iproc)
-c set_pdg_codes fills the niproc, parton_iproc, parton_pdg and parton_pdg_uborn
-         call set_pdg_codes(iproc,pd,nFKSprocess,i)
+c set_pdg_codes fills the niproc, parton_iproc, parton_pdg and
+c parton_pdg_uborn [Do only for the contributions that were already
+c available as part of the input -- NOT the ones that are created
+c through the call to separate_flavour_config(), since that will
+c overwrite the relevant information.]
+         if (i.le.icontr_orig) call set_pdg_codes(iproc,pd,nFKSprocess,i)
          if (separate_flavour_configs .and. ipr(i).eq.0) then
             call separate_flavour_config(i) ! this increases icontr
          endif
@@ -1733,18 +1754,15 @@
 c include the weight also in the 'wgt' array that contain the
 c coefficients for PDF and scale computations. 
       use weight_lines
+      use mint_module
       implicit none
       include 'orders.inc'
-      include 'mint.inc'
       integer orders(nsplitorders)
       integer i,j,iamp
       logical virt_found
       double precision bias
       character*7 event_norm
       common /event_normalisation/event_norm
-      double precision           virt_wgt_mint(0:n_ave_virt),
-     &                           born_wgt_mint(0:n_ave_virt)
-      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
 c Set the bias_wgt to 1 in case we do not have to do any biassing
       if (event_norm(1:4).ne.'bias') then
          do i=1,icontr
@@ -1993,11 +2011,11 @@
 c wgts() array to include the weights.
       use weight_lines
       use extra_weights
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'timing_variables.inc'
-      include 'FKSParams.inc'
       include 'genps.inc'
       integer i,kr,kf,iwgt_save,dd
       double precision xlum(maxscales),dlum,pi,mu2_r(maxscales),c_mu2_r
@@ -2076,11 +2094,11 @@
 c computations (ickkw.eq.-1).
       use weight_lines
       use extra_weights
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'timing_variables.inc'
-      include 'FKSParams.inc'
       include 'genps.inc'
       integer i,ks,kh,iwgt_save
       double precision xlum(maxscales),dlum,pi,mu2_r(maxscales)
@@ -2182,11 +2200,11 @@
 c wgts() array to include the weights.
       use weight_lines
       use extra_weights
+      use FKSParams
       implicit none
       include 'nexternal.inc'
       include 'run.inc'
       include 'timing_variables.inc'
-      include 'FKSParams.inc'
       include 'genps.inc'
       integer n,izero,i,nn
       parameter (izero=0)
@@ -2569,15 +2587,12 @@
       subroutine fill_mint_function(f)
 c Fills the function that is returned to the MINT integrator
       use weight_lines
+      use mint_module
       implicit none
       include 'nexternal.inc'
-      include 'mint.inc'
       include 'orders.inc'
       integer i,iamp,ithree,isix
       double precision f(nintegrals),sigint
-      double precision           virt_wgt_mint(0:n_ave_virt),
-     &                           born_wgt_mint(0:n_ave_virt)
-      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
       double precision virtual_over_born
       common /c_vob/   virtual_over_born
       sigint=0d0
@@ -2674,6 +2689,7 @@
 c configuration, while for the S-events also contributions from the
 c various FKS configurations can be summed together.
       use weight_lines
+      use mint_module
       implicit none
       include 'nexternal.inc'
       include 'genps.inc'
@@ -2686,9 +2702,6 @@
       integer iproc_save(fks_configs),eto(maxproc,fks_configs),
      &     etoi(maxproc,fks_configs),maxproc_found
       common/cproc_combination/iproc_save,eto,etoi,maxproc_found
-      logical               only_virt
-      integer         imode
-      common /c_imode/imode,only_virt
       call cpu_time(tBefore)
       if (icontr.eq.0) return
 c Find the contribution to sum all the S-event ones. This should be one
@@ -2835,21 +2848,15 @@
 c Fills the function that is returned to the MINT integrator. Depending
 c on the imode we should or should not include the virtual corrections.
       use weight_lines
+      use mint_module
       implicit none
       include 'nexternal.inc'
-      include 'mint.inc'
       include 'orders.inc'
       integer i,j,ict,iamp,ithree,isix
       double precision f(nintegrals),sigint,sigint1,sigint_ABS
      $     ,n1body_wgt,tmp_wgt,max_weight
-      double precision           virt_wgt_mint(0:n_ave_virt),
-     &                           born_wgt_mint(0:n_ave_virt)
-      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
       double precision virtual_over_born
       common /c_vob/   virtual_over_born
-      logical               only_virt
-      integer         imode
-      common /c_imode/imode,only_virt
       sigint=0d0
       sigint1=0d0
       sigint_ABS=0d0
@@ -3621,10 +3628,12 @@
      &                                            pjet,njet,jet)
                do i=1,NN
                   di(i)=sqrt(amcatnlo_fastjetdmergemax(i-1))
-                  if (i.gt.1.and.di(i).gt.di(i-1))then
-                     write(*,*)'Error in set_shower_scale_noshape'
-                     write(*,*)NN,i,di(i),di(i-1)
-                     stop
+                  if (i.gt.1)then
+                     if(di(i).gt.di(i-1))then
+                        write(*,*)'Error in set_shower_scale_noshape'
+                        write(*,*)NN,i,di(i),di(i-1)
+                        stop
+                     endif
                   endif
                enddo
                if(iSH.eq.1)shower_S_scale(iFKS)=di(NN)
@@ -5470,6 +5479,7 @@
 
       subroutine bornsoftvirtual(p,bsv_wgt,virt_wgt,born_wgt)
       use extra_weights
+      use mint_module
       implicit none
       include "genps.inc"
       include 'nexternal.inc'
@@ -5483,7 +5493,6 @@
       common /c_charges/particle_charge
       include "run.inc"
       include "fks_powers.inc"
-      include "mint.inc"
       double precision p(0:3,nexternal),bsv_wgt,born_wgt,avv_wgt
       double precision pp(0:3,nexternal)
       
@@ -5542,9 +5551,6 @@
       integer iminmax
       common/cExceptPSpoint/iminmax,ExceptPSpoint
 
-      double precision average_virtual(0:n_ave_virt,maxchannels)
-     $     ,virtual_fraction(maxchannels)
-      common/c_avg_virt/average_virtual,virtual_fraction
       double precision virtual_over_born
       common/c_vob/virtual_over_born
 
@@ -5871,12 +5877,22 @@
             tOLP=tOLP+(tAfter-tBefore)
             virtual_over_born=virt_wgt/born_wgt
             if (ickkw.ne.-1) then
-               virt_wgt=virt_wgt-average_virtual(0,ichan)*born_wgt
+               if (use_poly_virtual) then
+                  virt_wgt=virt_wgt-polyfit(0)*born_wgt
+               else
+                  virt_wgt=virt_wgt-average_virtual(0,ichan)*born_wgt
+               endif
                do iamp=1,amp_split_size
                   if (amp_split_virt(iamp).eq.0d0) cycle
+                  if (use_poly_virtual) then
+                     amp_split_virt(iamp)=amp_split_virt(iamp)-
+     $                    polyfit(iamp)
+     $                    *amp_split_born_for_virt(iamp)
+               else
                   amp_split_virt(iamp)=amp_split_virt(iamp)-
      $                 average_virtual(iamp,ichan)
      $                 *amp_split_born_for_virt(iamp)
+                  endif
                enddo
             endif
             if (abrv.ne.'virt') then
@@ -5896,12 +5912,21 @@
      $        amp_split_virt_save(1:amp_split_size)
       endif
       if (abrv(1:4).ne.'virt' .and. ickkw.ne.-1) then
-         avv_wgt=average_virtual(0,ichan)*born_wgt
-         do iamp=1, amp_split_size
-            if (amp_split_born_for_virt(iamp).eq.0d0) cycle
-            amp_split_avv(iamp)= average_virtual(iamp,ichan)
-     $           *amp_split_born_for_virt(iamp)
-         enddo
+         if (use_poly_virtual) then
+            avv_wgt=polyfit(0)*born_wgt
+            do iamp=1, amp_split_size
+               if (amp_split_born_for_virt(iamp).eq.0d0) cycle
+               amp_split_avv(iamp)= polyfit(iamp)
+     $              *amp_split_born_for_virt(iamp)
+            enddo
+         else
+            avv_wgt=average_virtual(0,ichan)*born_wgt
+            do iamp=1, amp_split_size
+               if (amp_split_born_for_virt(iamp).eq.0d0) cycle
+               amp_split_avv(iamp)= average_virtual(iamp,ichan)
+     $              *amp_split_born_for_virt(iamp)
+            enddo
+         endif
       endif
 
 c eq.(MadFKS.C.13)
@@ -6512,9 +6537,8 @@
       subroutine setfksfactor(match_to_shower)
       use weight_lines
       use extra_weights
+      use mint_module
       implicit none
-
-      include 'mint.inc'
       
       double precision CA,CF, PI
       parameter (CA=3d0,CF=4d0/3d0)
@@ -6529,7 +6553,7 @@
       logical softtest,colltest
       common/sctests/softtest,colltest
 
-      integer config_fks,i,j,fac1,fac2
+      integer config_fks,i,j,fac1,fac2,kchan
       logical match_to_shower
 
       double precision fkssymmetryfactor,fkssymmetryfactorBorn,
@@ -6917,11 +6941,11 @@
 c the grandmother corresponds (igranny) as well as the aunt (iaunt).
 c This information can be used to improve the phase-space
 c parametrisation.
+      use mint_module
       implicit none
       include 'genps.inc'
       include 'nexternal.inc'
       include 'nFKSconfigs.inc'
-      include 'mint.inc'
 c arguments
       integer nFKSprocess,iconf
       double precision mass_min(-nexternal:nexternal)
@@ -7110,3 +7134,21 @@
       endif
       return
       end
+
+
+      function ran2()
+!     Wrapper for the random numbers; needed for the NLO stuff
+      use mint_module
+      implicit none
+      double precision ran2,x,a,b
+      integer ii,jconfig
+      a=0d0                     ! min allowed value for x
+      b=1d0                     ! max allowed value for x
+      ii=0                      ! dummy argument of ntuple
+      jconfig=iconfig           ! integration channel (for off-set)
+      call ntuple(x,a,b,ii,jconfig)
+      ran2=x
+      return
+      end
+
+      

=== removed file 'Template/NLO/SubProcesses/fks_singular_paolo.f'
=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
--- Template/NLO/SubProcesses/genps_fks.f	2018-07-13 10:16:34 +0000
+++ Template/NLO/SubProcesses/genps_fks.f	2020-05-23 19:29:41 +0000
@@ -95,8 +95,8 @@
       double precision p(0:3,nexternal)
       integer itree(2,-max_branch:-1),iconf
       common /to_itree/itree,iconf
-      integer ndim
-      common/tosigint/ndim
+      integer nndim
+      common/tosigint/nndim
       double precision qmass_common(-nexternal:0),qwidth_common(
      &     -nexternal:0)
       common /c_qmass_qwidth/qmass_common,qwidth_common
@@ -110,7 +110,7 @@
          virtgranny=0d0
          return
       endif
-      call generate_momenta_conf(.true.,ndim,jac,xvar,granny_m2_red
+      call generate_momenta_conf(.true.,nndim,jac,xvar,granny_m2_red
      &     ,rat_xi_orig,itree,qmass_common,qwidth_common,p)
       if (jac.gt.0d0) then
          do j=0,3
@@ -152,8 +152,8 @@
       double precision p(0:3,nexternal)
       integer itree(2,-max_branch:-1),iconf
       common /to_itree/itree,iconf
-      integer ndim
-      common/tosigint/ndim
+      integer nndim
+      common/tosigint/nndim
       double precision qmass_common(-nexternal:0),qwidth_common(
      &     -nexternal:0)
       common /c_qmass_qwidth/qmass_common,qwidth_common
@@ -177,7 +177,7 @@
          virtgranny_red=0d0
          return
       endif
-      call generate_momenta_conf(.true.,ndim,jac,xvar,granny_m2_red
+      call generate_momenta_conf(.true.,nndim,jac,xvar,granny_m2_red
      &     ,rat_xi_orig,itree,qmass_common,qwidth_common,p)
       if (jac.gt.0d0) then
          do j=1,3
@@ -282,13 +282,14 @@
 
 
       
-      subroutine generate_momenta_conf_wrapper(ndim,jac,x,itree,qmass
+      subroutine generate_momenta_conf_wrapper(nndim,jac,x,itree,qmass
      $     ,qwidth,p)
+      use mint_module
       implicit none
       include 'nexternal.inc'
       include 'genps.inc'
       include 'nFKSconfigs.inc'
-      integer ndim
+      integer nndim
       double precision jac,x(99),p(0:3,nexternal)
       integer itree(2,-max_branch:-1),i,j
       double precision qmass(-nexternal:0),qwidth(-nexternal:0),del1
@@ -367,8 +368,6 @@
       common /to_abrv/ abrv
       character*10 shower_mc
       common /cMonteCarloType/shower_mc
-      logical              fixed_order,nlo_ps
-      common /c_fnlo_nlops/fixed_order,nlo_ps
 c
       write_granny(nFKSprocess)=.true.
       which_is_granny(nFKSprocess)=0
@@ -405,7 +404,7 @@
             only_event_phsp=.true.
             skip_event_phsp=.false.
             input_granny_m2=.false.
-            call generate_momenta_conf(input_granny_m2,ndim,jac,x
+            call generate_momenta_conf(input_granny_m2,nndim,jac,x
      $           ,granny_m2_red,rat_xi,itree,qmass,qwidth,p)
             if (nint(jac).eq.-222) return ! numerical inaccuracy: not
                                           ! even Born momenta generated.
@@ -540,7 +539,7 @@
 c a shower scale. The kinematics won't be used, because below jac=-222
 c and p(0,1)=-99d0.
             skip_event_phsp=.false. 
-            call generate_momenta_conf(input_granny_m2,ndim,jac,x
+            call generate_momenta_conf(input_granny_m2,nndim,jac,x
      $           ,granny_m2_red,rat_xi,itree,qmass,qwidth,p)
             jac=-222
             p(0,1)=-99d0
@@ -554,7 +553,7 @@
             do i=-1,1
                granny_m2_red(i)=-99d99
             enddo
-            call generate_momenta_conf(input_granny_m2,ndim,jac,x
+            call generate_momenta_conf(input_granny_m2,nndim,jac,x
      $           ,granny_m2_red,rat_xi,itree,qmass,qwidth,p)
 c In this case, we should not write the grandmother in the event file,
 c because the shower should not keep its inv. mass fixed.
@@ -592,7 +591,7 @@
             skip_event_phsp=.false.
             only_event_phsp=.true.
             granny_m2_red(0)=xmbe2inv
-            call generate_momenta_conf(input_granny_m2,ndim,jac,x
+            call generate_momenta_conf(input_granny_m2,nndim,jac,x
      $           ,granny_m2_red,rat_xi_orig,itree,qmass,qwidth,p)
 c multiply event jacobian by the numerically computed jacobian for the
 c derivative
@@ -602,7 +601,7 @@
             input_granny_m2=.false.
             only_event_phsp=.false.
             skip_event_phsp=.true.
-            call generate_momenta_conf(input_granny_m2,ndim,jac,x
+            call generate_momenta_conf(input_granny_m2,nndim,jac,x
      $           ,granny_m2_red,rat_xi,itree,qmass,qwidth,p)
          endif
          jac=jac*MC_sum_factor
@@ -610,7 +609,7 @@
          skip_event_phsp=.false.
          only_event_phsp =.false.
          input_granny_m2=.false.
-         call generate_momenta_conf(input_granny_m2,ndim,jac,x
+         call generate_momenta_conf(input_granny_m2,nndim,jac,x
      $        ,granny_m2_red,rat_xi,itree,qmass,qwidth,p)
       endif
       end
@@ -701,9 +700,6 @@
 c Masses of particles. Should be filled in setcuts.f
       double precision pmass(nexternal)
       common /to_mass/pmass
-c     For MINT:
-      integer ifold_energy,ifold_phi,ifold_yij
-      common /cifoldnumbers/ifold_energy,ifold_phi,ifold_yij
 c local
       integer i,j,nbranch,ns_channel,nt_channel,ionebody
      &     ,fksconfiguration,icountevts,imother,ixEi,ixyij,ixpi,isolsign
@@ -1000,14 +996,10 @@
       if (.not.only_event_phsp) xij_aor=(0.d0,0.d0)
 c
 c These will correspond to the vegas x's for the FKS variables xi_i,
-c y_ij and phi_i
+c y_ij and phi_i (changing this also requires changing folding parameters)
       ixEi=ndim-2
       ixyij=ndim-1
       ixpi=ndim
-c Set up the MINT folding:
-      ifold_energy=ixEi
-      ifold_phi=ixpi
-      ifold_yij=ixyij
 c
       imother=min(j_fks,i_fks)
       m_j_fks=pmass(j_fks)
@@ -2147,6 +2139,7 @@
       xp(3,2)=-xp(0,2)
 c
       xp(0,i_fks)=E_i_fks*(chy_lbst-shy_lbst*yijdir)
+      p_i_fks(0)=p_i_fks(0)*(chy_lbst-shy_lbst*yijdir)
       xpifksred(1)=sinth_i_fks*cosphi_i_fks    
       xpifksred(2)=sinth_i_fks*sinphi_i_fks    
       xpifksred(3)=chy_lbst*yijdir-shy_lbst

=== modified file 'Template/NLO/SubProcesses/handling_lhe_events.f'
--- Template/NLO/SubProcesses/handling_lhe_events.f	2017-06-15 19:31:40 +0000
+++ Template/NLO/SubProcesses/handling_lhe_events.f	2020-05-23 19:29:41 +0000
@@ -602,7 +602,7 @@
       endif
       write(ifile,'(a)')
      # '  </init>'
- 501  format(2(1x,i6),2(1x,e14.8),2(1x,i2),2(1x,i6),1x,i2,1x,i3)
+ 501  format(2(1x,i6),2(1x,e14.8),2(1x,i2),2(1x,i8),1x,i2,1x,i3)
  502  format(3(1x,e14.8),1x,i6)
 c
       return
@@ -665,8 +665,10 @@
 c     if event_id is zero or positive (that means that there was a call
 c     to write_lhef_header_banner) update it and write it
 c RF: don't use the event_id:
+      logical do_rwgt
       event_id = -99
 c
+      do_rwgt = do_rwgt_scale .or. do_rwgt_pdf .or. store_rwgt_info
       if (event_id.ge.0) then
          event_id=event_id+1
          if (event_id.le.9) then
@@ -720,7 +722,7 @@
      #                  PUP(1,I),PUP(2,I),PUP(3,I),PUP(4,I),PUP(5,I),
      #                  VTIMUP(I),SPINUP(I)
       enddo
-      if(buff(1:1).eq.'#' .and. (do_rwgt_scale .or. do_rwgt_pdf .or.
+      if(buff(1:1).eq.'#' .and. (do_rwgt .or.
      &     jwgtinfo.lt.0)) then
          write(ifile,'(a)') buff(1:len_trim(buff))
          read(buff,*)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,

=== modified file 'Template/NLO/SubProcesses/madfks_plot.f'
--- Template/NLO/SubProcesses/madfks_plot.f	2017-07-31 08:48:20 +0000
+++ Template/NLO/SubProcesses/madfks_plot.f	2020-05-23 19:29:41 +0000
@@ -148,14 +148,11 @@
 
       subroutine topout
       use extra_weights
+      use mint_module
       implicit none
       include "nexternal.inc"
       include 'run.inc'
       integer ii,jj,n,kk,nn
-      logical usexinteg,mint
-      common/cusexinteg/usexinteg,mint
-      integer itmax,ncall
-      common/citmax/itmax,ncall
       logical useitmax
       common/cuseitmax/useitmax
       real*8 xnorm
@@ -166,16 +163,10 @@
       common /for_applgrid/ iappl
       include "appl_common.inc"
 c
-      if(usexinteg.and..not.mint) then
-         xnorm=1.d0/float(itmax)
-      elseif(mint) then
-         xnorm=1.d0/float(ncall)
-      else
-         xnorm=1d0
-      endif
+      xnorm=1.d0/float(ncalls0)
       if(useitmax)xnorm=xnorm/float(itmax)
 c Normalization factor for the APPLgrid grids
-      if(iappl.ne.0) appl_norm_histo = 1d0 / dble(ncall*itmax)
+      if(iappl.ne.0) appl_norm_histo = 1d0 / dble(ncalls0*itmax)
       call analysis_end(xnorm)
 c Write the accumulated results to a file
       open (unit=34,file='scale_pdf_dependence.dat',status='unknown')

=== modified file 'Template/NLO/SubProcesses/makefile_fks_dir'
--- Template/NLO/SubProcesses/makefile_fks_dir	2017-09-28 15:13:09 +0000
+++ Template/NLO/SubProcesses/makefile_fks_dir	2020-05-23 19:29:41 +0000
@@ -25,42 +25,37 @@
 $(LIBDIR)libpdf.a $(LIBDIR)libcernlib.a
 
 # Files for all executables
-FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) 			\
-	 $(patsubst %.f,%.o,$(wildcard matrix_*.f)) 			\
-	 real_me_chooser.o chooser_functions.o				\
-	 genps_fks.o setcuts.o setscales.o veto_xsec.o		 	\
-	 $(patsubst %.f,%.o,$(wildcard b_sf_???.f)) 			\
-	 born.o sborn_sf.o extra_cnt_wrapper.o				\
-	 $(patsubst %.f,%.o,$(wildcard born_cnt_*.f))			\
-	 fks_Sij.o $(fastjetfortran_madfks) fks_singular.o		\
-	 montecarlocounter.o reweight_xsec.o boostwdir2.o       	\
-	 initcluster.o cluster.o splitorders_stuff.o 			\
-	 reweight.o get_color.o FKSParamReader.o iproc_map.o		\
-	 MC_integer.o $(reweight_xsec_events_pdf_dummy)			\
-	 $(applgrid_interface) weight_lines.o
+FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) $(patsubst	\
+	%.f,%.o,$(wildcard matrix_*.f)) FKSParams.o real_me_chooser.o	\
+	chooser_functions.o genps_fks.o setcuts.o setscales.o		\
+	veto_xsec.o $(patsubst %.f,%.o,$(wildcard b_sf_???.f))		\
+	born.o sborn_sf.o extra_cnt_wrapper.o $(patsubst		\
+	%.f,%.o,$(wildcard born_cnt_*.f)) fks_Sij.o			\
+	$(fastjetfortran_madfks) fks_singular.o montecarlocounter.o	\
+	reweight_xsec.o boostwdir2.o initcluster.o cluster.o		\
+	splitorders_stuff.o reweight.o get_color.o 			\
+	iproc_map.o MC_integer.o $(reweight_xsec_events_pdf_dummy)	\
+	$(applgrid_interface) weight_lines.o mint_module.o polfit.o
 
 # Files needed for mintFO & mintMC
 RUN= $(FO_ANALYSE) $(FILES) cuts.o pythia_unlops.o recluster.o 		\
 	fill_MC_mshell.o born_hel.o open_output_files.o			\
-	add_write_info.o BinothLHA.o madfks_plot.o
+	add_write_info.o madfks_plot.o
 
 # Files for mintFO
-MINTFO= $(RUN) driver_mintFO.o mint-integrator2.o handling_lhe_events.o
+MINTFO= $(RUN) driver_mintFO.o handling_lhe_events.o BinothLHA.o
 
 # Files for mintMC
-MINTMC= $(RUN) mint-integrator2.o driver_mintMC.o write_event.o		\
-	handling_lhe_events.o
+MINTMC= $(RUN) driver_mintMC.o handling_lhe_events.o write_event.o BinothLHA.o
 
 # Files for check_poles
-POLES= $(FILES) cuts.o pythia_unlops.o recluster.o check_poles.o	\
-       BinothLHA.o born_hel.o
+POLES= $(RUN) check_poles.o BinothLHA.o
 
 # Files for tests
-TEST= $(FILES) test_soft_col_limits.o BinothLHADummy.o cuts.o		\
-      pythia_unlops.o recluster.o
+TEST= $(RUN) test_soft_col_limits.o BinothLHADummy.o
 
 # Files for gensym
-SYM= $(FILES) BinothLHADummy.o symmetry_fks_v3.o write_ajob.o
+SYM= $(RUN) BinothLHADummy.o symmetry_fks_v3.o write_ajob.o
 
 # Files for event reweighting to get e.g. scale and PDF uncertainties
 RWGFILES=reweight_xsec_events.o handling_lhe_events.o			\
@@ -73,6 +68,8 @@
 
 %.o: %.f
 	$(FC) $(FFLAGS) -c -I. -I$(LIBDIR) $<
+%.o: %.f90
+	$(FC) $(FFLAGS) -c -I. -I$(LIBDIR) $<
 %.o: $(ANADIR)/%.f
 	$(FC) $(FFLAGS) -c -I. -I$(ANADIR) $(FO_INCLUDEPATHS) $<
 %.o: %.for
@@ -121,9 +118,13 @@
 	for dir in `ls -d V*`; do cd $$dir; make; cd ../; done
 
 reweight_xsec_events.o: weight_lines.o
-fks_singular.o: weight_lines.o
-driver_mintMC.o: weight_lines.o
-driver_mintFO.o: weight_lines.o
+fks_singular.o: weight_lines.o mint_module.o
+driver_mintMC.o: weight_lines.o mint_module.o
+driver_mintFO.o: weight_lines.o mint_module.o
+setcuts.o: mint_module.o
+genps_fks.o: mint_module.o
+symmetry_fks_v3.o: mint_module.o
+test_soft_col_limits.o: mint_module.o
 
 clean:
 	@rm -f *.o write_ajob_basic reweight_xsec_events madevent_* check_poles\

=== removed file 'Template/NLO/SubProcesses/mint-integrator2.f'
=== removed file 'Template/NLO/SubProcesses/mint.inc'
=== added file 'Template/NLO/SubProcesses/mint_module.f90'
--- Template/NLO/SubProcesses/mint_module.f90	1970-01-01 00:00:00 +0000
+++ Template/NLO/SubProcesses/mint_module.f90	2020-05-23 19:29:41 +0000
@@ -0,0 +1,1705 @@
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! MINT Integrator Package
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! Original version by Paolo Nason (for POWHEG (BOX))
+! Modified by Rikkert Frederix (for MadGraph5_aMC@NLO)
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!      subroutine mint(fun,ndim,ncalls0,itmax,imode,
+! ndim=number of dimensions
+! ncalls0=# of calls per iteration
+! itmax =# of iterations
+! fun(xx,www,ifirst): returns the function to be integrated multiplied by www;
+!                     xx(1:ndim) are the variables of integration
+!                     ifirst=0: normal behaviour
+! imode: integer flag
+!
+! imode=-1:
+! same as imode=0 as far as this routine is concerned, except for the
+! fact that a grid is read at the beginning (rather than initialized).
+! The return value of imode will be zero.
+!
+! imode=0:
+! When called with imode=0 the routine integrates the absolute value of
+! the function and sets up a grid xgrid(0:50,ndim) such that in each
+! ndim-1 dimensional slice (i.e. xgrid(m-1,n)<xx(n)<xgrid(m,n)) the
+! contribution of the integral is the same the array xgrid is setup at
+! this stage; ans and err are the integral and its error
+!
+! imode=1 (in fact #0)
+! When called with imode=1, the routine performs the integral of the
+! function fun using the grid xgrid. If some number in the array ifold,
+! (say, ifold(n)) is different from 1, it must be a divisor of 50, and
+! the 50 intervals xgrid(0:50,n) are grouped into ifold(n) groups, each
+! group containing 50/ifold(n) nearby intervals. For example, if
+! ifold(1)=5, the 50 intervals for the first dimension are divided in 5
+! groups of 10. The integral is then performed by folding on top of each
+! other these 5 groups. Suppose, for example, that we choose a random
+! point in xx(1) = xgrid(2,1)+x*(xgrid(3,1)-xgrid(2,1)), in the group of
+! the first 5 interval.  we sum the contribution of this point to the
+! contributions of points
+! xgrid(2+m*10,1)+x*(xgrid(3+m*10,1)-xgrid(2+m*10,1)), with m=1,...,4.
+! In the sequence of calls to the function fun, the call for the first
+! point is performed with ifirst=0, and that for all subsequent points
+! with ifirst=1, so that the function can avoid to compute quantities
+! that only depend upon dimensions that have ifold=1, and do not change
+! in each group of folded call. The values returned by fun in a sequence
+! of folded calls with ifirst=0 and ifirst=1 are not used. The function
+! itself must accumulate the values, and must return them when called
+! with ifirst=2.
+! 
+! Added the posibility to keep track of more than one integral:
+!
+! nintegrals=1 : the function that is used to update the grids. This is
+! the ABS cross section. If imode.eq.1, this does not contain the
+! virtual corrections because for them a separate maximum is kept using (5).
+! nintegrals=2 : the actual cross section. This includes virtual corrections.
+! nintegrals=3 : the cross section from the M_Virt/M_Born ratio alone:
+! this defines the average virtual that is added to each phase-space
+! point
+! nintegrals=4 : the cross section of the actual virtual minus the
+! average virtual. This is used to determine the fraction of phase-space
+! points for which we include the virtual.
+! nintegrals=5 : abs of 3
+! nintegrals=6 : born
+! nintegrals>6 : virtual and born order by order
+!
+
+module mint_module
+  use FKSParams ! contains use_poly_virtual
+  implicit none
+  integer, parameter, private :: nintervals=32    ! max number of intervals in the integration grids
+  integer, parameter, public  :: ndimmax=60       ! max number of dimensions of the integral
+  integer, parameter, public  :: n_ave_virt=10    ! max number of grids to set up to approx virtual
+  integer, parameter, public  :: nintegrals=26    ! number of integrals to keep track of
+  integer, parameter, private :: nintervals_virt=8! max number of intervals in the grids for the approx virtual
+  integer, parameter, private :: min_inter=4      ! minimal number of intervals
+  integer, parameter, private :: min_it0=4        ! minimal number of iterations in the mint step 0 phase
+  integer, parameter, private :: min_it1=5        ! minimal number of iterations in the mint step 1 phase
+  integer, parameter, private :: max_points=100000! maximum number of points to trow per iteration if not enough non-zero points can be found.
+  integer, parameter, public  :: maxchannels=20 ! set as least as large as in amcatnlo_run_interface
+  ! Note that the number of intervals in the integration grids, 'nintervals', cannot be arbitrarily large.
+  ! It should be equal to
+  !     nintervals = min_inter * 2^n,
+  ! where 'n' is an integer smaller than or equal to min(min_it0,min_it1).
+  !
+  ! The number of intergrals should be equal to
+  !     nintegrals=6+2*n_ave_virt
+  !
+
+! public variables 
+  integer, public :: ncalls0,ndim,itmax,imode,n_ord_virt,nchans,iconfig,ichan,ifold_energy,ifold_yij,ifold_phi
+  integer, dimension(ndimmax), public :: ifold
+  integer, dimension(maxchannels), public :: iconfigs
+  double precision, public :: accuracy,min_virt_fraction_mint,wgt_mult
+  double precision, dimension(0:n_ave_virt,maxchannels), public :: average_virtual
+  double precision, dimension(0:n_ave_virt), public :: virt_wgt_mint,born_wgt_mint,polyfit
+  double precision, dimension(maxchannels), public :: virtual_fraction
+  double precision, dimension(nintegrals,0:maxchannels), public :: ans,unc
+  logical :: only_virt,new_point,pass_cuts_check
+
+! private variables
+  character(len=13), parameter, dimension(nintegrals), private :: title=(/ &