Merge lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny

[Rikkert Frederix]

Hi guys,

Before getting super excited and doing lots of test for the 4FS, we should release this code first...

Note: I proposed for merging with the FKS_EW_granny branch -- I don't know if this is exactly the branch relevant to the 3.0.0 release.

Cheers,
Rik

[marco zaro (marco-zaro)]

have we enabled the NLO+PS for QCD only?

[Rikkert Frederix (frederix)]

It's not enabled. We should. Could you please do this?

[lp:~maddevelopers/mg5amcnlo/3.0.1]

942. By marco zaro on 2018-10-04

fixes in amcatnlo_run_interface.py to enable NLOPS for QCD-only
not fully working yet

[marco zaro (marco-zaro)]

I have pushed a fix, however, it is not complete yet (and the behaviour is quite strange)
for example, I cannot set fixed_order=OFF via the switches.
However, ./bin/generate_events aMC@NLO works as expected (and herwig is compiled and run too), however on the screen the printout one gets is
launch aMC@NLO
The following switches determine which programs are run:
/================== Description ===================|================= values =================|======== other options =========\
| 1. Type of perturbative computation | order = NLO | LO |
| 2. No MC@[N]LO matching / event generation | fixed_order = ON ⇐ ̶O̶F̶F̶ ̶ | ON |
| 3. Shower the generated events | shower = OFF ⇐ ̶H̶E̶R̶W̶I̶G̶6̶ ̶ | OFF|PYTHIA6Q|PYTHIA6PT |
| 4. Decay onshell particles | madspin = OFF | ON|onshell |
| 5. Add weights to events for new hypp. | reweight = OFF | ON|NLO|NLO_TREE|LO |
| 6. Run MadAnalysis5 on the events generated | madanalysis = Not Avail. | Please install module |
\==============================================================================================================================/
Either type the switch number (1 to 6) to change its setting,

which looks weird…
can anyone who knows better the business of these switches have a look?

thanks!
cheers,

Marco

> On 4 Oct 2018, at 15:07, Rikkert Frederix <email address hidden> wrote:
>
> It's not enabled. We should. Could you please do this?
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[lp:~maddevelopers/mg5amcnlo/3.0.1]

943. By Rikkert Frederix on 2018-10-04

fixed the switches

[Rikkert Frederix (frederix)]

Hi,

A small change means that the switches are okay now for the QCD stuff. However, for QED, the formatting looks off. I'm not sure we should be too worried about this, since it works correctly. If there is a simple fix, we might apply it.

Cheers,
Rikkert

[lp:~maddevelopers/mg5amcnlo/3.0.1]

944. By marco zaro on 2018-10-04

small fix

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rikkert,

The display is actually perfect for QED (see below):

The following switches determine which programs are run:
/================== Description ==================|=========== values ===========|================ other options ================\
| 1. Type of perturbative computation | order = NLO | LO |
| 2. No MC@[N]LO matching / event generation | fixed_order = ON | No NLO+PS available for EW correction |
\================================================================================================================================/
Either type the switch number (1 to 6) to change its setting,
Set any switch explicitly (e.g. type 'madspin=onshell' at the prompt)
Type 'help' for the list of all valid option
Type '0', 'auto', 'done' or just press enter when you are done.[60s to answer]

Maybe are you referring to the display in debug mode, which include the hidden line (madspin/...)
which are important to see (but only in debug mode). [Also for the hidden option of MadSpin]
In that case, formatting is indeed weird.
I have (slightly) improve it
1) passing the color to green (more consistent with debug coloring in general)
2) fixing the change of length due to the presence of the coloring
This is still not perfect since the length of each column does not take into account those hidden line.
and therefore they can go to overflow.

Cheers,

Olivier

> On 4 Oct 2018, at 17:16, Rikkert Frederix <email address hidden> wrote:
>
> Hi,
>
> A small change means that the switches are okay now for the QCD stuff. However, for QED, the formatting looks off. I'm not sure we should be too worried about this, since it works correctly. If there is a simple fix, we might apply it.
>
> Cheers,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[lp:~maddevelopers/mg5amcnlo/3.0.1]

945. By olivier-mattelaer on 2018-10-04

fix some question formatting for hidden line [debgug mode only]

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

I run (some of) the test this morning (actually 630 of thoses)
and this one sounds quite bad:

ERROR: test_generate_fks_ew (tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW)
check that the generate command works as expected.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/unit_tests/fks/test_extra_ew.py", line 112, in test_generate_fks_ew
    self.assertEqual(set(split), self.interface._fks_multi_proc['splitting_types'])
  File "/Users/omattelaer/Documents/workspace/3.0.1/madgraph/core/base_objects.py", line 69, in __getitem__
    self.is_valid_prop(name) #raise the correct error
  File "/Users/omattelaer/Documents/workspace/3.0.1/madgraph/core/base_objects.py", line 86, in is_valid_prop
    Valid property are %s""" % (name,self.__class__.__name__, self.keys())
PhysicsObjectError: splitting_types is not a valid property for this object: FKSMultiProcess

    Valid property are ['has_isr', 'use_numerical', 'amplitudes', 'born_processes', 'pdgs', 'ignore_six_quark_processes', 'collect_mirror_procs', 'process_definitions', 'OLP', 'real_amplitudes', 'init_lep_split', 'diagram_filter', 'ncores_for_proc_gen', 'has_fsr', 'loop_filter']

some other crashes like this one:

======================================================================
FAIL: testIO_test_wprod_fksew (tests.unit_tests.iolibs.test_export_fks.IOExportFKSTest)
target: SubProcesses/[P0.*\/.+\.(inc|f)]
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 257, in __wrapper
    testKeys=[(testGroup, newTestName)])
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 700, in runIOTests
    self.assertFileContains(open(file_path), goal)
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 389, in assertFileContains
    self.assertEqual(a,b)
AssertionError: ' CALL FFV1_2(W(1,4),W(1,2),GC_3,DCMPLX(ZERO),W(1,8))' != ' CALL FFV1_2(W(1,4),W(1,2),-GC_4,DCMPLX(ZERO),W(1,8))'

Which should be related to model optimization include in 2.6.3. I believe that this is fine (amy does not agree?)

Cheers,

Olivier

On 4 Oct 2018, at 23:24, Olivier Mattelaer <<email address hidden><mailto:<email address hidden>>> wrote:

Hi Rikkert,

The display is actually perfect for QED (see below):

The following switches determine which programs are run:
/================== Description ==================|=========== values ===========|================ other options ================\
| 1. Type of perturbative computation | order = NLO | LO |
| 2. No MC@[N]LO matching / event generation | fixed_order = ON | No NLO+PS available for EW correction |
\================================================================================================================================/
Either type the switch number (1 to 6) to change its setting,
Set any switch explicitly (e.g. type 'madspin=onshell' at the prompt)
T...

[lp:~maddevelopers/mg5amcnlo/3.0.1]

946. By marco zaro on 2018-10-05

fixed one test

[marco zaro (marco-zaro)]

Hi Olivier,
the test you mentioned below was simply not updated after a change I did.
Now it is passing.

Cheers,

Marco

> On 5 Oct 2018, at 09:48, Olivier Mattelaer <email address hidden> wrote:
>
> Hi,
>
> I run (some of) the test this morning (actually 630 of thoses)
> and this one sounds quite bad:
>
> ERROR: test_generate_fks_ew (tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW)
> check that the generate command works as expected.
> ----------------------------------------------------------------------
> Traceback (most recent call last):
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/unit_tests/fks/test_extra_ew.py", line 112, in test_generate_fks_ew
> self.assertEqual(set(split), self.interface._fks_multi_proc['splitting_types'])
> File "/Users/omattelaer/Documents/workspace/3.0.1/madgraph/core/base_objects.py", line 69, in __getitem__
> self.is_valid_prop(name) #raise the correct error
> File "/Users/omattelaer/Documents/workspace/3.0.1/madgraph/core/base_objects.py", line 86, in is_valid_prop
> Valid property are %s""" % (name,self.__class__.__name__, self.keys())
> PhysicsObjectError: splitting_types is not a valid property for this object: FKSMultiProcess
>
> Valid property are ['has_isr', 'use_numerical', 'amplitudes', 'born_processes', 'pdgs', 'ignore_six_quark_processes', 'collect_mirror_procs', 'process_definitions', 'OLP', 'real_amplitudes', 'init_lep_split', 'diagram_filter', 'ncores_for_proc_gen', 'has_fsr', 'loop_filter']
>
> some other crashes like this one:
>
> ======================================================================
> FAIL: testIO_test_wprod_fksew (tests.unit_tests.iolibs.test_export_fks.IOExportFKSTest)
> target: SubProcesses/[P0.*\/.+\.(inc|f)]
> ----------------------------------------------------------------------
> Traceback (most recent call last):
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 257, in __wrapper
> testKeys=[(testGroup, newTestName)])
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 700, in runIOTests
> self.assertFileContains(open(file_path), goal)
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 389, in assertFileContains
> self.assertEqual(a,b)
> AssertionError: ' CALL FFV1_2(W(1,4),W(1,2),GC_3,DCMPLX(ZERO),W(1,8))' != ' CALL FFV1_2(W(1,4),W(1,2),-GC_4,DCMPLX(ZERO),W(1,8))'
>
>
> Which should be related to model optimization include in 2.6.3. I believe that this is fine (amy does not agree?)
>
> Cheers,
>
> Olivier
>
> On 4 Oct 2018, at 23:24, Olivier Mattelaer <<email address hidden> <mailto:<email address hidden>><mailto:<email address hidden> <mailto:<email address hidden>>>> wrote:
>
> Hi Rikkert,
>
> The display is actually perfect for QED (see below):
>
>
> The following switches determine which programs are run:
> /================== Description ==================|=========== values ===========|================ other options ================\
> | 1. Type of perturbative computation | order = NLO | LO |
> | 2. No MC@[N]LO matching / e...

[Rikkert Frederix (frederix)]

Hi,

There is another issue that should be fixed. In particular, when requiring e.g.,

p p > j j QED=0 QCD=4 [QED]

the code also includes QCD corrections (i.e., NLO_1), since at the squared order level, they are allowed. To remove them, one needs to set QCD=2.

In my opinion, the QED=0 and QCD=4 should be used to define which LO contributions should be included, and the ones between the square brackets to determine if one should include QCD and/or QED corrections on top of the Born contributions. Hence,

p p > j j QED=0 QCD=4 [QED]
and
p p > j j QED=0 QCD=2 [QED]

should lead to exactly the same results.

Who knows how to fix this?

Cheers,
Rikkert

[Olivier Mattelaer (olivier-mattelaer)]

On that topic,

Would not be a good idea to put somewhere the plot corresponding to NLO_X LO_Y?
This could be very usefull to know what is computed and what is included by the code.
(what about for 3.0.2?)

Cheers,

Olivier

> On 5 Oct 2018, at 10:48, Rikkert Frederix <email address hidden> wrote:
>
> Hi,
>
> There is another issue that should be fixed. In particular, when requiring e.g.,
>
> p p > j j QED=0 QCD=4 [QED]
>
> the code also includes QCD corrections (i.e., NLO_1), since at the squared order level, they are allowed. To remove them, one needs to set QCD=2.
>
> In my opinion, the QED=0 and QCD=4 should be used to define which LO contributions should be included, and the ones between the square brackets to determine if one should include QCD and/or QED corrections on top of the Born contributions. Hence,
>
> p p > j j QED=0 QCD=4 [QED]
> and
> p p > j j QED=0 QCD=2 [QED]
>
> should lead to exactly the same results.
>
> Who knows how to fix this?
>
> Cheers,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[marco zaro (marco-zaro)]

ciao Rikkert,
I would prefer not to touch the order business now, as I fear to screw up other things with a quick patch.
I agree that it is not optimal, but we may want to upgrade it for the next release with more calm, possibly allowing the user to specify orders at the diagram or amplitude level in a more flexible way.
What is everybody’s idea here?

Cheers,

Marco

> On 5 Oct 2018, at 10:48, Rikkert Frederix <email address hidden> wrote:
>
> Hi,
>
> There is another issue that should be fixed. In particular, when requiring e.g.,
>
> p p > j j QED=0 QCD=4 [QED]
>
> the code also includes QCD corrections (i.e., NLO_1), since at the squared order level, they are allowed. To remove them, one needs to set QCD=2.
>
> In my opinion, the QED=0 and QCD=4 should be used to define which LO contributions should be included, and the ones between the square brackets to determine if one should include QCD and/or QED corrections on top of the Born contributions. Hence,
>
> p p > j j QED=0 QCD=4 [QED]
> and
> p p > j j QED=0 QCD=2 [QED]
>
> should lead to exactly the same results.
>
> Who knows how to fix this?
>
> Cheers,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[Rikkert Frederix (frederix)]

Hi,

Marco: I think you are right. Maybe it's better to wait with this for the next version.

Olivier: I agree with you that for 3.0.2 we should write a better default analysis, which knows about the separate orders, and has a few more plots than just the total cross section.

Cheers,
Rikkert

[Olivier Mattelaer (olivier-mattelaer)]

I was actually not thinking at the analysis level but more at the time of the generation of the diagram
to have the famous bullet plot such that it is simple to see which order are include.

Olivier

> On 5 Oct 2018, at 11:22, Rikkert Frederix <email address hidden> wrote:
>
> Hi,
>
> Marco: I think you are right. Maybe it's better to wait with this for the next version.
>
> Olivier: I agree with you that for 3.0.2 we should write a better default analysis, which knows about the separate orders, and has a few more plots than just the total cross section.
>
> Cheers,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[davide.pagani.85 (davide-pagani)]

Hi,

What Rik suggested is important, the example he is quoting is actually
an error that I did and and made loose quite some time to get what was
going on.

On the other hand, I think is more important to release as soon as
possible this new version and this error is actually due to a wrong
usage of the syntax.

So if it takes long, I would do it for next-to-next release.

-------------------------------------------------------

I think it would be very helpful to add the HwU analysis that Rik wrote
for the results in the paper.

Many times I start from that one or I use it to get fast results.

Do you agree?

Ciao

Davide

On 05/10/18 11:00, marco zaro wrote:
> ciao Rikkert,
> I would prefer not to touch the order business now, as I fear to screw up other things with a quick patch.
> I agree that it is not optimal, but we may want to upgrade it for the next release with more calm, possibly allowing the user to specify orders at the diagram or amplitude level in a more flexible way.
> What is everybody’s idea here?
>
> Cheers,
>
> Marco
>
>
>> On 5 Oct 2018, at 10:48, Rikkert Frederix <email address hidden> wrote:
>>
>> Hi,
>>
>> There is another issue that should be fixed. In particular, when requiring e.g.,
>>
>> p p > j j QED=0 QCD=4 [QED]
>>
>> the code also includes QCD corrections (i.e., NLO_1), since at the squared order level, they are allowed. To remove them, one needs to set QCD=2.
>>
>> In my opinion, the QED=0 and QCD=4 should be used to define which LO contributions should be included, and the ones between the square brackets to determine if one should include QCD and/or QED corrections on top of the Born contributions. Hence,
>>
>> p p > j j QED=0 QCD=4 [QED]
>> and
>> p p > j j QED=0 QCD=2 [QED]
>>
>> should lead to exactly the same results.
>>
>> Who knows how to fix this?
>>
>> Cheers,
>> Rikkert
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>

[Rikkert Frederix (frederix)]

Hi guys,

Any updates?

Best,
Rikkert

[lp:~maddevelopers/mg5amcnlo/3.0.1]

947. By Pagani <email address hidden> on 2018-10-26

General analysis with options of tagging spefic orders explained

[davide.pagani.85 (davide-pagani)]

Hi guys,

sorry for the naive question.
What does it mean merging 3.0.1 into FKS_EW_granny?

Isn't the latter an old version?

Ciao
Davide

On 23/10/18 21:51, Rikkert Frederix wrote:
> Hi guys,
>
> Any updates?
>
> Best,
> Rikkert
>

[Rikkert Frederix (frederix)]

FKS_EW_granny is the current public version ---there is no 3.0.0 branch. Hence, this request is to replace the current public version with the 3.0.1 branch.

Best,
Rikkert

[davide.pagani.85 (davide-pagani)]

Hi Marco,

I see

942. By marco zaro on 2018-10-04
fixes in amcatnlo_run_interface.py to enable NLOPS for QCD-only
not fully working yet

in which sense is not fully working?

Concerning the Fixed-Order part, the code is ok for me to be replaced.

Cheers
Davide

On 26/10/18 19:57, Rikkert Frederix wrote:
> FKS_EW_granny is the current public version ---there is no 3.0.0 branch. Hence, this request is to replace the current public version with the 3.0.1 branch.
>
> Best,
> Rikkert

[lp:~maddevelopers/mg5amcnlo/3.0.1]

948. By olivier-mattelaer on 2018-10-28

merge with 2.6.4

949. By olivier-mattelaer on 2018-10-28

merge with 2.6.4

[marco zaro (marco-zaro)]

Ciao Davide,
this was fixed in rev 943 by rikkert
cheers,

Marco

> On 28 Oct 2018, at 09:15, davide.pagani.85 <email address hidden> wrote:
>
> Hi Marco,
>
> I see
>
> 942. By marco zaro on 2018-10-04
> fixes in amcatnlo_run_interface.py to enable NLOPS for QCD-only
> not fully working yet
>
> in which sense is not fully working?
>
> Concerning the Fixed-Order part, the code is ok for me to be replaced.
>
>
> Cheers
> Davide
>
> On 26/10/18 19:57, Rikkert Frederix wrote:
>> FKS_EW_granny is the current public version ---there is no 3.0.0 branch. Hence, this request is to replace the current public version with the 3.0.1 branch.
>>
>> Best,
>> Rikkert
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[davide.pagani.85 (davide-pagani)]

Then I am ok for the merging.

Let me know how I can help you.

Ciao
Davide

On 29/10/18 09:21, marco zaro wrote:
> Ciao Davide,
> this was fixed in rev 943 by rikkert
> cheers,
>
> Marco
>
>> On 28 Oct 2018, at 09:15, davide.pagani.85 <email address hidden> wrote:
>>
>> Hi Marco,
>>
>> I see
>>
>> 942. By marco zaro on 2018-10-04
>> fixes in amcatnlo_run_interface.py to enable NLOPS for QCD-only
>> not fully working yet
>>
>> in which sense is not fully working?
>>
>> Concerning the Fixed-Order part, the code is ok for me to be replaced.
>>
>>
>> Cheers
>> Davide
>>
>> On 26/10/18 19:57, Rikkert Frederix wrote:
>>> FKS_EW_granny is the current public version ---there is no 3.0.0 branch. Hence, this request is to replace the current public version with the 3.0.1 branch.
>>>
>>> Best,
>>> Rikkert
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>

[lp:~maddevelopers/mg5amcnlo/3.0.1]

950. By olivier-mattelaer on 2018-11-09

UPdate IOtest

951. By olivier-mattelaer on 2018-11-09

last merge with 2.6.4

[Olivier Mattelaer (olivier-mattelaer)]

Hi Marco,

I realised that you removed that file:
tests/acceptance_tests/test_export_fks.py

Is this expected?

Olivier

[Olivier Mattelaer (olivier-mattelaer)]

Also it seems to have a bunch of (scary) print statement in
test_fks_base.py at line 246
Related to some bypass of tests...

Would be better to clean that right?

Olivier

[marco zaro (marco-zaro)]

Hi Olviier,
all IOtests for the NLO output are in
tests/unit_tests/iolibs/test_export_fks.py
feel free to move it to the acceptance tests if it takes too much

cheers,

Marco

> On 9 Nov 2018, at 11:31, Olivier Mattelaer <email address hidden> wrote:
>
> Hi Marco,
>
> I realised that you removed that file:
> tests/acceptance_tests/test_export_fks.py
>
> Is this expected?
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[Olivier Mattelaer (olivier-mattelaer)]

OK good and concerning the test_fks_base.py?
(print statement that some test are bypassed?)

Cheers,

Olivier

PS: I have found workaround for the lhapdf issue and Mojave issue that I pushed in 2.6.4 (and include those here obviously)

> On 9 Nov 2018, at 12:15, marco zaro <email address hidden> wrote:
>
> Hi Olviier,
> all IOtests for the NLO output are in
> tests/unit_tests/iolibs/test_export_fks.py
> feel free to move it to the acceptance tests if it takes too much
>
> cheers,
>
> Marco
>
>> On 9 Nov 2018, at 11:31, Olivier Mattelaer <email address hidden> wrote:
>>
>> Hi Marco,
>>
>> I realised that you removed that file:
>> tests/acceptance_tests/test_export_fks.py
>>
>> Is this expected?
>>
>> Olivier
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[lp:~maddevelopers/mg5amcnlo/3.0.1]

952. By marco zaro on 2018-11-09

commented sprint statament in test_fks_base

[marco zaro (marco-zaro)]

these are commented now (r952)
which LHAPDF issue?
cheers,

Marco

> On 9 Nov 2018, at 12:30, Olivier Mattelaer <email address hidden> wrote:
>
> OK good and concerning the test_fks_base.py?
> (print statement that some test are bypassed?)
>
>
> Cheers,
>
> Olivier
>
>
> PS: I have found workaround for the lhapdf issue and Mojave issue that I pushed in 2.6.4 (and include those here obviously)
>
>> On 9 Nov 2018, at 12:15, marco zaro <email address hidden> wrote:
>>
>> Hi Olviier,
>> all IOtests for the NLO output are in
>> tests/unit_tests/iolibs/test_export_fks.py
>> feel free to move it to the acceptance tests if it takes too much
>>
>> cheers,
>>
>> Marco
>>
>>> On 9 Nov 2018, at 11:31, Olivier Mattelaer <email address hidden> wrote:
>>>
>>> Hi Marco,
>>>
>>> I realised that you removed that file:
>>> tests/acceptance_tests/test_export_fks.py
>>>
>>> Is this expected?
>>>
>>> Olivier
>>> --
>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>>
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[Olivier Mattelaer (olivier-mattelaer)]

They were two lhapdf issues:

- it was not compiling on Mojave (the same problem will likely re-appear in a year with NO easy solution this time). The problem is that the yaml parser shipped with lhapdf seems to have issue on mac.

- the download of (any) PDF was crashing. This was due to the change of hepforge.
We were using the ofiicial command (lhapdf install) but that one does not work anymore and after discussing with them they do not seem interested to fix it.
I'm still trying to use that command, but if it fails then I guess the path of the set and do a "simple" wget on that one.

Cheers,

Olivier

PS: The status on the tests are now:
- one unnitest failing
- around 10 acceptance test failing (so far, this is still running)

> On 9 Nov 2018, at 12:42, marco zaro <email address hidden> wrote:
>
> these are commented now (r952)
> which LHAPDF issue?
> cheers,
>
> Marco
>
>> On 9 Nov 2018, at 12:30, Olivier Mattelaer <email address hidden> wrote:
>>
>> OK good and concerning the test_fks_base.py?
>> (print statement that some test are bypassed?)
>>
>>
>> Cheers,
>>
>> Olivier
>>
>>
>> PS: I have found workaround for the lhapdf issue and Mojave issue that I pushed in 2.6.4 (and include those here obviously)
>>
>>> On 9 Nov 2018, at 12:15, marco zaro <email address hidden> wrote:
>>>
>>> Hi Olviier,
>>> all IOtests for the NLO output are in
>>> tests/unit_tests/iolibs/test_export_fks.py
>>> feel free to move it to the acceptance tests if it takes too much
>>>
>>> cheers,
>>>
>>> Marco
>>>
>>>> On 9 Nov 2018, at 11:31, Olivier Mattelaer <email address hidden> wrote:
>>>>
>>>> Hi Marco,
>>>>
>>>> I realised that you removed that file:
>>>> tests/acceptance_tests/test_export_fks.py
>>>>
>>>> Is this expected?
>>>>
>>>> Olivier
>>>> --
>>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>>>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>>>
>>>
>>> --
>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>>
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[Olivier Mattelaer (olivier-mattelaer)]

Hi Marco,

The following test is also not passing (going trough ok with 2.6.4):
test_madspin_LOonly

Even though this is a madspin designed tests, it actually crash within the LOonly module.
If you just do
geneate p p > w+ [LOonly]
It will lead to a crash.

I have fixed all the other tests (but a couple of IOtests and i still have to run the short paralel test).

Since this is a beta version, we can also decide to not fix this in this version.

[marco zaro (marco-zaro)]

Hi Olivier,
it is strange, it seems that the LOonly mode is broken (it always returns zero)
I will have a look on monday, and let you know.

cheers,

Marco

> On 9 Nov 2018, at 21:10, Olivier Mattelaer <email address hidden> wrote:
>
> Hi Marco,
>
> The following test is also not passing (going trough ok with 2.6.4):
> test_madspin_LOonly
>
> Even though this is a madspin designed tests, it actually crash within the LOonly module.
> If you just do
> geneate p p > w+ [LOonly]
> It will lead to a crash.
>
> I have fixed all the other tests (but a couple of IOtests and i still have to run the short paralel test).
>
> Since this is a beta version, we can also decide to not fix this in this version.
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[marco zaro (marco-zaro)]

Hi Again,
actually, doing p p > e+ ve [LOonly] works ok, at LO, NLO, aMC@LO, aMC@NLO
So it may be something relate to the phase-space of 2->1?
the reason why p p > w+ [LOonly] returns zero is because the function passcuts always returns False, exiting here
c Make sure have reasonable 4-momenta
      if (p(0,1) .le. 0d0) then
         passcuts=.false.
         return
      endif

with P(0,1) = -99

So clearly something goes wrong in the phase space generation. Rikkert, anything that could be related to the granny stuff?

also note that
p p > w+ [QED] has one test failing (soft test for photon emission off the W), where a lot of NaN’s appear
p p > e+ ve [QCD] seems ok though

Cheers,

Marco

> On 9 Nov 2018, at 21:42, marco zaro <email address hidden> wrote:
>
> Hi Olivier,
> it is strange, it seems that the LOonly mode is broken (it always returns zero)
> I will have a look on monday, and let you know.
>
> cheers,
>
> Marco
>
>> On 9 Nov 2018, at 21:10, Olivier Mattelaer <email address hidden> wrote:
>>
>> Hi Marco,
>>
>> The following test is also not passing (going trough ok with 2.6.4):
>> test_madspin_LOonly
>>
>> Even though this is a madspin designed tests, it actually crash within the LOonly module.
>> If you just do
>> geneate p p > w+ [LOonly]
>> It will lead to a crash.
>>
>> I have fixed all the other tests (but a couple of IOtests and i still have to run the short paralel test).
>>
>> Since this is a beta version, we can also decide to not fix this in this version.
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[Olivier Mattelaer (olivier-mattelaer)]

Thanks,

Update on my side:
1) the IOTest suite are bugged (not critical for a release) I have ask Valentin to look at it
2) the (fast) parrelel test leads to their own discovery of bug:

2.1) e+ e- > t t~ [real=QCD]
crashes with
FKSProcessError: This process does not have any correction up to NLO in QED,QCD
[ test: test_short_check_generate_events_nlo_py6pt_fsr]
same for e+ e- > p p p [real=QCD] in test: test_short_check_eejjj_lo_lhapdf]

2.2) A division by zero occurs in the jet veto of rik [test_short_jet_veto_xsec]

p p > e+ ve [QCD]
generate_events NLO
                 set ickkw -1
                 set ptj 10

2.3) The value of some squared order amplitude does not match anymore [test_short_sqso ]
Will look at this one (maybe with Valentin).

Cheers,

Olivier

On 9 Nov 2018, at 22:00, marco zaro <<email address hidden><mailto:<email address hidden>>> wrote:

Hi Again,
actually, doing p p > e+ ve [LOonly] works ok, at LO, NLO, aMC@LO, aMC@NLO
So it may be something relate to the phase-space of 2->1?
the reason why p p > w+ [LOonly] returns zero is because the function passcuts always returns False, exiting here
c Make sure have reasonable 4-momenta
     if (p(0,1) .le. 0d0) then
        passcuts=.false.
        return
     endif

with P(0,1) = -99

So clearly something goes wrong in the phase space generation. Rikkert, anything that could be related to the granny stuff?

also note that
p p > w+ [QED] has one test failing (soft test for photon emission off the W), where a lot of NaN’s appear
p p > e+ ve [QCD] seems ok though

Cheers,

Marco

On 9 Nov 2018, at 21:42, marco zaro <<email address hidden><mailto:<email address hidden>>> wrote:

Hi Olivier,
it is strange, it seems that the LOonly mode is broken (it always returns zero)
I will have a look on monday, and let you know.

cheers,

Marco

On 9 Nov 2018, at 21:10, Olivier Mattelaer <<email address hidden><mailto:<email address hidden>>> wrote:

Hi Marco,

The following test is also not passing (going trough ok with 2.6.4):
test_madspin_LOonly

Even though this is a madspin designed tests, it actually crash within the LOonly module.
If you just do
geneate p p > w+ [LOonly]
It will lead to a crash.

I have fixed all the other tests (but a couple of IOtests and i still have to run the short paralel test).

Since this is a beta version, we can also decide to not fix this in this version.
--
https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

--
https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

--
https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[marco zaro (marco-zaro)]

Ciao Olivier,
thanks for looking at this.
Let me answer inline

On Fri, Nov 9, 2018 at 10:30 PM Olivier Mattelaer <
<email address hidden>> wrote:

> Thanks,
>
> Update on my side:
> 1) the IOTest suite are bugged (not critical for a release) I have ask
> Valentin to look at it
> 2) the (fast) parrelel test leads to their own discovery of bug:
>
> 2.1) e+ e- > t t~ [real=QCD]
> crashes with
> FKSProcessError: This process does not have any correction up to NLO in
> QED,QCD
> [ test: test_short_check_generate_events_nlo_py6pt_fsr]
> same for e+ e- > p p p [real=QCD] in test:
> test_short_check_eejjj_lo_lhapdf]
>
> for ee collisions, one has to tell the code
set include_lepton_initiated_processes True (there is a message printed on
screen).
The py6isr test should be skipped, as the has_isr flag is not updated
anymore in this branch (there were some print statement about this you
asked me to remove)

  If you want to include it, set the

  'include_lepton_initiated_processes' option to True

> 2.2) A division by zero occurs in the jet veto of rik
> [test_short_jet_veto_xsec]
>
> p p > e+ ve [QCD]
> generate_events NLO
> set ickkw -1
> set ptj 10
>
>
> the jet_veto is likely to be broken in this branch. There are some message
errors in the logs of the fortran executables about this.

> 2.3) The value of some squared order amplitude does not match anymore
> [test_short_sqso ]
> Will look at this one (maybe with Valentin).
>

I will check this right away

Thanks again!

cheers,

Marco

>
> Cheers,
>
> Olivier
>
> On 9 Nov 2018, at 22:00, marco zaro <<email address hidden><mailto:
> <email address hidden>>> wrote:
>
> Hi Again,
> actually, doing p p > e+ ve [LOonly] works ok, at LO, NLO, aMC@LO, aMC@NLO
> So it may be something relate to the phase-space of 2->1?
> the reason why p p > w+ [LOonly] returns zero is because the function
> passcuts always returns False, exiting here
> c Make sure have reasonable 4-momenta
> if (p(0,1) .le. 0d0) then
> passcuts=.false.
> return
> endif
>
> with P(0,1) = -99
>
> So clearly something goes wrong in the phase space generation. Rikkert,
> anything that could be related to the granny stuff?
>
> also note that
> p p > w+ [QED] has one test failing (soft test for photon emission off the
> W), where a lot of NaN’s appear
> p p > e+ ve [QCD] seems ok though
>
> Cheers,
>
> Marco
>
>
>
> On 9 Nov 2018, at 21:42, marco zaro <<email address hidden><mailto:
> <email address hidden>>> wrote:
>
> Hi Olivier,
> it is strange, it seems that the LOonly mode is broken (it always returns
> zero)
> I will have a look on monday, and let you know.
>
> cheers,
>
> Marco
>
> On 9 Nov 2018, at 21:10, Olivier Mattelaer <<email address hidden>
> <mailto:<email address hidden>>> wrote:
>
> Hi Marco,
>
> The following test is also not passing (going trough ok with 2.6.4):
> test_madspin_LOonly
>
> Even though this is a madspin designed tests, it actually crash within the
> LOonly module.
> If you just do
> geneate p p > w+ [LOonly]
> It will lead to a crash.
>
> I have fixed all the other tests (but a couple of IOtests and i still have
> to run the short ...

[Rikkert Frederix (frederix)]

Hi Olivier, Marco,

I looked into the following two issues:

1. The error in the phase-space for the process 'p p > w+ [LOonly=QCD]'.
This is related to the fact that in the 3.0.1 branch the final state W boson is the fks_j particle, while in the 2.6.4 branch it's one of the initial state quarks. However, the phase-space for 2->1 processes with j-fks the massive final state particle does not work. There are divisions by zero everywhere... I'm not sure if there is an easy fix here, apart from telling the code to use an initial state particle as fake j_fks. This will work for the LOonly 2->1 processes. On the other hand, we need to fix this phase-space at some point, since the same is used (and needed!) for the process 'p p > w+ [QED]'.

2. The error for the VETO X-SEC (ickkw=-1).
The whole split-amplitude business has not been implemented for the veto-xsec stuff. I suggest that we do not try to rush this through now, but leave the 3.0.1 code as beta and fix this at a later stage.

Cheers,
Rik

[lp:~maddevelopers/mg5amcnlo/3.0.1]

953. By marco zaro on 2018-11-12

fix for 2->1 processes

954. By marco zaro on 2018-11-12

unit test for previous push and fix for p p > w+ [LOonly=QCD]

955. By marco zaro on 2018-11-16

reverted change from revision 613 (related to split orders and WEIGHTED)
added a debug message

[marco zaro (marco-zaro)]

Hi all,
I have just noticed this (fks_singular, line 3818)
      call Qterms_reduced_timelike(j_type,i_type,ch_i,ch_j,t,z,Q)
it should be instead
      call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,t,z,Q)
(switch ch_i <--> ch_j), shouldn't it?

It may be that in the end there is no effect, since, for the relevant cases when either ch_i or ch_j ==0 and the other is not one has always Q=0

c g/a ->qqbar splitting
         Qterms(1) = 4d0 * TR * z*(1d0-z)**2
         Qterms(2) = 4d0 * dble(abs(col1)) * ch1**2 * z*(1d0-z)**2

      elseif ((abs(col1).eq.3 .and. col2.eq.8) .or.
     & (dabs(ch1).gt.0d0 .and. dabs(ch2).eq.0d0)) then
c q->q g/a splitting
         Qterms(1) = 0d0
         Qterms(2) = 0d0

Let me know

Cheers,

Marco

review: Needs Fixing

[marco zaro (marco-zaro)]

sorry, the piece of code i wanted to post is

      elseif ((abs(col1).eq.3 .and. col2.eq.8) .or.
     & (dabs(ch1).gt.0d0 .and. dabs(ch2).eq.0d0)) then
c q->q g/a splitting
         Qterms(1) = 0d0
         Qterms(2) = 0d0

      elseif ((col1.eq.8 .and. abs(col2).eq.3) .or.
     & (dabs(ch1).eq.0d0 .and. dabs(ch2).gt.0d0)) then
c q->g/a q splitting
         Qterms(1) = 0d0
         Qterms(2) = 0d0
      else

> On 21 Nov 2018, at 16:16, marco zaro <email address hidden> wrote:
>
> Review: Needs Fixing
>
> Hi all,
> I have just noticed this (fks_singular, line 3818)
> call Qterms_reduced_timelike(j_type,i_type,ch_i,ch_j,t,z,Q)
> it should be instead
> call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,t,z,Q)
> (switch ch_i <--> ch_j), shouldn't it?
>
> It may be that in the end there is no effect, since, for the relevant cases when either ch_i or ch_j ==0 and the other is not one has always Q=0
>
>
> c g/a ->qqbar splitting
> Qterms(1) = 4d0 * TR * z*(1d0-z)**2
> Qterms(2) = 4d0 * dble(abs(col1)) * ch1**2 * z*(1d0-z)**2
>
> elseif ((abs(col1).eq.3 .and. col2.eq.8) .or.
> & (dabs(ch1).gt.0d0 .and. dabs(ch2).eq.0d0)) then
> c q->q g/a splitting
> Qterms(1) = 0d0
> Qterms(2) = 0d0
>
> Let me know
>
> Cheers,
>
> Marco
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[Rikkert Frederix (frederix)]

Hi Marco,

I agree that the current way is wrong, but has no effect on the running of the code whatsoever. Even though, we should replace ch_i and ch_j in the call to the Qterms.
Note that, if there was 'true' bug here, the soft and/or collinear checks and the beginning of any run would not have passed.

Cheers,
Rik

[lp:~maddevelopers/mg5amcnlo/3.0.1]

956. By marco zaro on 2018-11-21

fix in a call to
Qterms_reduced_timelike
luckiliy harmless

[marco zaro (marco-zaro)]

Hi Rik,
thanks for the cross-check
I have pushed the fix into 3.0.1, rev 956
cheers,

Marco

> On 21 Nov 2018, at 16:24, Rikkert Frederix <email address hidden> wrote:
>
> Hi Marco,
>
> I agree that the current way is wrong, but has no effect on the running of the code whatsoever. Even though, we should replace ch_i and ch_j in the call to the Qterms.
> Note that, if there was 'true' bug here, the soft and/or collinear checks and the beginning of any run would not have passed.
>
> Cheers,
> Rik
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

I would like to put this quite high priority for this week (i.e. for me Wednesday-friday).

I have rerun the tests and found that some tests are crashing when run via the full script but are working if run alone. So I have run (the very slow script to determine which of the tests had such side effects and the culprint is:
tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW.test_generate_fks_2to1_no_finalstate_confs
which has a side effect that lead to the following crash:

FAIL: test_generate_fks_ew_dijet_qed2qcd2_qcd (tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW)
check that the generate command works as expected.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/unit_tests/fks/test_extra_ew.py", line 869, in test_generate_fks_ew_dijet_qed2qcd2_qcd
    self.assertEqual(real.fks_infos[0]['extra_cnt_index'], -1)
AssertionError: 0 != -1

Might be interesting that you take a look. But I think that we should release this version this week independently.

Any objection? comment?

Olivier

[lp:~maddevelopers/mg5amcnlo/3.0.1]

957. By marco zaro on 2018-12-17

tests/unit_tests/fks/test_extra_ew.py has been modified to avoid
side effects

958. By marco zaro on 2018-12-17

fix in amcatnlo_interface to prevent unphysical coupling orders

[marco zaro (marco-zaro)]

Hello Olivier,
I agree we should release 3.0.1 before the Christmas break.
concerning the side effect, i reproduced it, and it goes away if I do the following change
 class TestAMCatNLOEW(unittest.TestCase):
     """ a suite of extra tests for the ew stuff """

- interface = mgcmd.MasterCmd()
+ def setUp(self):
+ self.interface = mgcmd.MasterCmd()

I have pushed the fix into rev 957.

The other bogus behavious was related to the order settings in t-channel single top, which we discussed in a separate email.
I have added an error message that should warn the user to specify the coupling orders in that case (rev 958)

Cheers,

Marco

> On 17 Dec 2018, at 09:23, Olivier Mattelaer <email address hidden> wrote:
>
> Hi,
>
> I would like to put this quite high priority for this week (i.e. for me Wednesday-friday).
>
> I have rerun the tests and found that some tests are crashing when run via the full script but are working if run alone. So I have run (the very slow script to determine which of the tests had such side effects and the culprint is:
> tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW.test_generate_fks_2to1_no_finalstate_confs
> which has a side effect that lead to the following crash:
>
> FAIL: test_generate_fks_ew_dijet_qed2qcd2_qcd (tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW)
> check that the generate command works as expected.
> ----------------------------------------------------------------------
> Traceback (most recent call last):
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/unit_tests/fks/test_extra_ew.py", line 869, in test_generate_fks_ew_dijet_qed2qcd2_qcd
> self.assertEqual(real.fks_infos[0]['extra_cnt_index'], -1)
> AssertionError: 0 != -1
>
> Might be interesting that you take a look. But I think that we should release this version this week independently.
>
> Any objection? comment?
>
> Olivier
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[Rikkert Frederix (frederix)]

Hi Olivier,

Thanks! We should really release it before the break. It will remain a beta version and exist along side the main branch, right? I think that would still be better and removing the beta tag.

best,
Rikkert

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

> Thanks! We should really release it before the break. It will remain a beta version and exist along side the main branch, right?

Correct.

> I think that would still be better and removing the beta tag.

It would not make sense to remove the "beta" tag and keep two official version.
And I do not think that we are ready to keep only this version.

Cheers,

Olivier

> On 18 Dec 2018, at 10:33, Rikkert Frederix <email address hidden> wrote:
>
> Hi Olivier,
>
> Thanks! We should really release it before the break. It will remain a beta version and exist along side the main branch, right? I think that would still be better and removing the beta tag.
>
> best,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[Rikkert Frederix (frederix)]

Sorry, I forgot a "NOT" in my suggestion :-D
Hence I agree with not removing the beta tag.

Cheers,
Rik

[lp:~maddevelopers/mg5amcnlo/3.0.1]

959. By olivier-mattelaer on 2018-12-19

fixing some additional issue with some test --in particular a real bug for the reweighting--

[Olivier Mattelaer (olivier-mattelaer)]

I have re-run all the tests today (and fix a couple of small stuffs). Except a couple that we have decide to ignore above they all go trough now. So I will make the release tomorrow.

Cheers,

Olivier

review: Approve

[Rikkert Frederix (frederix)]

Thanks a lot, Olivier.

Cheers,
Rik

[marco zaro (marco-zaro)]

Thanks a lot Olivier!
great that we finally made it!

cheers,

Marco

> On 19 Dec 2018, at 23:42, Rikkert Frederix <email address hidden> wrote:
>
> Thanks a lot, Olivier.
>
> Cheers,
> Rik
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

[lp:~maddevelopers/mg5amcnlo/3.0.1]

960. By olivier-mattelaer on 2018-12-20

UpdateNotes/Version update

[Olivier Mattelaer (olivier-mattelaer)]

Done!

[davide.pagani.85 (davide-pagani)]

Thanks a lot!

Cheers
Davide

On 20/12/18 09:47, Olivier Mattelaer wrote:
> Done!

[Olivier Mattelaer (olivier-mattelaer)]

By the way,

Do you know if we can close the following bug report:
https://bugs.launchpad.net/mg5amcnlo/+bug/1779831

If not we should at least comment on it (this is one was forbidding ATLAS to use the code)

Cheers,

Olivier

[Stefano Frixione (stefano-frixione)]

I have no idea.
Perhaps one should double check with the ATLAS people;
I've never understood how genser works, and their relationship
with experiments.
Cheers, Stefano.

On Fri, 21 Dec 2018, Olivier Mattelaer wrote:

> By the way,
>
> Do you know if we can close the following bug report:
> https://bugs.launchpad.net/mg5amcnlo/+bug/1779831
>
> If not we should at least comment on it (this is one was forbidding
> ATLAS to use the code)
>
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2018-04-25 06:47:06 +0000
+++ MadSpin/decay.py	2018-12-20 08:18:55 +0000
@@ -3963,14 +3963,23 @@
                 elif mode=='full':
                     stdin_text="5 0 0 0 \n"  # before closing, write down the seed 
                     external = self.calculator[('full',path)]
-                    external.stdin.write(stdin_text)
+                    try:
+                        external.stdin.write(stdin_text)
+                    except Exception:
+                        continue
                     ranmar_state=external.stdout.readline()
                     ranmar_file=pjoin(path,'ranmar_state.dat')
                     ranmar=open(ranmar_file, 'w')
                     ranmar.write(ranmar_state)
                     ranmar.close()
-                    external.stdin.close()
-                    external.stdout.close()
+                    try:
+                        external.stdin.close()
+                    except Exception:
+                        continue
+                    try:
+                        external.stdout.close()
+                    except Exception:
+                        continue
                     external.terminate()
                     del external
         else:

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2018-03-11 12:10:19 +0000
+++ MadSpin/interface_madspin.py	2018-12-20 08:18:55 +0000
@@ -35,6 +35,7 @@
 import madgraph.interface.madgraph_interface as mg_interface
 import madgraph.interface.master_interface as master_interface
 import madgraph.interface.madevent_interface as madevent_interface
+import madgraph.interface.common_run_interface as common_run_interface
 import madgraph.interface.reweight_interface as rwgt_interface
 import madgraph.various.misc as misc
 import madgraph.iolibs.files as files
@@ -53,6 +54,68 @@
 cmd_logger = logging.getLogger('cmdprint2') # -> print
 
 
+class MadSpinOptions(banner.ConfigFile):
+    
+    def default_setup(self):
+
+        self.add_param("max_weight", -1)
+        self.add_param('curr_dir', os.path.realpath(os.getcwd()))
+        self.add_param('Nevents_for_max_weigth', 0)
+        self.add_param("max_weight_ps_point", 400)
+        self.add_param('BW_cut', -1)
+        self.add_param('nb_sigma', 0.)
+        self.add_param('ms_dir', '')
+        self.add_param('max_running_process', 100)
+        self.add_param('onlyhelicity', False)
+        self.add_param('spinmode', "madspin", allowed=['madspin','none','onshell'])
+        self.add_param('use_old_dir', False, comment='should be use only for faster debugging')
+        self.add_param('run_card', '' , comment='define cut for spinmode==none. Path to run_card to use')
+        self.add_param('fixed_order', False, comment='to activate fixed order handling of counter-event')
+        self.add_param('seed', 0, comment='control the seed of madspin')
+        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'])
+        
+        
+    ############################################################################
+    ##  Special post-processing of the options                                ## 
+    ############################################################################
+    def post_set_ms_dir(self, value, change_userdefine, raiseerror):
+        """ special handling for set ms_dir """
+        
+        self.__setitem__('curr_dir', value, change_userdefine=change_userdefine)
+        
+    ############################################################################
+    def post_set_seed(self, value, change_userdefine, raiseerror):
+        """ special handling for set seed """
+        
+        random.seed(value)
+
+    ############################################################################        
+    def post_set_run_card(self, value, change_userdefine, raiseerror):
+        """ special handling for set run_card """
+        
+        if value == 'default':
+            self.run_card = None
+        elif os.path.isfile(value):
+            self.run_card = banner.RunCard(value)
+            
+        args = value.split()
+        if  len(args) >2:
+            if not self.options['run_card']:
+                self.run_card =  banner.RunCardLO()
+                self.run_card.remove_all_cut()
+            self.run_card[args[0]] = ' '.join(args[1:])
+            
+        
+    ############################################################################
+    def post_fixed_order(self, value, change_userdefine, raiseerror):
+        """ special handling for set fixed_order """
+        
+        if value:
+            logger.warning('Fix order madspin fails to have the correct scale information. This can bias the results!')
+            logger.warning('Not all functionalities of MadSpin handle this mode correctly (only onshell mode so far).')
+
 
 class MadSpinInterface(extended_cmd.Cmd):
     """Basic interface for madspin"""
@@ -77,22 +140,7 @@
         self.mode = "madspin" # can be flat/bridge change the way the decay is done.
                               # note amc@nlo does not support bridge.
         
-        self.options = {'max_weight': -1, 
-                        'curr_dir': os.path.realpath(os.getcwd()),
-                        'Nevents_for_max_weigth': 0,
-                        'max_weight_ps_point': 400,
-                        'BW_cut':-1,
-                        'nb_sigma':0,
-                        'ms_dir':None,
-                        'max_running_process':100,
-                        'onlyhelicity': False,
-                        'spinmode': "madspin",
-                        'use_old_dir': False, #should be use only for faster debugging
-                        'run_card': None, # define cut for spinmode==none.
-                        'fixed_order': False # to activate fixed order handling of counter-event
-                        }
-        
-
+        self.options = MadSpinOptions()
         
         self.events_file = None
         self.decay_processes = {}
@@ -107,7 +155,21 @@
             logger.info("Extracting the banner ...")
             self.do_import(event_path)
             
-            
+    
+    def setup_for_pure_decay(self):
+        """this is for spinmode=None -> simple decay
+           We go here if they are no banner.
+           -> this requires that a command import model appears in the card!
+        """
+        
+        logger.info("Setup the code for pure decay mode")
+        self.proc_option = []
+        self.final_state_full = ''
+        self.final_state_compact = ''
+        self.prod_branches = ''
+        self.final_state = set()
+        
+     
     def do_import(self, inputfile):
         """import the event file"""
         
@@ -126,6 +188,7 @@
         if not os.path.exists(inputfile):
             if inputfile.endswith('.gz'):
                 if not os.path.exists(inputfile[:-3]):
+                    misc.sprint(os.getcwd(), os.listdir('.'), inputfile, os.path.exists(inputfile), os.path.exists(inputfile[:-3]))
                     raise self.InvalidCmd('No such file or directory : %s' % inputfile)
                 else: 
                     inputfile = inputfile[:-3]
@@ -134,14 +197,22 @@
             else: 
                 raise self.InvalidCmd('No such file or directory : %s' % inputfile)
 
+        self.inputfile = inputfile
+        if self.options['spinmode'] == 'none' and \
+           (self.options['input_format'] not in ['lhe','auto'] or 
+             (self.options['input_format'] == 'auto' and '.lhe'  not in inputfile[-7:])):  
+            self.banner = banner.Banner()
+            self.setup_for_pure_decay()
+            return  
+        
         if inputfile.endswith('.gz'):
             misc.gunzip(inputfile)
             inputfile = inputfile[:-3]
-
         # 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:
             self.events_file = None
@@ -163,7 +234,6 @@
                 self.options['nb_sigma'] = N_sigma
             if self.options['BW_cut'] == -1:
                 self.options['BW_cut'] = float(self.banner.get_detail('run_card', 'bwcutoff'))
-
         else:
             if not self.options['Nevents_for_max_weigth']:
                 self.options['Nevents_for_max_weigth'] = 75
@@ -273,7 +343,7 @@
         
 
         #Check the param_card
-        if not bypass_check:
+        if not (bypass_check or self.options['input_format'] in ['hepmc', 'lhe_no_banner']):
             if not hasattr(self.banner, 'param_card'):
                 self.banner.charge_card('slha')
             param_card = check_param_card.ParamCard(card)
@@ -395,42 +465,10 @@
         
         args = self.split_arg(line)
         self.check_set(args)
-        
-        if args[0] not in ['ms_dir', 'run_card']:
-            args[1] = args[1].lower()
-        
-        if args[0] in  ['max_weight', 'BW_effect','ms_dir', 'spinmode']:
-            self.options[args[0]] = args[1]
-            if args[0] == 'ms_dir':
-                self.options['curr_dir'] = self.options['ms_dir']
-        elif args[0] == 'seed':
-            random.seed(int(args[1]))
-            self.seed = int(args[1])
-        elif args[0] == 'BW_cut':
-            self.options[args[0]] = float(args[1])
-        elif args[0] in ['onlyhelicity', 'use_old_dir']:
-            self.options[args[0]] = banner.ConfigFile.format_variable(args[1], bool, args[0])
-        elif args[0] in ['run_card']:
-            if args[1] == 'default':
-                self.options['run_card'] = None
-            elif os.path.isfile(args[1]):
-                self.options['run_card'] = banner.RunCard(args[1])
-            elif  len(args) >2:
-                if not self.options['run_card']:
-                    self.options['run_card'] =  banner.RunCardLO()
-                    self.options['run_card'].remove_all_cut()
-                self.options['run_card'][args[1]] = args[2]
-        elif args[0] in 'fixed_order':
-            if args[1].lower() in ['t', 'true']:
-                self.options['fixed_order'] = True
-                logger.warning('Fix order madspin fails to have the correct scale information. This can bias the results!')
-                logger.warning('Not all functionality of MadSpin handle this mode correctly (only onshell mode so far).')
-            else:
-                self.options['fixed_order'] = False
-            logger.info('fixed_order options set to %s', self.options['fixed_order'])
-        else:
-            self.options[args[0]] = int(args[1])
-    
+
+        self.options[args[0]] = ' '.join(args[1:])
+        
+
     def complete_set(self,  text, line, begidx, endidx):
         
 
@@ -559,17 +597,17 @@
 
         model_line = self.banner.get('proc_card', 'full_model_line')
 
-        if not self.seed:
-            self.seed = random.randint(0, int(30081*30081))
-            self.do_set('seed %s' % self.seed)
-            logger.info('Will use seed %s' % self.seed)
-            self.history.insert(0, 'set seed %s' % self.seed)
+        if not self.options['seed']:
+            self.options['seed'] = random.randint(0, int(30081*30081))
+            #self.do_set('seed %s' % self.seed)
+            logger.info('Will use seed %s' % self.options['seed'])
+            self.history.insert(0, 'set seed %s' % self.options['seed'])
 
-        if self.seed > 30081*30081: # can't use too big random number
-            msg = 'Random seed too large ' + str(self.seed) + ' > 30081*30081'
+        if self.options['seed'] > 30081*30081: # can't use too big random number
+            msg = 'Random seed too large ' + str(self.options['seed']) + ' > 30081*30081'
             raise Exception, msg
 
-        self.options['seed'] = self.seed
+        #self.options['seed'] = self.seed
         text = '%s\n' % '\n'.join([ line for line in self.history if line])
         self.banner.add_text('madspin' , text)
         
@@ -678,9 +716,9 @@
         
         # NOW we have all the information available for RUNNING
         
-        if self.seed:
+        if self.options['seed']:
             #seed is specified need to use that one:
-            open(pjoin(self.options['ms_dir'],'seeds.dat'),'w').write('%s\n'%self.seed)
+            open(pjoin(self.options['ms_dir'],'seeds.dat'),'w').write('%s\n'%self.options['seed'])
             #remove all ranmar_state
             for name in misc.glob(pjoin('*', 'SubProcesses','*','ranmar_state.dat'), 
                                                         self.options['ms_dir']):
@@ -741,30 +779,61 @@
             for name in misc.glob("decay_*_*", self.path_me):
                 shutil.rmtree(name)
 
-        self.events_file.close()
-        orig_lhe = lhe_parser.EventFile(self.events_file.name)
-        
+        if self.events_file:
+            self.events_file.close()
+            filename = self.events_file.name
+        else:
+            filename = self.inputfile
+
+        if self.options['input_format'] == 'auto':
+            if '.lhe' in filename :
+                self.options['input_format']  = 'lhe'
+            elif '.hepmc' in filename:
+                self.options['input_format']  = 'hepmc'
+            else:
+                raise Exception, "fail to recognized input format automatically"
+                
+        if self.options['input_format'] in ['lhe', 'lhe_no_banner']:
+            orig_lhe = lhe_parser.EventFile(filename)
+            if self.options['input_format'] == 'lhe_no_banner':
+                orig_lhe.allow_empty_event = True
+                
+        elif self.options['input_format'] in ['hepmc']:
+            import madgraph.various.hepmc_parser as hepmc_parser
+            orig_lhe = hepmc_parser.HEPMC_EventFile(filename)
+            logger.info("Parsing input event to know how many decay to generate. This can takes few minuts.")
+        else:
+            raise Exception
+            
         to_decay = collections.defaultdict(int)
         nb_event = 0
+ 
         for event in orig_lhe:
             nb_event +=1
             for particle in event:
                 if particle.status == 1 and particle.pdg in asked_to_decay:
                     # final state and tag as to decay
                     to_decay[particle.pdg] += 1
+            if self.options['input_format'] == 'hepmc' and nb_event == 250:
+                currpos = orig_lhe.tell()
+                filesize = orig_lhe.getfilesize()
+                for key in to_decay:
+                    to_decay[key] *= 1.05 * filesize/ currpos 
+                    # 1.05 to avoid accidental coincidence with nevents
+                break
 
         # Handle the banner of the output file
-        if not self.seed:
-            self.seed = random.randint(0, int(30081*30081))
-            self.do_set('seed %s' % self.seed)
-            logger.info('Will use seed %s' % self.seed)
-            self.history.insert(0, 'set seed %s' % self.seed)
+        if not self.options['seed']:
+            self.options['seed'] = random.randint(0, int(30081*30081))
+            #self.do_set('seed %s' % self.seed)
+            logger.info('Will use seed %s' % self.options['seed'])
+            self.history.insert(0, 'set seed %s' % self.options['seed'])
 
-        if self.seed > 30081*30081: # can't use too big random number
-            msg = 'Random seed too large ' + str(self.seed) + ' > 30081*30081'
+        if self.options['seed'] > 30081*30081: # can't use too big random number
+            msg = 'Random seed too large ' + str(self.options['seed']) + ' > 30081*30081'
             raise Exception, msg
 
-        self.options['seed'] = self.seed
+        #self.options['seed'] = self.options['seed']
         
         text = '%s\n' % '\n'.join([ line for line in self.history if line])
         self.banner.add_text('madspin' , text)
@@ -780,7 +849,7 @@
             for pdg, nb_needed in to_decay.items():
                 #check if a splitting is needed
                 if nb_needed == nb_event:
-                    evt_decayfile[pdg] = self.generate_events(pdg, nb_needed, mg5)
+                    evt_decayfile[pdg] = self.generate_events(pdg, min(nb_needed,100000), mg5)
                 elif nb_needed %  nb_event == 0:
                     nb_mult = nb_needed // nb_event
                     part = self.model.get_particle(pdg)
@@ -788,64 +857,107 @@
                     if name not in self.list_branches:
                         continue
                     elif len(self.list_branches[name]) == nb_mult:
-                        evt_decayfile[pdg] = self.generate_events(pdg, nb_event, mg5)
-                    else:
-                        evt_decayfile[pdg] = self.generate_events(pdg, nb_needed, mg5, cumul=True)
+                        evt_decayfile[pdg] = self.generate_events(pdg, min(nb_event,100000), mg5)
+                    else:
+                        evt_decayfile[pdg] = self.generate_events(pdg, min(nb_needed,100000), mg5, cumul=True)
+                elif self.options['cross_section']:
+                    #cross-section hard-coded -> allow 
+                    part = self.model.get_particle(pdg)
+                    name = part.get_name()
+                    
+                    if name not in self.list_branches:
+                        continue
+                    else:
+                        try:
+                            evt_decayfile[pdg] = self.generate_events(pdg, min(nb_needed,100000), mg5, cumul=True)
+                        except common_run_interface.ZeroResult:
+                            logger.warning("Branching ratio is zero for this particle. Not decaying it")
+                            del to_decay[pdg]                    
                 else:
                     part = self.model.get_particle(pdg)
                     name = part.get_name()
                     if name not in self.list_branches or len(self.list_branches[name]) == 0:
                         continue
-                    raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
+                    raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed. One workaround is to force the final cross-section manually.")
                     
                      
-                
-        
         # Compute the branching ratio.
-        br = 1
-        for (pdg, event_files) in evt_decayfile.items():
-            if not event_files:
-                continue
-            totwidth = float(self.banner.get('param', 'decay', abs(pdg)).value)
-            if to_decay[pdg] == nb_event:
-                # Exactly one particle of this type to decay by event
-                pwidth = sum([event_files[k].cross for k in event_files])
-                if pwidth > 1.01 * totwidth:
-                    logger.critical("Branching ratio larger than one for %s " % pdg) 
-                br *= pwidth / totwidth
-            elif to_decay[pdg] % nb_event == 0:
-                # More than one particle of this type to decay by event
-                # Need to check the number of event file to check if we have to 
-                # make separate type of decay or not.
-                nb_mult = to_decay[pdg] // nb_event
-                if nb_mult == len(event_files):
-                    for k in event_files:
-                        pwidth = event_files[k].cross
-                        if pwidth > 1.01 * totwidth:
-                            logger.critical("Branching ratio larger than one for %s " % pdg)                       
-                        br *= pwidth / totwidth
-                    br *= math.factorial(nb_mult)
-                else:
-                    pwidth = sum(event_files[k].cross for k in event_files)
+        if not self.options['cross_section']:
+            br = 1
+            for (pdg, event_files) in evt_decayfile.items():
+                if not event_files:
+                    continue
+                totwidth = float(self.banner.get('param', 'decay', abs(pdg)).value)
+                if to_decay[pdg] == nb_event:
+                    # Exactly one particle of this type to decay by event
+                    pwidth = sum([event_files[k].cross for k in event_files])
                     if pwidth > 1.01 * totwidth:
                         logger.critical("Branching ratio larger than one for %s " % pdg) 
-                    br *= (pwidth / totwidth)**nb_mult
-            else:
-                raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
+                    br *= pwidth / totwidth
+                elif to_decay[pdg] % nb_event == 0:
+                    # More than one particle of this type to decay by event
+                    # Need to check the number of event file to check if we have to 
+                    # make separate type of decay or not.
+                    nb_mult = to_decay[pdg] // nb_event
+                    if nb_mult == len(event_files):
+                        for k in event_files:
+                            pwidth = event_files[k].cross
+                            if pwidth > 1.01 * totwidth:
+                                logger.critical("Branching ratio larger than one for %s " % pdg)                       
+                            br *= pwidth / totwidth
+                        br *= math.factorial(nb_mult)
+                    else:
+                        pwidth = sum(event_files[k].cross for k in event_files)
+                        if pwidth > 1.01 * totwidth:
+                            logger.critical("Branching ratio larger than one for %s " % pdg) 
+                        br *= (pwidth / totwidth)**nb_mult
+                else:
+                    raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
+        else:
+            br = 1
         self.branching_ratio = br
         self.efficiency = 1
-        self.cross, self.error = self.banner.get_cross(witherror=True)
+        try:
+            self.cross, self.error = self.banner.get_cross(witherror=True)
+        except:
+            if self.options['input_format'] != 'lhe':
+                self.cross, self.error = 0, 0
         self.cross *= br
         self.error *= br
         
-        
         # modify the cross-section in the init block of the banner
-        self.banner.scale_init_cross(self.branching_ratio)
-                    
-        
+        if not self.options['cross_section']:
+            self.banner.scale_init_cross(self.branching_ratio)
+        else:
+            
+            if self.options['input_format'] in ['lhe_no_banner','hepmc'] and 'init' not in self.banner:
+                self.cross = sum(self.options['cross_section'].values())
+                self.error = 0
+                self.branching_ratio = 1
+            else:  
+                self.banner.modify_init_cross(self.options['cross_section'])
+                new_cross, new_error =   self.banner.get_cross(witherror=True)
+                self.branching_ratio = new_cross / self.cross
+                self.cross = new_cross   
+                self.error = new_error
+
         # 3. Merge the various file together.
-        output_lhe = lhe_parser.EventFile(orig_lhe.name.replace('.lhe', '_decayed.lhe.gz'), 'w')
-        self.banner.write(output_lhe, close_tag=False)
+        if self.options['input_format'] == 'hepmc':
+            name = orig_lhe.name.replace('.hepmc', '_decayed.lhe')
+            if not name.endswith('.gz'):
+                name = '%s.gz' % name
+            
+            output_lhe = lhe_parser.EventFile(name, 'w')
+        else:
+            name = orig_lhe.name.replace('.lhe', '_decayed.lhe')
+            if not name.endswith('.gz'):
+                name = '%s.gz' % name
+            output_lhe = lhe_parser.EventFile(name, 'w')
+        try:
+            self.banner.write(output_lhe, close_tag=False)
+        except Exception:
+            if self.options['input_format'] == 'lhe':
+                raise
         
         # initialise object which store not use event due to wrong helicity
         bufferedEvents_decay = {}
@@ -856,8 +968,9 @@
         start = time.time()
         counter = 0
         orig_lhe.seek(0)
+
         for event in orig_lhe:
-            if counter and counter % 1000 == 0 and float(str(counter)[1:]) ==0:
+            if counter and counter % 100 == 0 and float(str(counter)[1:]) ==0:
                 print "decaying event number %s [%s s]" % (counter, time.time()-start)
             counter +=1
             
@@ -866,10 +979,17 @@
             particles = [p for p in event if int(p.status) == 1.0]
             random.shuffle(particles)
             ids = [particle.pid for particle in particles]
+            br = 1 #br for that particular events (for special/weighted case)
+            hepmc_output = lhe_parser.Event() #for hepmc case: collect the decay particle
             for i,particle in enumerate(particles):
+                #misc.sprint(i, particle.pdg, particle.pid)
+                #misc.sprint(self.final_state, evt_decayfile)
                 # check if we need to decay the particle 
-                if particle.pdg not in self.final_state or particle.pdg not in evt_decayfile:
-                    continue # nothing to do for this particle
+                if self.final_state and particle.pdg not in self.final_state:
+                    continue # nothing to do for this particle
+                if particle.pdg not in evt_decayfile:
+                    continue # nothing to do for this particle
+                
                 # check how the decay need to be done
                 nb_decay = len(evt_decayfile[particle.pdg])
                 if nb_decay == 0:
@@ -888,14 +1008,28 @@
                     cumul = 0
                     for j,events in evt_decayfile[particle.pdg].items():
                         cumul += events.cross
-                        if r < cumul:
+                        if r <= cumul:
                             decay_file = events
                             decay_file_nb = j
-                        else:
                             break
-                        
+                    else:
+                        # security for numerical accuracy issue... (unlikely but better safe)
+                        if (cumul-tot)/tot < 1e-5:
+                            decay_file = events
+                            decay_file_nb = j
+                        else:
+                            misc.sprint(j,cumul, events.cross, tot, (tot-cumul)/tot)
+                            raise Exception
+                
+                if self.options['new_wgt'] == 'BR':
+                    tot_width = float(self.banner.get('param', 'decay', abs(pdg)).value)
+                    if tot_width:
+                        br = decay_file.cross / tot_width
                 # ok start the procedure
-                helicity = particle.helicity
+                if hasattr(particle,'helicity'):
+                    helicity = particle.helicity
+                else:
+                    helicity = 9
                 bufferedEvents = bufferedEvents_decay[particle.pdg][decay_file_nb]
                 
                 # now that we have the file to read. find the associate event
@@ -910,8 +1044,8 @@
                         except StopIteration:
                             # check how far we are
                             ratio = counter / nb_event 
-                            needed = 1.05 * to_decay[particle.pdg]/ratio - counter
-                            needed = min(1000, max(needed, 1000))
+                            needed = 1.05 * to_decay[particle.pdg] - counter
+                            needed = min(100000, max(needed, 6000))
                             with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
                                 new_file = self.generate_events(particle.pdg, needed, mg5, [decay_file_nb])
                             evt_decayfile[particle.pdg].update(new_file)
@@ -927,14 +1061,36 @@
                             # only add to the buffering if the buffer is not too large
                             bufferedEvents[helicity].append(decay)
                 # now that we have the event make the merge
-                particle.add_decay(decay)
+                if self.options['input_format'] != 'hepmc':
+                    particle.add_decay(decay)
+                else:
+                    if len(hepmc_output) == 0:
+                        hepmc_output.append(lhe_parser.Particle(event=hepmc_output))
+                        hepmc_output[0].color2 = 0
+                    decayed_particle = lhe_parser.Particle(particle, hepmc_output)
+                    hepmc_output.append(decayed_particle)
+                    decayed_particle.add_decay(decay)
             # change the weight associate to the event
-            event.wgt *= self.branching_ratio
-            wgts = event.parse_reweight()
-            for key in wgts:
-                wgts[key] *= self.branching_ratio
-            # all particle have been decay if needed
-            output_lhe.write(str(event))
+            if self.options['new_wgt'] == 'cross-section':
+                event.wgt *= self.branching_ratio
+                br = self.branching_ratio
+            else:
+                event.wgt *= br
+                
+            if self.options['input_format'] != 'hepmc':
+                wgts = event.parse_reweight()
+                for key in wgts:
+                    wgts[key] *= br
+                # all particle have been decay if needed
+                output_lhe.write(str(event))
+            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:
+                raise Exception
         output_lhe.write('</LesHouchesEvents>\n')        
                     
     
@@ -971,7 +1127,12 @@
            cumul allow to merge all the definition in one run (add process)
                  to generate events according to cross-section
         """
-        
+        if not hasattr(self, 'me_int'):
+            self.me_int = {}
+            
+        
+        
+        nb_event = int(nb_event) # in case of hepmc request the nb_event is not an integer
         if cumul:
             width = 0.
         else:   
@@ -1003,31 +1164,36 @@
                 options = dict(mg5.options)
                 if self.options['ms_dir']:
                     # we are in gridpack mode -> create it
-                    me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
-                                            decay_dir), options=options)
-                    me5_cmd.options["automatic_html_opening"] = False
-                    me5_cmd.options["madanalysis5_path"] = None
-                    me5_cmd.options["madanalysis_path"] = None
-                    me5_cmd.allow_notification_center = False
-                    try:
-                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
-                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
-                    except Exception,error:
-                        logger.debug(error)
-                        pass            
+                    if decay_dir in self.me_int:
+                        me5_cmd = self.me_int[decay_dir]
+                    else:
+                        me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
+                                                decay_dir), options=options)
+                        me5_cmd.options["automatic_html_opening"] = False
+                        me5_cmd.options["madanalysis5_path"] = None
+                        me5_cmd.options["madanalysis_path"] = None
+                        me5_cmd.allow_notification_center = False
+                        try:
+                            os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
+                            os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
+                        except Exception,error:
+                            logger.debug(error)
+                            pass 
+                        self.me_int[decay_dir] = me5_cmd
+
                     if self.options["run_card"]:
-                        run_card = self.options["run_card"]
+                        run_card = self.run_card
                     else:
                         run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))                        
                     
-                    run_card["iseed"] = self.seed
+                    run_card["iseed"] = self.options['seed']
                     run_card['gridpack'] = True
                     run_card['systematics_program'] = 'False'
                     run_card['use_syst'] = False
                     run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
                     param_card = self.banner['slha']
                     open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
-                    self.seed += 1
+                    self.options['seed'] += 1
                     # actually creation
                     me5_cmd.exec_cmd("generate_events run_01 -f")
                     if output_width:
@@ -1041,26 +1207,35 @@
                     misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir)
             
             # Now generate the events
-
             if not self.options['ms_dir']:
-                me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
-                                                decay_dir), options=mg5.options)
-                me5_cmd.options["automatic_html_opening"] = False
-                me5_cmd.options["automatic_html_opening"] = False
-                me5_cmd.options["madanalysis5_path"] = None
-                me5_cmd.options["madanalysis_path"] = None
-                me5_cmd.allow_notification_center = False
-                try:
-                    os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
-                    os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
-                except Exception,error:
-                    logger.debug(error)
-                    pass                    
+                if decay_dir in self.me_int:
+                        me5_cmd = self.me_int[decay_dir]
+                else:
+                    me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
+                                                    decay_dir), options=mg5.options)
+                    me5_cmd.options["automatic_html_opening"] = False
+                    me5_cmd.options["automatic_html_opening"] = False
+                    me5_cmd.options["madanalysis5_path"] = None
+                    me5_cmd.options["madanalysis_path"] = None
+                    me5_cmd.allow_notification_center = False
+                    try:
+                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
+                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
+                    except Exception,error:
+                        logger.debug(error)
+                        pass                 
+                    self.me_int[decay_dir] = me5_cmd
                 if self.options["run_card"]:
-                    run_card = 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)
+                # Handle the banner of the output file
+                if not self.seed:
+                    self.seed = random.randint(0, int(30081*30081))
+                    self.do_set('seed %s' % self.seed)
+                    logger.info('Will use seed %s' % self.seed)
+                    self.history.insert(0, 'set seed %s' % self.seed)
                 run_card["iseed"] = self.seed
                 run_card["systematics_program"] = 'None'
                 run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
@@ -1140,8 +1315,7 @@
                 if particle.status == 1 and particle.pdg in asked_to_decay:
                     # final state and tag as to decay
                     to_decay[particle.pdg] += 1
-        #misc.sprint(to_decay)
-        #misc.sprint("import the mode -> temporary with logging")   
+
         with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
             mg5 = self.mg5cmd
             if not self.model:
@@ -1436,10 +1610,6 @@
         full_event = lhe_parser.Event(str(production))
         full_event = full_event.add_decays(decays)
         full_me = self.calculate_matrix_element(full_event)
-        #misc.sprint(full_event)
-        #misc.sprint([p.pdg for p in production])
-        #misc.sprint([p.pdg for p in full_event])
-        #misc.sprint(full_me, production_me, decay_me)
         return full_event, full_me/(production_me*decay_me)
         
         
@@ -1457,7 +1627,6 @@
             init = (-init[0],)
             final = tuple(-i for i in final)
             tag = (init, final)
-            misc.sprint([k for k in  self.all_me.keys() if len(k[0])==1])
             orig_order = self.all_me[tag]['order']
             
         pdir = self.all_me[tag]['pdir']
@@ -1520,7 +1689,6 @@
                 
             rwgt_interface.ReweightInterface.get_LO_definition_from_NLO()        
         
-        
 
 
 if __name__ == '__main__':

=== modified file 'MadSpin/madspin'
--- MadSpin/madspin	2013-11-12 23:52:05 +0000
+++ MadSpin/madspin	2018-12-20 08:18:55 +0000
@@ -49,6 +49,7 @@
 if len(args) == 0:
     args = ''
 
+
 import subprocess
 
 # Check if optimize mode is (and should be) activated
@@ -114,23 +115,6 @@
     pass
 import MadSpin.interface_madspin as interface
 
-try:
-    cmd_line = interface.MadSpinInterface()
-    cmd_line.use_rawinput = True
-    cmd_line.cmdloop()
-except:
-    pass
-try:
-    cmd_line.exec_cmd('quit all', printcmd=False)
-    readline.set_history_length(100)
-    if not os.path.exists(os.path.join(os.environ['HOME'], '.mg5')):
-        os.mkdir(os.path.join(os.environ['HOME'], '.mg5'))
-    readline.write_history_file(history_file)
-except Exception, error:
-    pass
-
-sys.exit(0)
-
 # Call the cmd interface main loop
 try:
     if options.file or args:
@@ -143,12 +127,14 @@
             cmd_line = interface.MasterCmdWeb()
             cmd_line.debug_output = os.path.join(os.path.dirname(input_file),'generation.log')
             cmd_line.use_rawinput = False
-            cmd_line.run_cmd('import ' + input_file)
+            cmd_line.import_command_file(input_file)
             cmd_line.run_cmd('quit')
         else:
+            print "using input+file", input_file
             cmd_line = interface.MadSpinInterface()
             cmd_line.use_rawinput = False
-            cmd_line.run_cmd('import ' + input_file)
+            cmd_line.haspiping = False
+            cmd_line.import_command_file(input_file)
             cmd_line.run_cmd('quit')
     else:
         # Interactive mode
@@ -160,8 +146,20 @@
             cmd_line = interface.MasterCmdWeb()
             cmd_line.cmdloop()
         else:
-            cmd_line = interface.MasterCmd(mgme_dir = options.mgme_dir)
-            cmd_line.cmdloop()
+            try:
+                cmd_line = interface.MadSpinInterface()
+                cmd_line.use_rawinput = True
+                cmd_line.cmdloop()
+            except:
+                pass
+        try:
+            cmd_line.exec_cmd('quit all', printcmd=False)
+            readline.set_history_length(100)
+            if not os.path.exists(os.path.join(os.environ['HOME'], '.mg5')):
+                os.mkdir(os.path.join(os.environ['HOME'], '.mg5'))
+            readline.write_history_file(history_file)
+        except Exception, error:
+            pass
 except KeyboardInterrupt:
     print 'writting history and quit on KeyboardInterrupt' 
     pass

=== modified file 'MadSpin/src/driver.f'
--- MadSpin/src/driver.f	2017-08-03 18:49:25 +0000
+++ MadSpin/src/driver.f	2018-12-20 08:18:55 +0000
@@ -340,10 +340,10 @@
            if (jac.lt.0d0) then
              counter2=counter2+1 
              counter3=counter3+1 
-c             if (counter3.gt.500) then
-c               write(*,*) "500_pts_failed_stop_executation"
-c               stop
-c             endif
+             if (counter3.gt.500*8*100) then
+               write(*,*) "500_pts_failed_stop_executation"
+               stop
+             endif
              if (counter2.ge.8) then ! use another topology to generate PS points
                do k=1,n_max_cg
                  amp2(k)=0d0
@@ -960,7 +960,7 @@
       include 'nexternal_prod.inc'
       !include 'run.inc'
 c arguments
-      double precision jac,x(36),p(0:3,nexternal), p_prod(0:3,nexternal)
+      double precision jac,x(36),p(0:3,nexternal), p_prod(0:3,nexternal_prod)
       integer itree(2,-nexternal:-1)
       double precision qmass(-nexternal:0),qwidth(-nexternal:0)
 c common

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2018-04-25 21:20:23 +0000
+++ Template/LO/Cards/run_card.dat	2018-12-20 08:18:55 +0000
@@ -12,7 +12,7 @@
 #   mind the format:   value    = variable     ! comment             *
 #                                                                    *
 #   To display more options, you can type the command:               *
-*      update full_run_card                                          *
+#      update full_run_card                                          *
 #*********************************************************************
 #
 #*******************                                                 

=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2018-04-20 08:38:20 +0000
+++ Template/LO/SubProcesses/genps.f	2018-12-20 08:18:55 +0000
@@ -41,7 +41,7 @@
       integer mincfig,maxcfig                  !Range of configurations
       integer invar
       double precision wgt                     !(input and output)
-      double precision x(maxdim),p(maxdim)     !x,p (output) [p(0:3,nexternal)]
+      double precision x(*),p(*)     !x,p (output) [p(0:3,nexternal)]
 c
 c     Local
 c
@@ -104,7 +104,7 @@
 c
       integer iconfig,mincfig,maxcfig,invar
       double precision p1(0:3,nexternal+1)
-      double precision x(maxinvar)
+      double precision x(*)
       double precision wgt
 c
 c     Local
@@ -191,8 +191,6 @@
 c
       include 'configs.inc'
       data firsttime/.true./
-      integer isym(0:100)
-c      data isym /2,1,5,27,42,47,0,0,0,0,0/
       data jfig/1/
 c-----
 c  Begin Code
@@ -215,11 +213,6 @@
          write(*,'(a,12e10.3)') ' Masses:',(m(i),i=1,nparticles)
       endif                          !First_time
 
-      if (.false.) then
-         iconfig = isym(jfig)
-         jfig = jfig+1
-         if (jfig .gt. isym(0)) jfig=1      
-      endif
       this_config = iconfig             !Pass iconfig to amplitude routine
 C
 C     Get fraction of beam energy if pdf's are used
@@ -279,19 +272,18 @@
          cm_rap=.5d0*dlog((p0+p3)/(p0-p3))
          set_cm_rap=.true.
 c        Set shat
-         s(-nbranch) = m2**2+2*xbk(1)*ebeam(1) *
-     $                 (ebeam(2)+sqrt(ebeam(2)**2-m2**2))
+         s(-nbranch) = x(ndim)*stot         
       elseif (abs(lpp(2)) .ge. 1) then
          call sample_get_x(sjac,x(ndim),ndim,mincfig,0d0,1d0)
          xbk(2) = x(ndim)
+
 c        Set CM rapidity for use in the rap() function
          p0=ebeam(1)+xbk(2)*ebeam(2)
          p3=sqrt(ebeam(1)**2-m1**2)-xbk(2)*ebeam(2)
          cm_rap=.5d0*dlog((p0+p3)/(p0-p3))
          set_cm_rap=.true.
 c        Set shat
-         s(-nbranch) = m1**2+2*(ebeam(1)+sqrt(ebeam(1)**2-m1**2))
-     $                 * xbk(2)*ebeam(2)
+         s(-nbranch) =  x(ndim)*stot         
       else
 c        Set CM rapidity for use in the rap() function
          p0=ebeam(1) + ebeam(2)
@@ -576,7 +568,7 @@
             if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
             if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
             if (mass_ion(1).ge.0d0) m1 = mass_ion(1)
-            if (mass_ion(2).ge.0d0) m1 = mass_ion(2)
+            if (mass_ion(2).ge.0d0) m2 = mass_ion(2)
             if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1
             if(ebeam(2).lt.m2.and.lpp(2).ne.9) ebeam(2)=m2
             pi1(0)=ebeam(1)
@@ -598,7 +590,6 @@
          maxcfig=iconfig
          call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
          maxwgt=0d0
-c         write(*,'(a,12i4)') 'Summing configs',(isym(i),i=1,isym(0))
          nparticles   = nexternal
          nfinal       = nparticles-nincoming
          nbranch      = nparticles-2
@@ -617,8 +608,6 @@
 c
 c         do i=1,mapconfig(0)
          do i=mincfig,maxcfig
-c         do k=1,isym(0)
-c            i = isym(k)
             write(*,'(15i4)') i,(minvar(j,i),j=1,ndim)
             do j=1,ndim
                ipole = minvar(j,i)
@@ -663,7 +652,7 @@
 c      double precision spole(-max_branch:0),swidth(-max_branch:0)
       double precision jac,pswgt
       integer nbranch
-      double precision x(40) ! ja 3/2/11 21->40 after strange segfault
+      double precision x(*) ! ja 3/2/11 21->40 after strange segfault
 c
 c     Local
 c

=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2018-03-11 12:29:51 +0000
+++ Template/LO/SubProcesses/reweight.f	2018-12-20 08:18:55 +0000
@@ -1,3 +1,21 @@
+c     for cross-checking change in this file. 
+c     here is a minimal list of process that we have to 
+c     test
+c    
+c     SM
+c     ----
+c     p p > t t~ (up to 2jet)
+c     p p > w+ (up to 3 jet)
+c     p p > j j w+ (ordering seems important) 
+c     p p > z t t~ j j  (no MLM needed)
+c
+c
+c     HEFT
+c     ----
+c     p p > h j b b~
+c     q q > a a g q q
+c     g g > h g q q
+
       double precision function gamma(q0)
 c**************************************************
 c   calculates the branching probability
@@ -226,13 +244,15 @@
       integer iddgluon, iddother, idgluon, idother
       logical isqcd
       external isqcd
+      integer get_color
+      external get_color 
 
       idmo=ipdg(imo)
       idda1=ipdg(ida1)
       idda2=ipdg(ida2)
 
       if (btest(mlevel,4)) then
-        write(*,*) ' updating ipart for: ',ida1,ida2,' -> ',imo
+        write(*,*) 'updating ipart for: ',ida1,ida2,' -> ',imo
       endif
 
       if (btest(mlevel,4)) then
@@ -274,7 +294,6 @@
      $        ' (',ipdg(imo),')'
          return
       endif        
-
 c     FS clustering
 c     Transmit parton PDG code for parton vertex
       if(isjet(idmo)) then
@@ -353,14 +372,26 @@
 c     quark -> gluon-quark or Z-quark or h-quark or W-quark
         ipart(1,imo)=ipart(1,ida2)
         ipart(2,imo)=0
-      else
+      else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda1)).eq.3.and.get_color(idda2).eq.1) then
+c       exotic q > q' Scalar         
+        ipart(1,imo)=ipart(1,ida1)
+        ipart(2,imo)=0
+      else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda2)).eq.3.and.get_color(idda1).eq.1) then
+c       exotic q > Scalar q'
+        ipart(1,imo)=ipart(1,ida2)
+        ipart(2,imo)=0
+      else if (get_color(idmo).eq.1) then
 c     Color singlet
          ipart(1,imo)=ipart(1,ida1)
          ipart(2,imo)=ipart(1,ida2)
+      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"
+         stop 3
       endif
       
       if (btest(mlevel,4)) then
-        write(*,*) ' -> ',(ipart(i,imo),i=1,2),' (',ipdg(imo),')'
+        write(*,*) 'XY -> ',(ipart(i,imo),i=1,2),' (',ipdg(imo),')'
       endif
 
       return
@@ -522,7 +553,8 @@
       common/to_stot/stot,m1,m2
 
 C   local variables
-      integer i, j, idi, idj, k
+      integer i, j, idi, idj, k,m
+      integer get_color
       real*8 PI
       parameter( PI = 3.14159265358979323846d0 )
       integer iforest(2,-max_branch:-1,lmaxconfigs)
@@ -545,6 +577,7 @@
       logical chclusold,fail,increasecode
       save chclusold
       integer tmpindex
+      integer pdgm, pdgid1, pdgid2
 
       logical isqcd,isjet,isparton,cluster,isjetvx,is_octet
       integer ifsno
@@ -676,6 +709,7 @@
 c     increasecode gives whether we should increase jcode at next vertex
       increasecode=.false.
       do n=1,nexternal-2
+c         write(*,*) 'QCD jet status (before n= ',n,'):',(iqjets(i),i=3,nexternal)
         do i=1,2 ! index of the child in the interaction
           do j=1,2 ! j index of the beam
             if(idacl(n,i).eq.ibeam(j))then
@@ -746,8 +780,16 @@
 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)
-                 if (ipart(2,imocl(n)).eq.0) then ! q a > q case
-                    iqjets(ipart(1,imocl(n)))=0
+                 ! 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
+                     continue
+                 elseif (ipart(2,imocl(n)).eq.0) then ! q a > q case
+                       iqjets(ipart(1,imocl(n)))=0
                  else ! octet. want to be sure that both are tagged as jet before removing one
                     ! this prevent that both are removed in case of g > g h , g > q1 q2, q1 > a q1.
                     ! at least one of the two should be kept as jet
@@ -791,7 +833,57 @@
            goodjet(imocl(n))=
      $          (isjet(ipdgcl(imocl(n),igraphs(1),iproc)).and.
      $          goodjet(idacl(n,1)).and.goodjet(idacl(n,2)))
-        endif
+            
+c       check case with g > g g
+c       where the hardest gluon is not goodjet but the other is.
+c       in that case change ipart(1,) of the mother gluon
+c             pure QCD jet
+c             need to take care of the following case:
+c                                      tttttttttt
+c                      gggggggggggggggg  
+c             ggggggggg                tttttttttt 
+c                      gggggggg
+c                        
+c             in that case the up gluon can be tag as the hardest one
+c             but this one is also lead to no QCD one.
+c             so in that case we have to change ipart(1) to the sofest gluon
+           pdgm   = ipdgcl(imocl(n),igraphs(1),iproc)
+           pdgid1 = ipdgcl(idacl(n,1),igraphs(1),iproc)
+           pdgid2 = ipdgcl(idacl(n,2),igraphs(1),iproc)
+           if (is_octet(pdgm).and.is_octet(pdgid1).and.is_octet(pdgid2))then
+c              write(*,*) 'pure QCD vertex (2)'
+c              write(*,*) pdgm , '>', pdgid1,' ', pdgid2
+c              write(*,*) 'ipart', ipart(1,imocl(n)), ipart(2,imocl(n))
+c              write(*,*) 'id', imocl(n), idacl(n,1),idacl(n,2)
+c              write(*,*) 'ipart',ipart(1,imocl(n)),'/', ipart(2,imocl(n)), ipart(1,idacl(n,1)),'/', ipart(2,idacl(n,1)),
+c     $             ipart(1,idacl(n,2)),'/', ipart(2,idacl(n,2))
+c              write(*,*) 'googjet', goodjet(idacl(n,1)),goodjet(idacl(n,2))
+              if (ipart(1,imocl(n)).eq.ipart(1, idacl(n,1))) then
+                  if (.not.goodjet(idacl(n,1)).and.goodjet(idacl(n,2))) then
+c                     write(*,*) 'ggg with hard jet set a QED the second jet lead to', ipart(1,idacl(n,2)), ipart(2,idacl(n,2))
+                     do m =n_max_cl,n,-1
+                        if(ipart(1,m).eq.ipart(1,imocl(n)).and.ipart(2,m).eq.ipart(2,imocl(n)))then
+                           ipart(1,m) = ipart(1,idacl(n,2))
+                           ipart(2,m) = ipart(2,idacl(n,2))
+                        endif
+                     enddo
+                  endif
+               else
+                  if (.not.goodjet(idacl(n,2)).and.goodjet(idacl(n,1))) then
+c                     write(*,*) 'ggg with hard jet set a QED the second jet lead to', ipart(1,idacl(n,1)), ipart(2,idacl(n,1))
+                     do m =n_max_cl,n,-1
+                        if(ipart(1,m).eq.ipart(1,imocl(n)).and.ipart(2,m).eq.ipart(2,imocl(n)))then
+                           ipart(1,m) = ipart(1,idacl(n,1))
+                           ipart(2,m) = ipart(2,idacl(n,1))
+                        endif
+                     enddo
+                  endif
+               endif
+            endif
+
+
+
+       endif
       enddo
 
       if (btest(mlevel,4))then
@@ -836,8 +928,8 @@
             endif
          enddo
          njetstore(iconfig)=njets
-         if (btest(mlevel,4))
-     $        write(*,*) 'Storing jets: ',(iqjetstore(i,iconfig),i=1,njets)
+        if (btest(mlevel,4))
+     $       write(*,*) 'Storing jets: ',(iqjetstore(i,iconfig),i=1,njets)
 c     Recluster without requiring chcluster
          goto 100
       else

=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2017-06-22 08:03:29 +0000
+++ Template/LO/SubProcesses/unwgt.f	2018-12-20 08:18:55 +0000
@@ -513,6 +513,10 @@
 
       double precision stot,m1,m2
       common/to_stot/stot,m1,m2
+
+      INTEGER IMIRROR
+      INTEGER IPROC
+      COMMON/TO_MIRROR/IMIRROR, IPROC
 c
 c     Data
 c
@@ -641,6 +645,13 @@
             pb(4,isym(j,jsym))=pmass(j)
          enddo
       endif
+
+      if (IMIRROR.eq.2.and.pmass(1).ne.pmass(2)) then
+c        Note that in this context isym(1,jsym) should never be "2" since the mass differ 
+         pb(4,isym(1,jsym))=pmass(2)
+         pb(4,isym(2,jsym))=pmass(1)
+      endif
+
 c
 c     Add info on resonant mothers
 c

=== added file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f'
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f	1970-01-01 00:00:00 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f	2018-12-20 08:18:55 +0000
@@ -0,0 +1,880 @@
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+      integer i,l
+      character*6 cc(2)
+      data cc/'|T@NLO','|T@LO '/
+
+c     Also specific perturbative orders can be directly plotted, adding for examples
+c     the following further entries in the variable data
+c     $        ' |T@QCD4QED0 ',
+c     $        ,' |T@QCD2QED2 ',' |T@QCD0QED4 ',' |T@QCD6QED0 '
+c     $        ,' |T@QCD4QED2 ',' |T@QCD2QED4 ',' |T@QCD0QED6 '
+c     $        ,' |T@QCD8QED0 ',' |T@QCD6QED2 ',' |T@QCD4QED4 '
+c     $        ,' |T@QCD2QED6 ',' |T@QCD0QED8 '
+c     
+c     See also line 376 in this file 
+      call HwU_inithist(nwgt,weights_info)
+      do i=1,2
+      l=(i-1)*55
+c transverse momenta
+      call HwU_book(l+ 1,'total rate                 '//cc(i),1,0.5d0,1.5d0)
+      call HwU_book(l+ 2,'1st charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+ 3,'2nd charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+ 4,'3rd charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+ 5,'4th charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+ 6,'electron log pT            '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+ 7,'positron log pT            '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+ 8,'mu-plus log pT             '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+ 9,'mu-minus log pT            '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+14,'epem log pT                '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+15,'mupmum log pT              '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+16,'1st OSSF lep pair log pT   '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+17,'2nd OSSF lep pair log pT   '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+18,'epve log pT                '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+19,'mumvm log pT               '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+20,'1st SF lep-neu pair log pT '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+21,'2nd SF lep-neu pair log pT '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+22,'epemmupmum log pT          '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+23,'epvemumvm log pT           '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+24,'top quark log pT           '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+25,'anti-top quark log pT      '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+26,'ttbar log pT               '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+27,'log missing Et (neutrinos) '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+28,'1st higgs log pT           '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+29,'2nd higgs log pT           '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+30,'higgs-pair log pT          '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+31,'1st V-boson log pT         '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+32,'2nd V-boson log pT         '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+33,'3rd V-boson log pT         '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+34,'1st jet log pT             '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+35,'2nd jet log pT             '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+36,'3rd jet log pT             '//cc(i),25,-0.2d0,3.8d0)
+c
+c invariant masses
+      call HwU_book(l+37,'epem log M                 '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+38,'mupmum log M               '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+39,'1st OSSF lep pair log M    '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+40,'2nd OSSF lep pair log M    '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+41,'epve log M                 '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+42,'mumvm log M                '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+43,'1st SF lep-neu pair log M  '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+44,'2nd SF lep-neu pair log M  '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+45,'epemmupmum log M           '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+46,'epvemumvm log M            '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+47,'ttbar log M                '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+48,'higgs-pair log M           '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+49,'j1j2 log M                 '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+50,'j2j3 log M                 '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+51,'j1j3 log M                 '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+52,'j1j2j3 log M               '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+53,'3w log M                   '//cc(i),25,-0.2d0,3.8d0)
+c
+c HT
+      call HwU_book(l+54,'log HT (partons)           '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+55,'log HT (reconstructed)     '//cc(i),25,-0.2d0,3.8d0)
+      enddo
+      return
+      end
+
+
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(dummy)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      double precision dummy
+      call HwU_write_file
+      return                
+      end
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      implicit none
+      include 'nexternal.inc'
+      include 'cuts.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      double precision p(0:4,nexternal)
+      double precision wgts(*)
+      integer ibody
+      integer i,j,k,kk,l
+      double precision www,pQCD(0:3,nexternal),palg,rfj,sycut,yjmax
+     $     ,pjet(0:3,nexternal),tmp,ptlep(4),ptem,ptep,ptmm,ptmp,ptepem
+     $     ,ptmpmm,ptepve,ptmmvm,ptepemmpmm,ptepvemmvm,ptt,ptat,pttt
+     $     ,etmiss,pth(2),pthh,ptv(3),ptjet(3),Mepem,Mmpmm,Mepve,Mmmvm
+     $     ,Mepemmpmm,Mepvemmvm,Mtt,Mhh,Mj1j2,Mj1j3,Mj2j3,Mj1j2j3,Mvvv
+     $     ,HTparton,HTreco,p_reco(0:4,nexternal)
+      integer nQCD,jet(nexternal),njet,itop,iatop,iem,iep,imp,imm,ive
+     $     ,ivm,iv1,iv2,iv3,ih1,ih2,il,ipdg_reco(nexternal)
+      double precision getptv4,getptv4_2,getptv4_4,getinvm4_2
+     $     ,getinvm4_3,getinvm4_4,l10
+      external getptv4,getptv4_2,getptv4_4,getinvm4_2,getinvm4_3
+     $     ,getinvm4_4,l10
+      integer orders_tag_plot
+      common /corderstagplot/ orders_tag_plot
+c First, try to recombine photons with leptons      
+      if (.not.quarkphreco) then
+         write (*,*) 'quark-photon recombination is turned off. '/
+     $        /'Do need it'
+         stop
+      endif
+      if (.not. lepphreco) then
+         write (*,*) 'lepton-photon recombination is turned off. '
+     $        //'Do need it.'
+         stop
+      endif
+      call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
+     $                       p, iPDG, p_reco, iPDG_reco)
+
+c Put all (light) QCD partons(+photon) in momentum array for jet clustering.
+      nQCD=0
+      do j=nincoming+1,nexternal
+         if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 .or.
+     $        ipdg_reco(j).eq.22)then
+            nQCD=nQCD+1
+            do i=0,3
+               pQCD(i,nQCD)=p_reco(i,j)
+            enddo
+         endif
+      enddo
+      
+C---CLUSTER THE EVENT
+      palg  = jetalgo
+      rfj   = jetradius
+      sycut = ptj
+      yjmax = etaj
+      
+c******************************************************************************
+c     call FASTJET to get all the jets
+c
+c     INPUT:
+c     input momenta:               pQCD(0:3,nexternal), energy is 0th component
+c     number of input momenta:     nQCD
+c     radius parameter:            rfj
+c     minumum jet pt:              sycut
+c     jet algorithm:               palg, 1.0=kt, 0.0=C/A, -1.0 = anti-kt
+c
+c     OUTPUT:
+c     jet momenta:                             pjet(0:3,nexternal), E is 0th cmpnt
+c     the number of jets (with pt > SYCUT):    njet
+c     the jet for a given particle 'i':        jet(i),   note that this is
+c     the particle in pQCD, which doesn't necessarily correspond to the particle
+c     label in the process
+c
+      call amcatnlo_fastjetppgenkt_etamax(pQCD,nQCD,rfj,sycut,yjmax,palg
+     $     ,pjet,njet,jet)
+c
+c******************************************************************************
+
+c Look for the other physics objects
+      itop=0
+      iatop=0
+      iem=0
+      iep=0
+      iem=0
+      iep=0
+      imp=0
+      imm=0
+      ive=0
+      ivm=0
+      iv1=0
+      iv2=0
+      iv3=0
+      ih1=0
+      ih2=0
+      do i=1,nexternal
+         if (ipdg_reco(i).eq.6) then
+            itop=i
+         elseif(ipdg_reco(i).eq.-6) then
+            iatop=i
+         elseif(ipdg_reco(i).eq.11) then
+            iem=i
+         elseif(ipdg_reco(i).eq.13) then
+            imm=i
+         elseif(ipdg_reco(i).eq.-11) then
+            iep=i
+         elseif(ipdg_reco(i).eq.-13) then
+            imp=i
+         elseif(abs(ipdg_reco(i)).eq.12) then
+            ive=i
+         elseif(abs(ipdg_reco(i)).eq.14) then
+            ivm=i
+         elseif(abs(ipdg_reco(i)).eq.24 .or. ipdg_reco(i).eq.23) then
+            if (iv1.eq.0) then
+               iv1=i
+            else
+               if (iv2.eq.0) then
+                  iv2=i
+               else
+                  if(iv3.eq.0) then
+                     iv3=i
+                  else
+                     write (*,*) 'too many vector bosons'
+                     stop
+                  endif
+               endif
+            endif
+         elseif(abs(ipdg_reco(i)).eq.25) then
+            if (ih1.eq.0) then
+               ih1=i
+            else
+               if (ih2.eq.0) then
+                  ih2=i
+               else
+                  write (*,*) 'too many higgs bosons'
+                  stop
+               endif
+            endif
+         endif
+      enddo
+      if (itop.ne.0) ptt=getptv4(p_reco(0,itop))
+      if (iatop.ne.0) ptat=getptv4(p_reco(0,iatop))
+      if (itop.ne.0 .and. iatop.ne.0) then
+         pttt=getptv4_2(p_reco(0,itop),p_reco(0,iatop))
+         Mtt=getinvm4_2(p_reco(0,itop),p_reco(0,iatop))
+      endif
+      if (iem.ne.0) ptem=getptv4(p_reco(0,iem))
+      if (iep.ne.0) ptep=getptv4(p_reco(0,iep))
+      if (imm.ne.0) ptmm=getptv4(p_reco(0,imm))
+      if (imp.ne.0) ptmp=getptv4(p_reco(0,imp))
+      if (iem.ne.0 .and. iep.ne.0) then
+         ptepem=getptv4_2(p_reco(0,iem),p_reco(0,iep))
+         Mepem=getinvm4_2(p_reco(0,iem),p_reco(0,iep))
+      endif
+      if (imm.ne.0 .and. imp.ne.0) then
+         ptmpmm=getptv4_2(p_reco(0,imm),p_reco(0,imp))
+         Mmpmm=getinvm4_2(p_reco(0,imm),p_reco(0,imp))
+      endif
+      if (iep.ne.0 .and. ive.ne.0) then
+         ptepve=getptv4_2(p_reco(0,iep),p_reco(0,ive))
+         Mepve=getinvm4_2(p_reco(0,iep),p_reco(0,ive))
+      endif
+      if (imm.ne.0 .and. ivm.ne.0) then
+         ptmmvm=getptv4_2(p_reco(0,imm),p_reco(0,ivm))
+         Mmmvm=getinvm4_2(p_reco(0,imm),p_reco(0,ivm))
+      endif
+      if (iem.ne.0 .and. iep.ne.0 .and. imm.ne.0 .and. imp.ne.0) then
+         ptepemmpmm=getptv4_4(p_reco(0,iem),p_reco(0,iep),p_reco(0,imm),p_reco(0,imp))
+         Mepemmpmm=getinvm4_4(p_reco(0,iem),p_reco(0,iep),p_reco(0,imm),p_reco(0,imp))
+      endif      
+      if (ive.ne.0 .and. iep.ne.0 .and. imm.ne.0 .and. ivm.ne.0) then
+         ptepvemmvm=getptv4_4(p_reco(0,iep),p_reco(0,ive),p_reco(0,imm),p_reco(0,ivm))
+         Mepvemmvm=getinvm4_4(p_reco(0,iep),p_reco(0,ive),p_reco(0,imm),p_reco(0,ivm))
+      endif
+
+      il=0
+      if (iem.ne.0) then
+         il=il+1
+         ptlep(il)=ptem
+      endif
+      if (iep.ne.0) then
+         il=il+1
+         ptlep(il)=ptep
+      endif
+      if (imm.ne.0) then
+         il=il+1
+         ptlep(il)=ptmm
+      endif
+      if (imp.ne.0) then
+         il=il+1
+         ptlep(il)=ptmp
+      endif
+      if (il.gt.1) then
+         do i=1,il-1
+            do j=1,il-i
+               if (ptlep(j).lt.ptlep(j+1)) then
+                  tmp=ptlep(j)
+                  ptlep(j)=ptlep(j+1)
+                  ptlep(j+1)=tmp
+               endif
+            enddo
+         enddo
+      endif
+
+      
+c missing Et
+      if (ive.ne.0 .and. ivm.ne.0) then
+         etmiss=getptv4_2(p_reco(0,ivm),p_reco(0,ive))
+      elseif (ive.ne.0 .or. ivm.ne.0) then
+         etmiss=getptv4(p_reco(0,ive+ivm))
+      endif
+      
+      if (ih1.ne.0)pth(1)=getptv4(p_reco(0,ih1))
+      if (ih2.ne.0)pth(2)=getptv4(p_reco(0,ih2))
+      if (ih1.ne.0 .and. ih2.ne.0) then
+         pthh=getptv4_2(p_reco(0,ih1),p_reco(0,ih2))
+         Mhh=getinvm4_2(p_reco(0,ih1),p_reco(0,ih2))
+c     order the higgs bosons (if there are 2)
+         if (pth(1).lt.pth(2)) then
+            tmp=pth(1)
+            pth(1)=pth(2)
+            pth(2)=tmp
+         endif
+      endif
+      if (iv1.ne.0) ptv(1)=getptv4(p_reco(0,iv1))
+      if (iv2.ne.0) ptv(2)=getptv4(p_reco(0,iv2))
+      if (iv3.ne.0) ptv(3)=getptv4(p_reco(0,iv3))
+      if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) then
+         Mvvv=getinvm4_3(p_reco(0,iv1),p_reco(0,iv2),p_reco(0,iv3))
+c     order the vector bosons (if there are 3)
+         do i=1,2
+            do j=1,3-i
+               if (ptv(j).lt.ptv(j+1)) then
+                  tmp=ptv(j)
+                  ptv(j)=ptv(j+1)
+                  ptv(j+1)=tmp
+               endif
+            enddo
+         enddo
+      elseif (iv1.ne.0 .and. iv2.ne.0) then
+c     order the vector bosons (if there are 2)
+         if (ptv(1).lt.ptv(2)) then
+            tmp=ptv(1)
+            ptv(1)=ptv(2)
+            ptv(2)=tmp
+         endif
+      endif
+      do i=1,njet
+         ptjet(i)=getptv4(pjet(0,i))
+      enddo
+      if (njet.ge.2) then
+         Mj1j2=getinvm4_2(pjet(0,1),pjet(0,2))
+      endif
+      if (njet.ge.3) then
+         Mj2j3=getinvm4_2(pjet(0,2),pjet(0,3))
+         Mj1j3=getinvm4_2(pjet(0,1),pjet(0,3))
+         Mj1j2j3=getinvm4_3(pjet(0,1),pjet(0,2),pjet(0,3))
+      endif
+c
+      HTparton=0d0
+      HTreco=0d0
+      do i=1,nexternal
+         HTparton=HTparton+getptv4(p(0,i))
+         if (abs(ipdg_reco(i)).gt.5 .and. ipdg_reco(i).ne.21 .and.
+     $        ipdg_reco(i).ne.22 .and. abs(ipdg_reco(i)).ne.12 .and.
+     $        abs(ipdg_reco(i)).ne.14 .and. abs(ipdg_reco(i)).ne.16)
+     $        then
+            HTreco=HTreco+getptv4(p_reco(0,i))
+         endif
+      enddo
+      do i=1,njet
+         HTreco=HTreco+getptv4(pjet(0,i))
+      enddo
+      if (ive.ne.0 .or. ivm.ne.0) HTreco=HTreco+etmiss
+         
+      do i=1,2
+         l=(i-1)*55
+         if (ibody.ne.3 .and.i.eq.2) cycle
+
+c         How to tag orders (QCD+QED*100)
+c
+c         if (i.eq. 3.and.orders_tag_plot.ne.4) cycle
+c         if (i.eq. 4.and.orders_tag_plot.ne.202) cycle
+c         if (i.eq. 5.and.orders_tag_plot.ne.400) cycle
+c         if (i.eq. 6.and.orders_tag_plot.ne.6) cycle
+c         if (i.eq. 7.and.orders_tag_plot.ne.204) cycle
+c         if (i.eq. 8.and.orders_tag_plot.ne.402) cycle
+c         if (i.eq. 9.and.orders_tag_plot.ne.600) cycle
+c         if (i.eq.10.and.orders_tag_plot.ne.8) cycle
+c         if (i.eq.11.and.orders_tag_plot.ne.206) cycle
+c         if (i.eq.12.and.orders_tag_plot.ne.404) cycle
+c         if (i.eq.13.and.orders_tag_plot.ne.602) cycle
+c         if (i.eq.14.and.orders_tag_plot.ne.800) cycle
+
+
+
+c total rate         
+         call HwU_fill(l+ 1,1d0,wgts)
+c transverse momenta         
+         if (il .ge.1) call HwU_fill(l+ 2,l10(ptlep(1)),wgts)
+         if (il .ge.2) call HwU_fill(l+ 3,l10(ptlep(2)),wgts)
+         if (il .ge.3) call HwU_fill(l+ 4,l10(ptlep(3)),wgts)
+         if (il .ge.4) call HwU_fill(l+ 5,l10(ptlep(4)),wgts)
+         if (iem.ne.0) call HwU_fill(l+ 6,l10(ptem),wgts)
+         if (iep.ne.0) call HwU_fill(l+ 7,l10(ptep),wgts)
+         if (imp.ne.0) call HwU_fill(l+ 8,l10(ptmp),wgts)
+         if (imm.ne.0) call HwU_fill(l+ 9,l10(ptmm),wgts)
+         if (iep.ne.0 .and. iem.ne.0) call HwU_fill(l+14,l10(ptepem),wgts)
+         if (imp.ne.0 .and. imm.ne.0) call HwU_fill(l+15,l10(ptmpmm),wgts)
+         if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0) then
+            if (ptepem.gt.ptmpmm) then
+               call HwU_fill(l+16,l10(ptepem),wgts)
+               call HwU_fill(l+17,l10(ptmpmm),wgts)
+            else 
+              call HwU_fill(l+17,l10(ptepem),wgts)
+               call HwU_fill(l+16,l10(ptmpmm),wgts)
+            endif
+         elseif (iep.ne.0 .and. iem.ne.0) then
+            call HwU_fill(l+17,l10(ptepem),wgts)
+         elseif (imp.ne.0 .and. imm.ne.0) then
+            call HwU_fill(l+17,l10(ptmpmm),wgts)
+         endif
+         if (iep.ne.0 .and. ive.ne.0) call HwU_fill(l+18,l10(ptepve),wgts)
+         if (imm.ne.0 .and. ivm.ne.0) call HwU_fill(l+19,l10(ptmmvm),wgts)
+         if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0) then
+            if (ptepve.gt.ptmmvm) then
+               call HwU_fill(l+20,l10(ptepve),wgts)
+               call HwU_fill(l+21,l10(ptmmvm),wgts)
+            else
+               call HwU_fill(l+21,l10(ptepve),wgts)
+               call HwU_fill(l+20,l10(ptmmvm),wgts)
+            endif
+         elseif (iep.ne.0 .and. ive.ne.0) then
+            call HwU_fill(l+20,l10(ptepve),wgts)
+         elseif (imm.ne.0 .and. ivm.ne.0) then
+            call HwU_fill(l+20,l10(ptmmvm),wgts)
+         endif
+         if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0)
+     $        call HwU_fill(l+22,l10(ptepemmpmm),wgts)
+         if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0)
+     $        call HwU_fill(l+23,l10(ptepvemmvm),wgts)
+         if (itop.ne.0) call HwU_fill(l+24,l10(ptt),wgts)
+         if (iatop.ne.0) call HwU_fill(l+25,l10(ptat),wgts)
+         if (itop.ne.0 .and. iatop.ne.0) call HwU_fill(l+26,l10(pttt),wgts)
+         if (ive.ne.0 .or. ivm .ne.0) call HwU_fill(l+27,l10(etmiss),wgts)
+         if (ih1.ne.0) call HwU_fill(l+28,l10(pth(1)),wgts)
+         if (ih2.ne.0) call HwU_fill(l+29,l10(pth(2)),wgts)
+         if (ih1.ne.0 .and. ih2.ne.0) call HwU_fill(l+30,l10(pthh),wgts)
+         if (iv1.ne.0) call HwU_fill(l+31,l10(ptv(1)),wgts)
+         if (iv2.ne.0) call HwU_fill(l+32,l10(ptv(2)),wgts)
+         if (iv3.ne.0) call HwU_fill(l+33,l10(ptv(3)),wgts)
+         if (njet.ge.1) call HwU_fill(l+34,l10(ptjet(1)),wgts)
+         if (njet.ge.2) call HwU_fill(l+35,l10(ptjet(2)),wgts)
+         if (njet.ge.3) call HwU_fill(l+36,l10(ptjet(3)),wgts)
+c invariant masses
+         if (iep.ne.0 .and. iem.ne.0) call HwU_fill(l+37,l10(Mepem),wgts)
+         if (imp.ne.0 .and. imm.ne.0) call HwU_fill(l+38,l10(Mmpmm),wgts)
+         if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0) then
+            if (Mepem.gt.Mmpmm) then
+               call HwU_fill(l+39,l10(Mepem),wgts)
+               call HwU_fill(l+40,l10(Mmpmm),wgts)
+            else
+               call HwU_fill(l+40,l10(Mepem),wgts)
+               call HwU_fill(l+39,l10(Mmpmm),wgts)
+            endif
+         elseif (iep.ne.0 .and. iem.ne.0) then
+            call HwU_fill(l+39,l10(Mepem),wgts)
+         elseif (imp.ne.0 .and. imm.ne.0) then
+            call HwU_fill(l+39,l10(Mmpmm),wgts)
+         endif
+         if (iep.ne.0 .and. ive.ne.0) call HwU_fill(l+41,l10(Mepve),wgts)
+         if (imm.ne.0 .and. ivm.ne.0) call HwU_fill(l+42,l10(Mmmvm),wgts)
+         if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0) then
+            if (Mepve.gt.Mmmvm) then
+               call HwU_fill(l+43,l10(Mepve),wgts)
+               call HwU_fill(l+44,l10(Mmmvm),wgts)
+            else
+               call HwU_fill(l+44,l10(Mepve),wgts)
+               call HwU_fill(l+43,l10(Mmmvm),wgts)
+            endif
+         elseif (iep.ne.0 .and. ive.ne.0) then
+            call HwU_fill(l+43,l10(Mepve),wgts)
+         elseif (imm.ne.0 .and. ivm.ne.0) then
+            call HwU_fill(l+43,l10(Mmmvm),wgts)
+         endif
+         if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0)
+     $        call HwU_fill(l+45,l10(Mepemmpmm),wgts)
+         if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0)
+     $        call HwU_fill(l+46,l10(Mepvemmvm),wgts)
+         if (itop.ne.0 .and. iatop.ne.0) call HwU_fill(l+47,l10(Mtt),wgts)
+         if (ih1.ne.0 .and. ih2.ne.0) call HwU_fill(l+48,l10(Mhh),wgts)
+         if (njet.ge.2) call HwU_fill(l+49,l10(Mj1j2),wgts)
+         if (njet.ge.3) call HwU_fill(l+50,l10(Mj2j3),wgts)
+         if (njet.ge.3) call HwU_fill(l+51,l10(Mj1j3),wgts)
+         if (njet.ge.3) call HwU_fill(l+52,l10(Mj1j2j3),wgts)
+         if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) call HwU_fill(l+53,l10(Mvvv),wgts)
+c HT
+         call HwU_fill(l+54,l10(HTparton),wgts)
+         call HwU_fill(l+55,l10(HTreco),wgts)
+      enddo
+
+ 999  return      
+      end
+
+
+      double precision function l10(var)
+      implicit none
+      double precision var
+      if (var.gt.0) then
+         l10=log10(var)
+      else
+         l10=-99d99
+      endif
+      return
+      end
+      
+      function getrapidity(en,pl)
+      implicit none
+      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
+      parameter (tiny=1.d-8)
+c
+      xplus=en+pl
+      xminus=en-pl
+      if(xplus.gt.tiny.and.xminus.gt.tiny)then
+        if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny )then
+          y=0.5d0*log( xplus/xminus )
+        else
+          y=sign(1.d0,pl)*1.d8
+        endif
+      else
+        y=sign(1.d0,pl)*1.d8
+      endif
+      getrapidity=y
+      return
+      end
+
+
+      function getpseudorap(en,ptx,pty,pl)
+      implicit none
+      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
+      parameter (tiny=1.d-5)
+c
+      pt=sqrt(ptx**2+pty**2)
+      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
+        eta=sign(1.d0,pl)*1.d8
+      else
+        th=atan2(pt,pl)
+        eta=-log(tan(th/2.d0))
+      endif
+      getpseudorap=eta
+      return
+      end
+
+
+      function getinvm(en,ptx,pty,pl)
+      implicit none
+      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
+      parameter (tiny=1.d-5)
+c
+      tmp=en**2-ptx**2-pty**2-pl**2
+      if(tmp.gt.0.d0)then
+        tmp=sqrt(tmp)
+      elseif(tmp.gt.-tiny)then
+        tmp=0.d0
+      else
+        write(*,*)'Attempt to compute a negative mass'
+        stop
+      endif
+      getinvm=tmp
+      return
+      end
+
+
+      function getdelphi(ptx1,pty1,ptx2,pty2)
+      implicit none
+      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
+      parameter (tiny=1.d-5)
+c
+      pt1=sqrt(ptx1**2+pty1**2)
+      pt2=sqrt(ptx2**2+pty2**2)
+      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
+        tmp=ptx1*ptx2+pty1*pty2
+        tmp=tmp/(pt1*pt2)
+        if(abs(tmp).gt.1.d0+tiny)then
+          write(*,*)'Cosine larger than 1'
+          stop
+        elseif(abs(tmp).ge.1.d0)then
+          tmp=sign(1.d0,tmp)
+        endif
+        tmp=acos(tmp)
+      else
+        tmp=1.d8
+      endif
+      getdelphi=tmp
+      return
+      end
+
+
+      function getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
+      implicit none
+      real*8 getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
+     # getpseudorap,getdelphi
+c
+      deta=getpseudorap(en1,ptx1,pty1,pl1)-
+     #     getpseudorap(en2,ptx2,pty2,pl2)
+      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
+      getdr=sqrt(dphi**2+deta**2)
+      return
+      end
+
+
+      function getdry(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
+      implicit none
+      real*8 getdry,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
+     # getrapidity,getdelphi
+c
+      deta=getrapidity(en1,pl1)-
+     #     getrapidity(en2,pl2)
+      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
+      getdry=sqrt(dphi**2+deta**2)
+      return
+      end
+
+
+      function getptv(p)
+      implicit none
+      real*8 getptv,p(5)
+c
+      getptv=sqrt(p(1)**2+p(2)**2)
+      return
+      end
+
+
+      function getpseudorapv(p)
+      implicit none
+      real*8 getpseudorapv,p(5)
+      real*8 getpseudorap
+c
+      getpseudorapv=getpseudorap(p(4),p(1),p(2),p(3))
+      return
+      end
+
+
+      function getrapidityv(p)
+      implicit none
+      real*8 getrapidityv,p(5)
+      real*8 getrapidity
+c
+      getrapidityv=getrapidity(p(4),p(3))
+      return
+      end
+
+
+      function getdrv(p1,p2)
+      implicit none
+      real*8 getdrv,p1(5),p2(5)
+      real*8 getdr
+c
+      getdrv=getdr(p1(4),p1(1),p1(2),p1(3),
+     #             p2(4),p2(1),p2(2),p2(3))
+      return
+      end
+
+
+      function getinvmv(p)
+      implicit none
+      real*8 getinvmv,p(5)
+      real*8 getinvm
+c
+      getinvmv=getinvm(p(4),p(1),p(2),p(3))
+      return
+      end
+
+
+      function getdelphiv(p1,p2)
+      implicit none
+      real*8 getdelphiv,p1(5),p2(5)
+      real*8 getdelphi
+c
+      getdelphiv=getdelphi(p1(1),p1(2),
+     #                     p2(1),p2(2))
+      return
+      end
+
+
+      function getptv4(p)
+      implicit none
+      real*8 getptv4,p(0:3)
+c
+      getptv4=sqrt(max(p(1)**2+p(2)**2,0d0))
+      return
+      end
+
+      function getptv4_2(p1,p2)
+      implicit none
+      real*8 getptv4_2,p1(0:3),p2(0:3),psum(0:3)
+      integer i
+      do i=0,3
+         psum(i)=p1(i)+p2(i)
+      enddo
+      getptv4_2=sqrt(max(psum(1)**2+psum(2)**2,0d0))
+      return
+      end
+
+      function getptv4_4(p1,p2,p3,p4)
+      implicit none
+      real*8 getptv4_4,p1(0:3),p2(0:3),p3(0:3),p4(0:3),psum(0:3)
+      integer i
+      do i=0,3
+         psum(i)=p1(i)+p2(i)+p3(i)+p4(i)
+      enddo
+      getptv4_4=sqrt(max(psum(1)**2+psum(2)**2,0d0))
+      return
+      end
+
+
+      function getpseudorapv4(p)
+      implicit none
+      real*8 getpseudorapv4,p(0:3)
+      real*8 getpseudorap
+c
+      getpseudorapv4=getpseudorap(p(0),p(1),p(2),p(3))
+      return
+      end
+
+
+      function getrapidityv4(p)
+      implicit none
+      real*8 getrapidityv4,p(0:3)
+      real*8 getrapidity
+c
+      getrapidityv4=getrapidity(p(0),p(3))
+      return
+      end
+
+
+      function getdrv4(p1,p2)
+      implicit none
+      real*8 getdrv4,p1(0:3),p2(0:3)
+      real*8 getdr
+c
+      getdrv4=getdr(p1(0),p1(1),p1(2),p1(3),
+     #              p2(0),p2(1),p2(2),p2(3))
+      return
+      end
+
+
+      function getinvm4(p)
+      implicit none
+      real*8 getinvm4,p(0:3)
+      real*8 getinvm
+c
+      getinvm4=getinvm(p(0),p(1),p(2),p(3))
+      return
+      end
+
+      function getinvm4_2(p1,p2)
+      implicit none
+      real*8 getinvm4_2,p1(0:3),p2(0:3),p(0:3)
+      real*8 getinvm
+      integer i
+      do i=0,3
+         p(i)=p1(i)+p2(i)
+      enddo
+      getinvm4_2=getinvm(p(0),p(1),p(2),p(3))
+      return
+      end
+
+      function getinvm4_3(p1,p2,p3)
+      implicit none
+      real*8 getinvm4_3,p1(0:3),p2(0:3),p3(0:3),p(0:3)
+      real*8 getinvm
+      integer i
+      do i=0,3
+         p(i)=p1(i)+p2(i)+p3(i)
+      enddo
+      getinvm4_3=getinvm(p(0),p(1),p(2),p(3))
+      return
+      end
+
+      function getinvm4_4(p1,p2,p3,p4)
+      implicit none
+      real*8 getinvm4_4,p1(0:3),p2(0:3),p3(0:3),p4(0:3),p(0:3)
+      real*8 getinvm
+      integer i
+      do i=0,3
+         p(i)=p1(i)+p2(i)+p3(i)+p4(i)
+      enddo
+      getinvm4_4=getinvm(p(0),p(1),p(2),p(3))
+      return
+      end
+
+
+      function getdelphiv4(p1,p2)
+      implicit none
+      real*8 getdelphiv4,p1(0:3),p2(0:3)
+      real*8 getdelphi
+c
+      getdelphiv4=getdelphi(p1(1),p1(2),
+     #                      p2(1),p2(2))
+      return
+      end
+
+
+      function getcosv4(q1,q2)
+      implicit none
+      real*8 getcosv4,q1(0:3),q2(0:3)
+      real*8 xnorm1,xnorm2,tmp
+c
+      if(q1(0).lt.0.d0.or.q2(0).lt.0.d0)then
+        getcosv4=-1.d10
+        return
+      endif
+      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
+      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
+      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
+        tmp=-1.d10
+      else
+        tmp=q1(1)*q2(1)+q1(2)*q2(2)+q1(3)*q2(3)
+        tmp=tmp/(xnorm1*xnorm2)
+        if(abs(tmp).gt.1.d0.and.abs(tmp).le.1.001d0)then
+          tmp=sign(1.d0,tmp)
+        elseif(abs(tmp).gt.1.001d0)then
+          write(*,*)'Error in getcosv4',tmp
+          stop
+        endif
+      endif
+      getcosv4=tmp
+      return
+      end
+
+
+
+      function getmod(p)
+      implicit none
+      double precision p(0:3),getmod
+
+      getmod=sqrt(p(1)**2+p(2)**2+p(3)**2)
+
+      return
+      end
+
+
+
+      subroutine getperpenv4(q1,q2,qperp)
+c Normal to the plane defined by \vec{q1},\vec{q2}
+      implicit none
+      real*8 q1(0:3),q2(0:3),qperp(0:3)
+      real*8 xnorm1,xnorm2
+      integer i
+c
+      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
+      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
+      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
+        do i=1,4
+          qperp(i)=-1.d10
+        enddo
+      else
+        qperp(1)=q1(2)*q2(3)-q1(3)*q2(2)
+        qperp(2)=q1(3)*q2(1)-q1(1)*q2(3)
+        qperp(3)=q1(1)*q2(2)-q1(2)*q2(1)
+        do i=1,3
+          qperp(i)=qperp(i)/(xnorm1*xnorm2)
+        enddo
+        qperp(0)=1.d0
+      endif
+      return
+      end
+
+
+
+
+      subroutine getwedge(p1,p2,pout)
+      implicit none
+      real*8 p1(0:3),p2(0:3),pout(0:3)
+
+      pout(1)=p1(2)*p2(3)-p1(3)*p2(2)
+      pout(2)=p1(3)*p2(1)-p1(1)*p2(3)
+      pout(3)=p1(1)*p2(2)-p1(2)*p2(1)
+      pout(0)=0d0
+
+      return
+      end
+

=== modified file 'Template/NLO/MCatNLO/Makefile_MadFKS'
--- Template/NLO/MCatNLO/Makefile_MadFKS	2016-02-24 00:19:50 +0000
+++ Template/NLO/MCatNLO/Makefile_MadFKS	2018-12-20 08:18:55 +0000
@@ -1,7 +1,7 @@
 -include ../../Source/make_opts
 
 DEBUG=
-FF=$(FC) $(FFLAGS) $(DEBUG)
+FF=$(FC) $(FFLAGS) $(DEBUG) -std=legacy
 CPP=$(CXX) $(CXXFLAGS) $(DEBUG) $(INCLOPTION)
 CC=$(CXX) $(CFLAGS) $(DEBUG) $(INCLOPTION)
 

=== modified file 'Template/NLO/SubProcesses/MC_integer.f'
--- Template/NLO/SubProcesses/MC_integer.f	2015-02-09 18:56:09 +0000
+++ Template/NLO/SubProcesses/MC_integer.f	2018-12-20 08:18:55 +0000
@@ -68,6 +68,7 @@
       common/to_readgrid/flat_grid                !Tells if grid read from file
       if (this_dim.lt.1.or.this_dim.gt.maxdim) then
          write (*,*) 'Increase maxdim in MC_integer.f',maxdim,this_dim
+         stop 1
       endif
 c Set the number of intervales for all the dimensions to zero the very
 c first time this subroutine is called
@@ -283,12 +284,8 @@
       enddo
 c
 c Reset the accumulated results because we start new iteration.
-      do this_dim=1,maxdim
-         do i=0,nintervals(this_dim)
-            acc(i,this_dim)=0d0
-            ncall(i,this_dim)=0
-         enddo
-      enddo
+      call empty_MC_integer
+c
       return
  999  write (*,*) 'Cannot open "grid.MC_integer" file'
       stop

=== modified file 'Template/NLO/SubProcesses/analysis_lhe.f'
--- Template/NLO/SubProcesses/analysis_lhe.f	2017-06-19 16:05:33 +0000
+++ Template/NLO/SubProcesses/analysis_lhe.f	2018-12-20 08:18:55 +0000
@@ -7,9 +7,10 @@
       implicit none
       integer nwgt
       character*(*) weights_info(*)
+      logical passed_unwgt
       integer nwgts,nevents
       double precision sum_of_wgts
-      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
+      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
       logical lopen
       nwgts=nwgt
       sum_of_wgts=0d0
@@ -23,9 +24,10 @@
       subroutine analysis_end(xnorm)
       implicit none
       double precision xnorm
+      logical passed_unwgt
       integer nwgts,nevents
       double precision sum_of_wgts
-      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
+      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
       logical lopen
       integer                                   npoints
       double precision            cross_section
@@ -55,9 +57,10 @@
       double precision p(0:4,nexternal)
       double precision wgts(*)
       integer ibody
+      logical passed_unwgt
       integer nwgts,nevents
       double precision sum_of_wgts
-      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
+      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
 c
       integer i,j,npart
       double precision ran2
@@ -95,23 +98,33 @@
       integer                                   npoints
       double precision            cross_section
       common /for_FixedOrder_lhe/ cross_section,npoints
+      double precision,allocatable :: wwgts(:)
 c --- do the partial unweighting (don't do it when to close to singular
 c --- region)
+      if (.not. allocated(wwgts)) then
+         allocate(wwgts(nwgts))
+      endif
+      passed_unwgt=.true.
       twgt=abs(cross_section)*FO_LHE_weight_ratio !/npoints
       if(abs(wgts(1)).lt.abs(twgt) .and. xi_i_fks_ev.gt.1d-3 .and.
      $     1d0-y_ij_fks_ev.gt.1d-2)then
          R = ran2()*abs(twgt)
          if (R.gt.abs(wgts(1)))then
+            passed_unwgt=.false.
             return
          else
             do i=2,nwgts
-               wgts(i) = wgts(i)*abs(twgt/wgts(1))
+               wwgts(i) = wgts(i)*abs(twgt/wgts(1))
             enddo
-            wgts(1) = sign(twgt,wgts(1))
+            wwgts(1) = sign(twgt,wgts(1))
          endif
+      else
+         do i=1,nwgts
+            wwgts(i)=wgts(i)
+         enddo
       endif
 c --- accumulate the total weights of all the events in the event file
-      sum_of_wgts=sum_of_wgts+wgts(1)
+      sum_of_wgts=sum_of_wgts+wwgts(1)
 
 c --- fill the multi-weight common blocks with the scale and PDF
 c --- variation weights
@@ -122,12 +135,12 @@
                do ii=1,nint(scalevarF(0))
                   do jj=1,nint(scalevarR(0))
                      i_wgt=i_wgt+1
-                     wgtxsecmu(jj,ii,kk)= wgts(i_wgt)
+                     wgtxsecmu(jj,ii,kk)= wwgts(i_wgt)
                   enddo
                enddo
             else
                i_wgt=i_wgt+1
-               wgtxsecmu(1,1,kk)= wgts(i_wgt)
+               wgtxsecmu(1,1,kk)= wwgts(i_wgt)
             endif
          enddo
       endif
@@ -136,11 +149,11 @@
             if (lpdfvar(nn)) then
                do n=0,nmemPDF(nn)
                   i_wgt=i_wgt+1
-                  wgtxsecPDF(n,nn) = wgts(i_wgt)
+                  wgtxsecPDF(n,nn) = wwgts(i_wgt)
                enddo
             else
                i_wgt=i_wgt+1
-               wgtxsecPDF(0,nn) = wgts(i_wgt)
+               wgtxsecPDF(0,nn) = wwgts(i_wgt)
             endif
          enddo
       endif
@@ -187,7 +200,7 @@
 
 c --- write the event
       call write_lhef_event(41,
-     &     npart,IDPRUP,wgts(1),0d0,0d0,0d0,
+     &     npart,IDPRUP,wwgts(1),0d0,0d0,0d0,
      &     IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
  201  format(a9,1x,i1,4(1x,i2),2(1x,d14.8),2x,i2,2(1x,i2),5(1x,d14.8))
@@ -196,12 +209,15 @@
 c This we can use for the event grouping!      
       subroutine HwU_add_points
       implicit none
+      logical passed_unwgt
       integer nwgts,nevents
       double precision sum_of_wgts
-      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
-      nevents=nevents+1
-      write (41,'(a)') '</eventgroup>'
-      write (41,'(a)') '<eventgroup>'
+      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
+      if (passed_unwgt) then
+         nevents=nevents+1
+         write (41,'(a)') '</eventgroup>'
+         write (41,'(a)') '<eventgroup>'
+      endif
       end
 
 c Dummy routines

=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2017-07-31 08:48:20 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2018-12-20 08:18:55 +0000
@@ -18,7 +18,7 @@
       include 'genps.inc'
       include 'nexternal.inc'
       include 'nFKSconfigs.inc'
-      double precision p(0:3, nexternal), prambo(0:3, nexternal)
+      double precision p(0:3, nexternal), prambo(0:3,100)
       double precision p_born(0:3,nexternal-1)
       common/pborn/p_born
       double precision pswgt
@@ -48,7 +48,7 @@
       include 'coupl.inc'
       include 'q_es.inc'
       integer nsqso,MLResArrayDim
-      double precision pmass(nexternal), pmass_rambo(nexternal)
+      double precision pmass(nexternal), pmass_rambo(100)
       integer nfail
       logical first_time
       data first_time/.TRUE./

=== modified file 'Template/NLO/SubProcesses/combine_plots_FO.sh'
--- Template/NLO/SubProcesses/combine_plots_FO.sh	2013-05-25 12:52:41 +0000
+++ Template/NLO/SubProcesses/combine_plots_FO.sh	2018-12-20 08:18:55 +0000
@@ -1,6 +1,4 @@
 #!/bin/bash
-TD=td_mac_intel
-PSOPEN=open
 
 # $? is the value of last executed command. A call to this function
 # after a failure will cause the program to quit the script
@@ -15,126 +13,74 @@
 fi
 }
 
-thisdir=`pwd`
-if [ -f $thisdir/MADatNLO.top ]
-then
-rm -f $thisdir/MADatNLO*.top
-fi
-make read40
-EXENAME=$thisdir/read40
+function combine {
+counter=0
+i=1
 echo -n "" > dir
-counterp=0
-for p in P* ; do
-    cd $thisdir
-    cd $p
-    rm -f MADatNLO*
-    rm -f *read40*.out
-    echo -n "" > dir
-    counter=0
-    i=1
-    for g in $* ; do
-        if [ ! -f $g/MADatNLO*.top ]
-        then
-          echo "No topdrawer file in $p/$g"
-	else
-	counter=$(($counter + 1))
-	echo $g"/MADatNLO.top" >> dir
-	if [[ $(($counter % 40)) == 0 ]]; then
-            echo -n "" > dir2
-	    echo "40" > dir2
-	    cat dir >> dir2
-	    $EXENAME < dir2
-            teststatus Failure in step 1
-	    rm -f dir
-	    rm -f dir2
-	    mv fort.88 MADatNLO$i\.top 
-	    mv read40.out "S1read40_"$i\.out
-	    i=$(($i + 1))
-	    echo -n "" > dir
-        fi
-	fi
-    done
-    if [[ $(($counter % 40)) != 0 ]] ; then
-	echo $(($counter % 40)) > dir2
+for plot_file in "${@:2}" ; do
+    counter=$(($counter + 1))
+    echo $plot_file >> dir
+    if [[ $(($counter % 40)) == 0 ]]; then
+        echo -n "" > dir2
+	echo "40" > dir2
 	cat dir >> dir2
-	while [[ $(($counter % 40)) -ne 0 ]]; do
-	    counter=$[$counter-1]
-	done
-	$EXENAME < dir2
-	teststatus Failure in step 2
+	./read40 < dir2
+        teststatus Failure in step 1
 	rm -f dir
 	rm -f dir2
-	mv fort.88 MADatNLO$i\.top 
-	mv read40.out "S2read40_"$i\.out
-    fi
-	
-    counter=0
-    for m in MADatNLO*.top ; do
-	counter=$[$counter+1]
-	echo $m >> dir
-    done
-    if [[ $counter -gt 1 ]] ; then
-	echo $counter > dir2
-	cat dir >> dir2
-	$EXENAME < dir2
-	teststatus Failure in step 3
-	rm -f dir2
-	mv fort.88 MADatNLO.top 
-	mv read40.out S3read40.out
-    else
-	mv MADatNLO1.top MADatNLO.top
-    fi
-    rm -f dir
-    cd $thisdir
-    counterp=$(($counterp + 1))
-    echo $p"/MADatNLO.top" >> dir
-    echo $counterp $p "done"
-done
-
-cd $thisdir
-make split40
-./split40
-
-i=0
-#pdir* are created by split
-for p in pdir* ; do
-    i=$(($i+1))
-    $EXENAME < $p    
-    teststatus Failure in step 4
-    mv fort.88 MADatNLO$i\.top
-    mv read40.out "S4read40_"$i\.out
-    rm $p
-    echo $i $p "done"
-done
-
-if [[ $i -gt 1 ]] ; then
-    echo $i > dir2
-    for p in MADatNLO*.top ; do
-	echo $p >>dir2
-    done
-    $EXENAME < dir2
-    teststatus Failure in step 5
+	mv fort.88 $1_$i\.top 
+	mv read40.out "S1read40_"$1_$i\.out
+	i=$(($i + 1))
+	echo -n "" > dir
+    fi
+done
+if [[ $(($counter % 40)) != 0 ]] ; then
+    echo $(($counter % 40)) > dir2
+    cat dir >> dir2
+    ./read40 < dir2
+    teststatus Failure in step 2
     rm -f dir
     rm -f dir2
-    mv fort.88 MADatNLO.top 
-    mv read40.out S5read40.out
+    mv fort.88 $1_$i\.top 
+    mv read40.out "S2read40_"$1_$i\.out
+fi
+}
+
+make read40
+rm -f MADatNLO*.top
+
+if [ "$#" -gt 64000 ]; then
+    combine "MADatNLO_A" "$@"
+    combine "MADatNLO_B" MADatNLO_A*.top
+    combine "MADatNLO_C" MADatNLO_B*.top
+    combine "MADatNLO_D" MADatNLO_C*.top
+    mv MADatNLO_D_1.top MADatNLO.top
+elif [ "$#" -gt 1600 ]; then
+    combine "MADatNLO_A" "$@"
+    combine "MADatNLO_B" MADatNLO_A*.top
+    combine "MADatNLO_C" MADatNLO_B*.top
+    mv MADatNLO_C_1.top MADatNLO.top
+elif [ "$#" -gt 40 ]; then
+    combine "MADatNLO_A" "$@"
+    combine "MADatNLO_B" MADatNLO_A*.top
+    mv MADatNLO_B_1.top MADatNLO.top
+elif [ "$#" -gt 1 ]; then
+    combine "MADatNLO_A" "$@"
+    mv MADatNLO_A_1.top MADatNLO.top
 else
-    mv MADatNLO1.top MADatNLO.top
+    cp $1 MADatNLO.top
 fi
 
-if [ -f $thisdir/read40.errors ]
+
+if [ -f read40.errors ]
 then
-rm -f $thisdir/read40.errors
+rm -f read40.errors
 fi
-find . -name S"*"read40"*".out -exec fgrep -il "READ40 ERROR READ40 ERROR" {} \; > $thisdir/read40.errors
-
-if [ ! -s $thisdir/read40.errors ] ; then
+
+find . -name S"*"read40"*".out -exec fgrep -il "READ40 ERROR READ40 ERROR" {} \; > read40.errors
+
+if [ ! -s read40.errors ] ; then
     echo 'no errors found'
 else
     echo 'ERRORS! see read40.errors for details'
 fi
-
-#./NLO_Born3.py
-#$TD MADatNLO_combined.top "POSTSCR,ORIENT=3"
-#$PSOPEN MADatNLO_combined.ps
-##td2pdf MADatNLO_con

=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2018-04-16 14:08:47 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2018-12-20 08:18:55 +0000
@@ -48,7 +48,7 @@
       include "mint.inc"
       integer nhits_in_grids(maxchannels)
       real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
-     $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
+     $     ,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
@@ -402,9 +402,10 @@
       double precision xx(ndimmax),vegas_wgt,f(nintegrals),jac,p(0:3
      $     ,nexternal),rwgt,vol,sig,x(99),MC_int_wgt
       integer ifl,nFKS_born,nFKS_picked,iFKS,nFKS_min,iamp
-     $     ,nFKS_max,izero,ione,itwo,mohdr
+     $     ,nFKS_max,izero,ione,itwo,mohdr,i,iran_picked
       parameter (izero=0,ione=1,itwo=2,mohdr=-100)
-      logical passcuts,passcuts_nbody,passcuts_n1body,sum
+      logical passcuts,passcuts_nbody,passcuts_n1body,sum,firsttime
+      data firsttime/.true./
       external passcuts
       integer             ini_fin_fks(maxchannels)
       common/fks_channels/ini_fin_fks
@@ -431,6 +432,16 @@
       common /for_applgrid/ iappl
       double precision       wgt_ME_born,wgt_ME_real
       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 (firsttime) then
+         firsttime=.false.
+         call setup_ini_fin_fks_map(ini_fin_fks_map)
+         write (*,*) 'initial-final FKS maps:'
+         write (*,*) 0 ,':',ini_fin_fks_map(0,:)
+         write (*,*) 1 ,':',ini_fin_fks_map(1,:)
+         write (*,*) 2 ,':',ini_fin_fks_map(2,:)
+      endif
       if (ifl.ne.0) then
          write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
          stop 1
@@ -453,20 +464,9 @@
       if (ickkw.eq.-1) H1_factor_virt=0d0
       if (ickkw.eq.3) call set_FxFx_scale(0,p)
       call update_vegas_x(xx,x)
-      call get_MC_integer(max(ini_fin_fks(ichan),1),fks_configs
-     $     ,nFKS_picked,vol)
-      if (ini_fin_fks(ichan).eq.1) then
-         do while (fks_j_d(nFKS_picked).le.nincoming) 
-            call get_MC_integer(ini_fin_fks(ichan),fks_configs
-     $           ,nFKS_picked,vol)
-         enddo
-      elseif (ini_fin_fks(ichan).eq.2) then
-         do while (fks_j_d(nFKS_picked).gt.nincoming) 
-            call get_MC_integer(ini_fin_fks(ichan),fks_configs
-     $           ,nFKS_picked,vol)
-         enddo
-      endif
-
+      call get_MC_integer(max(ini_fin_fks(ichan),1)
+     $     ,ini_fin_fks_map(ini_fin_fks(ichan),0),iran_picked,vol)
+      nFKS_picked=ini_fin_fks_map(ini_fin_fks(ichan),iran_picked)
       
 c The nbody contributions
       if (abrv.eq.'real') goto 11
@@ -503,14 +503,15 @@
       nbody=.false.
       if (sum) then
          nFKS_min=1
-         nFKS_max=fks_configs
+         nFKS_max=ini_fin_fks_map(ini_fin_fks(ichan),0)
          MC_int_wgt=1d0
       else
-         nFKS_min=nFKS_picked
-         nFKS_max=nFKS_picked
+         nFKS_min=iran_picked
+         nFKS_max=iran_picked
          MC_int_wgt=1d0/vol
       endif
-      do iFKS=nFKS_min,nFKS_max
+      do i=nFKS_min,nFKS_max
+         iFKS=ini_fin_fks_map(ini_fin_fks(ichan),i)
          calculatedBorn=.false. 
          ! MZ this is a temporary fix for processes without
          ! soft singularities associated to the initial state
@@ -566,11 +567,11 @@
 c Importance sampling for FKS configurations
       if (sum) then
          call get_wgt_nbody(sig)
-         call fill_MC_integer(max(ini_fin_fks(ichan),1),nFKS_picked
+         call fill_MC_integer(max(ini_fin_fks(ichan),1),iran_picked
      $        ,abs(sig))
       else
          call get_wgt_no_nbody(sig)
-         call fill_MC_integer(max(ini_fin_fks(ichan),1),nFKS_picked
+         call fill_MC_integer(max(ini_fin_fks(ichan),1),iran_picked
      $        ,abs(sig)*vol)
       endif
 
@@ -595,6 +596,35 @@
       return
       end
       
+      subroutine setup_ini_fin_FKS_map(ini_fin_FKS_map)
+      implicit none
+      include 'nexternal.inc'
+      include 'nFKSconfigs.inc'
+      include 'fks_info.inc'
+      integer ini_fin_FKS_map(0:2,0:fks_configs),iFKS
+      ini_fin_FKS_map(0,0)=0
+      ini_fin_FKS_map(1,0)=0
+      ini_fin_FKS_map(2,0)=0
+      do iFKS=1,fks_configs
+         ini_fin_FKS_map(0,0)=ini_fin_FKS_map(0,0)+1
+         ini_fin_FKS_map(0,ini_fin_FKS_map(0,0))=iFKS
+         if (fks_j_d(iFKS).le.nincoming .and.
+     $       fks_j_d(iFKS).gt.0) then
+            ini_fin_FKS_map(2,0)=ini_fin_FKS_map(2,0)+1
+            ini_fin_FKS_map(2,ini_fin_FKS_map(2,0))=iFKS
+         elseif (fks_j_d(iFKS).gt.nincoming .and.
+     $           fks_j_d(iFKS).le.nexternal) then
+            ini_fin_FKS_map(1,0)=ini_fin_FKS_map(1,0)+1
+            ini_fin_FKS_map(1,ini_fin_FKS_map(1,0))=iFKS
+         else
+            write (*,*) 'ERROR in setup_ini_fin_FKS_map',fks_j_d(iFKS)
+     $           ,nincoming,iFKS
+            stop 1
+         endif
+      enddo
+      return
+      end
+
       subroutine get_born_nFKSprocess(nFKS_in,nFKS_out)
       implicit none
       include 'nexternal.inc'

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2018-04-16 14:08:47 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2018-12-20 08:18:55 +0000
@@ -59,7 +59,7 @@
       common /event_normalisation/event_norm
 c For MINT:
       real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
-     $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
+     $     ,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)
@@ -842,7 +842,7 @@
          firsttime=.false.
 c Determines the proc_map that sets which FKS configuration can be
 c summed explicitly and which by MC-ing.
-         call setup_proc_map(sum,proc_map)
+         call setup_proc_map(sum,proc_map,ini_fin_fks)
 c For the S-events, we can combine processes when they give identical
 c processes at the Born. Make sure we check that we get indeed identical
 c IRPOC's
@@ -869,16 +869,6 @@
          if (ickkw.eq.3) call set_FxFx_scale(0,p)
          call update_vegas_x(xx,x)
          call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
-         if (ini_fin_fks.eq.1) then
-            do while (fks_j_d(proc_map(proc_map(0,1),1)).le.nincoming) 
-               call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
-            enddo
-         elseif (ini_fin_fks.eq.2) then
-            do while (fks_j_d(proc_map(proc_map(0,1),1)).gt.nincoming) 
-               call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
-            enddo
-         endif
-
 c The nbody contributions
          if (abrv.eq.'real') goto 11
          nbody=.true.
@@ -1043,7 +1033,7 @@
       end
 
 
-      subroutine setup_proc_map(sum,proc_map)
+      subroutine setup_proc_map(sum,proc_map,ini_fin_fks)
 c Determines the proc_map that sets which FKS configuration can be
 c summed explicitly and which by MC-ing.
       implicit none
@@ -1056,7 +1046,7 @@
       logical found_ini1,found_ini2,found_fnl
       integer proc_map(0:fks_configs,0:fks_configs)
      $     ,j_fks_proc(fks_configs),i_fks_pdg_proc(fks_configs)
-     $     ,j_fks_pdg_proc(fks_configs),i,sum,j
+     $     ,j_fks_pdg_proc(fks_configs),i,sum,j,ini_fin_fks
       integer              nFKSprocess
       common/c_nFKSprocess/nFKSprocess
       INTEGER              IPROC
@@ -1097,6 +1087,8 @@
 c them in the process map
          do nFKSprocess=1,fks_configs
             call fks_inc_chooser()
+            if (ini_fin_fks.eq.1 .and. j_fks.le.nincoming) cycle
+            if (ini_fin_fks.eq.2 .and. j_fks.gt.nincoming) cycle
             if (need_color_links.or.need_charge_links) then
                proc_map(0,0)=proc_map(0,0)+1
                proc_map(proc_map(0,0),0)=proc_map(proc_map(0,0),0)+1
@@ -1156,6 +1148,8 @@
 c gluons splitting
          do nFKSprocess=1,fks_configs
             call fks_inc_chooser()
+            if (ini_fin_fks.eq.1 .and. j_fks.le.nincoming) cycle
+            if (ini_fin_fks.eq.2 .and. j_fks.gt.nincoming) cycle
             if (.not.(need_color_links.or.need_charge_links)) then
                if (j_fks.eq.1 .and. found_ini1) then
                   do i=1,proc_map(0,0)

=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2018-04-16 14:08:47 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2018-12-20 08:18:55 +0000
@@ -1115,8 +1115,6 @@
       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
-      double precision  xi_i_hat_ev,xi_i_hat_cnt(-2:2)
-      common /cxi_i_hat/xi_i_hat_ev,xi_i_hat_cnt
       integer            i_fks,j_fks
       common/fks_indices/i_fks,j_fks
       double precision   xi_i_fks_cnt(-2:2)
@@ -1249,7 +1247,7 @@
       endif
 
 c f_* multiplication factors for real-emission, soft counter, ... etc.       
-      prefact=1d0/xi_i_hat_ev/(1-y_ij_fks_ev)
+      prefact=xinorm_ev/xi_i_fks_ev/(1-y_ij_fks_ev)
       f_r=prefact*jac_ev*enhance_real*fkssymmetryfactor*vegas_wgt
       f_MC_S=f_r
       f_MC_H=f_r
@@ -1736,11 +1734,16 @@
 c coefficients for PDF and scale computations. 
       use weight_lines
       implicit none
-      integer i,j
+      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,born_wgt_mint
+      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
@@ -1749,6 +1752,7 @@
          enddo
          return
       endif
+      virt_found=.false.
 c loop over all contributions
       do i=1,icontr
          if (itype(i).eq.1) then
@@ -1771,9 +1775,15 @@
          do j=1,3
             wgt(j,i)=wgt(j,i)*bias_wgt(i)
          enddo
-         if (itype(i).eq.14) then
-            virt_wgt_mint=virt_wgt_mint*bias_wgt(i)
-            born_wgt_mint=born_wgt_mint*bias_wgt(i)
+         if (itype(i).eq.14 .and. .not.virt_found) then
+            virt_found=.true.
+            virt_wgt_mint(0)=virt_wgt_mint(0)*bias_wgt(i)
+            born_wgt_mint(0)=born_wgt_mint(0)*bias_wgt(i)
+            do iamp=1,amp_split_size
+               call amp_split_pos_to_orders(iamp, orders)
+               virt_wgt_mint(iamp)=virt_wgt_mint(iamp)*bias_wgt(i)
+               born_wgt_mint(iamp)=born_wgt_mint(iamp)*bias_wgt(i)
+            enddo
          endif
       enddo
       return
@@ -3805,7 +3815,7 @@
       if (iextra_cnt.gt.0)
      1    call extra_cnt(p_born, iextra_cnt, ans_extra_cnt)
       call AP_reduced(j_type,i_type,ch_j,ch_i,t,z,ap)
-      call Qterms_reduced_timelike(j_type,i_type,ch_i,ch_j,t,z,Q)
+      call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,t,z,Q)
       wgt=0d0
       do iord = 1, nsplitorders
          if (.not.split_type(iord) .or. (iord.ne.qed_pos .and.

=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
--- Template/NLO/SubProcesses/genps_fks.f	2018-01-29 05:41:33 +0000
+++ Template/NLO/SubProcesses/genps_fks.f	2018-12-20 08:18:55 +0000
@@ -2557,7 +2557,7 @@
       integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
       double precision cBW_mass(-1:1,-nexternal:-1),
      &     cBW_width(-1:1,-nexternal:-1)
-      double precision s_mass(-nexternal:-1)
+      double precision s_mass(-nexternal:nexternal)
       common/to_phase_space_s_channel/s_mass
 c Common block with granny information
       logical granny_is_res

=== modified file 'Template/NLO/SubProcesses/makefile'
--- Template/NLO/SubProcesses/makefile	2017-06-15 16:22:23 +0000
+++ Template/NLO/SubProcesses/makefile	2018-12-20 08:18:55 +0000
@@ -7,9 +7,6 @@
 # Files for the read40 combiner of top drawer files
 READ40=read40.o
 
-# Files for the split40 splitter of top drawer files
-SPLIT40=split40.o
-
 # Files to collect all the events in the separate integration channels into a single file
 COLLECT_EVENTS=collect_events.o handling_lhe_events.o fill_MC_mshell.o
 
@@ -21,9 +18,6 @@
 read40: $(READ40)
 	$(FC) $(LDFLAGS) -o read40 $(READ40)
 
-split40: $(SPLIT40)
-	$(FC) $(LDFLAGS) -o split40 $(SPLIT40)
-
 collect_events: $(COLLECT_EVENTS) $(LIBS)
 	$(FC) -o collect_events $(COLLECT_EVENTS) $(LINKLIBS) $(LDFLAGS)
 	rm handling_lhe_events.o

=== modified file 'Template/NLO/SubProcesses/mint-integrator2.f'
--- Template/NLO/SubProcesses/mint-integrator2.f	2017-07-31 08:48:20 +0000
+++ Template/NLO/SubProcesses/mint-integrator2.f	2018-12-20 08:18:55 +0000
@@ -77,7 +77,7 @@
      $     ,ans(nintegrals,0:maxchannels),unc(nintegrals,0:maxchannels)
      $     ,ans3(nintegrals,3),unc3(nintegrals,3),ans_l3(nintegrals)
      $     ,unc_l3(nintegrals),chi2_l3(nintegrals),vol_chan,error_virt
-      real * 8 xint_virt(maxchannels),ymax_virt(maxchannels)
+      real * 8 xint_virt(maxchannels),ymax_virt(0:maxchannels)
      $     ,ans_chan(0:maxchannels)
       real * 8 x(ndimmax),vol
       real * 8 xacc(0:nintervals,ndimmax,maxchannels)
@@ -216,7 +216,7 @@
             ans_chan(0)=ans_chan(0)+ans(1,kchan)
          enddo
       endif
-      cross_section=ans_chan(0)
+      cross_section=ans_chan(0) * wgt_mult
       nit=0
       nit_included=0
       do i=1,nintegrals
@@ -1031,8 +1031,8 @@
       integer ndim,imode
       include "mint.inc"
       real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)
-     $     ,ymax(nintervals,ndimmax,maxchannels),ymax_virt(maxchannels)
-     $     ,x(ndimmax)
+     $     ,ymax(nintervals,ndimmax,maxchannels)
+     $     ,ymax_virt(0:maxchannels),x(ndimmax)
       real * 8 dx(ndimmax),xx(ndimmax),vol_chan,dummy
       integer icell(ndimmax),ncell(ndimmax),ncell_virt
       integer ifold(ndimmax),kfold(ndimmax)

=== modified file 'Template/NLO/SubProcesses/montecarlocounter.f'
--- Template/NLO/SubProcesses/montecarlocounter.f	2018-04-16 14:08:47 +0000
+++ Template/NLO/SubProcesses/montecarlocounter.f	2018-12-20 08:18:55 +0000
@@ -1,3231 +1,3259 @@
-      subroutine set_mc_matrices
-      implicit none
-      include "genps.inc"
-      include 'nexternal.inc'
-      include "born_nhel.inc"
-c Nexternal is the number of legs (initial and final) al NLO, while max_bcol
-c is the number of color flows at Born level
-      integer i,j,k,l,k0,mothercol(2),i1(2)
-      integer idup(nexternal-1,maxproc)
-      integer mothup(2,nexternal-1,maxproc)
-      integer icolup(2,nexternal-1,max_bcol)
-      include 'born_leshouche.inc'
-      integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
-      common /MC_info/ ipartners,colorflow
-      integer i_fks,j_fks
-      common/fks_indices/i_fks,j_fks
-      integer fksfather
-      common/cfksfather/fksfather
-      logical notagluon,found
-      common/cnotagluon/notagluon
-      integer nglu,nsngl
-      logical isspecial,isspecial0
-      common/cisspecial/isspecial
-      logical spec_case
-
-      include 'orders.inc'
-      logical split_type(nsplitorders) 
-      common /c_split_type/split_type
-      double precision particle_charge(nexternal)
-      common /c_charges/particle_charge
-
-c
-      logical is_leading_cflow(max_bcol)
-      integer num_leading_cflows
-      common/c_leading_cflows/is_leading_cflow,num_leading_cflows
-      include 'born_conf.inc'
-      include 'born_coloramps.inc'
-c
-      ipartners(0)=0
-      do i=1,nexternal-1
-         colorflow(i,0)=0
-      enddo
-
-C What follows is true for QCD-type splittings.
-C For QED-type splittings, ipartner is simply all the charged particles
-C in the event except for FKSfather. In this case, all the born color
-C flows are allowed
-
-c ipartners(0): number of particles that can be colour or anticolour partner 
-c   of the father, the Born-level particle to which i_fks and j_fks are 
-c   attached. If one given particle is the colour/anticolour partner of
-c   the father in more than one colour flow, it is counted only once
-c   in ipartners(0)
-c ipartners(i), 1<=i<=nexternal-1: the label (according to Born-level
-c   labelling) of the i^th colour partner of the father
-c
-c colorflow(i,0), 1<=i<=nexternal-1: number of colour flows in which
-c   the particle ipartners(i) is a colour partner of the father
-c colorflow(i,j): the actual label (according to born_leshouche.inc)
-c   of the j^th colour flow in which the father and ipartners(i) are
-c   colour partners
-c
-c Example: in the process q(1) qbar(2) -> g(3) g(4), the two color flows are
-c
-c j=1    i    icolup(1)    icolup(2)       j=2    i    icolup(1)    icolup(2)
-c        1      500           0                   1      500           0
-c        2       0           501                  2       0           501
-c        3      500          502                  3      502          501
-c        4      502          501                  4      500          502
-c
-c and if one fixes for example fksfather=3, then ipartners = 3
-c while colorflow =  0  0  0                                 1
-c                    1  1  0                                 4
-c                    2  1  2                                 2
-c                    1  2  0
-c
-
-      fksfather=min(i_fks,j_fks)
-
-c isspecial will be set equal to .true. only if the father is a gluon,
-c and another gluon will be found which is connected to it by both
-c colour and anticolour
-      isspecial=.false.
-c
-      if (split_type(qcd_pos)) then
-        ! identify the color partners 
-c consider only leading colour flows
-        num_leading_cflows=0
-        do i=1,max_bcol
-          is_leading_cflow(i)=.false.
-          do j=1,mapconfig(0)
-            if(icolamp(i,j,1))then
-               is_leading_cflow(i)=.true.
-               num_leading_cflows=num_leading_cflows+1
-               exit
-            endif
-          enddo
-        enddo
-c
-        do i=1,max_bcol
-          if(.not.is_leading_cflow(i))cycle
-c Loop over Born-level colour flows
-          isspecial0=.false.
-c nglu and nsngl are the number of gluons (except for the father) and of 
-c colour singlets in the Born process, according to the information 
-c stored in ICOLUP
-          nglu=0
-          nsngl=0
-          mothercol(1)=ICOLUP(1,fksfather,i)
-          mothercol(2)=ICOLUP(2,fksfather,i)
-          notagluon=(mothercol(1).eq.0 .or. mothercol(2).eq.0)
-c
-          do j=1,nexternal-1
-c Loop over Born-level particles; j is the possible colour partner of father,
-c and whether this is the case is determined inside this loop
-            if (j.ne.fksfather) then
-c Skip father (it cannot be its own colour partner)
-               if(ICOLUP(1,j,i).eq.0.and.ICOLUP(2,j,i).eq.0)
-     #           nsngl=nsngl+1
-               if(ICOLUP(1,j,i).ne.0.and.ICOLUP(2,j,i).ne.0)
-     #           nglu=nglu+1
-               if ( (j.le.nincoming.and.fksfather.gt.nincoming) .or.
-     #              (j.gt.nincoming.and.fksfather.le.nincoming) ) then
-c father and j not both in the initial or in the final state -- connect
-c colour (1) with colour (i1(1)), and anticolour (2) with anticolour (i1(2))
-                  i1(1)=1
-                  i1(2)=2
-               else
-c father and j both in the initial or in the final state -- connect
-c colour (1) with anticolour (i1(2)), and anticolour (2) with colour (i1(1))
-                  i1(1)=2
-                  i1(2)=1
-               endif
-               do l=1,2
-c Loop over colour and anticolour of father
-                  found=.false.
-                  if( ICOLUP(i1(l),j,i).eq.mothercol(l) .and.
-     &                ICOLUP(i1(l),j,i).ne.0 ) then
-c When ICOLUP(i1(l),j,i) = mothercol(l), the colour (if i1(l)=1) or
-c the anticolour (if i1(l)=2) of particle j is connected to the
-c colour (if l=1) or the anticolour (if l=2) of the father
-                     k0=-1
-                     do k=1,ipartners(0)
-c Loop over previously-found colour/anticolour partners of father
-                        if(ipartners(k).eq.j)then
-                           if(found)then
-c Safety measure: if this condition is met, it means that there exist
-c k1 and k2 such that ipartners(k1)=ipartners(k2). This is thus a bug,
-c since ipartners() is the list of possible partners of father, where each
-c Born-level particle must appears at most once
-                              write(*,*)'Error #1 in set_matrices'
-                              write(*,*)i,j,l,k
-                              stop
-                           endif
-                           found=.true.
-                           k0=k
-                        endif
-                     enddo
-                     if (.not.found) then
-                        ipartners(0)=ipartners(0)+1
-                        ipartners(ipartners(0))=j
-                        k0=ipartners(0)
-                     endif
-c At this point, k0 is the k0^th colour/anticolour partner of father.
-c Therefore, ipartners(k0)=j
-                     if(k0.le.0.or.ipartners(k0).ne.j)then
-                        write(*,*)'Error #2 in set_matrices'
-                        write(*,*)i,j,l,k0,ipartners(k0)
-                        stop
-                     endif
-                     spec_case=l.eq.2 .and. colorflow(k0,0).ge.1 .and.
-     &                    colorflow(k0,colorflow(k0,0)).eq.i 
-                     if (.not.spec_case)then
-c Increase by one the number of colour flows in which the father is
-c (anti)colour-connected with its k0^th partner (according to the
-c list defined by ipartners)
-                        colorflow(k0,0)=colorflow(k0,0)+1
-c Store the label of the colour flow thus found
-                        colorflow(k0,colorflow(k0,0))=i
-                     elseif (spec_case)then
-c Special case: father and ipartners(k0) are both gluons, connected
-c by colour AND anticolour: the number of colour flows was overcounted
-c by one unit, so decrease it
-                         if( notagluon .or.
-     &                       ICOLUP(i1(1),j,i).eq.0 .or.
-     &                       ICOLUP(i1(2),j,i).eq.0 )then
-                            write(*,*)'Error #3 in set_matrices'
-                            write(*,*)i,j,l,k0,i1(1),i1(2)
-                            stop
-                         endif
-                         colorflow(k0,colorflow(k0,0))=i
-                         isspecial0=.true.
-                     endif
-                  endif
-               enddo
-            endif
-          enddo
-          if( ((nglu+nsngl).gt.(nexternal-2)) .or.
-     #       (isspecial0.and.(nglu+nsngl).ne.(nexternal-2)) )then
-           write(*,*)'Error #4 in set_matrices'
-           write(*,*)isspecial0,nglu,nsngl
-           stop
-          endif
-          isspecial=isspecial.or.isspecial0
-        enddo
-
-      else if (split_type(qed_pos)) then
-        ! do nothing, the partner will be assigned at run-time 
-        ! (it is kinematics-dependent)
-        continue
-      endif
-      call check_mc_matrices
-      return
-      end
-
-
-
-      subroutine check_mc_matrices
-      implicit none
-      include "nexternal.inc"
-      include "born_nhel.inc"
-c      include "fks.inc"
-      integer fks_j_from_i(nexternal,0:nexternal)
-     &     ,particle_type(nexternal),pdg_type(nexternal)
-      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
-      integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
-      common /MC_info/ ipartners,colorflow
-      integer i,j,ipart,iflow,ntot,ithere(1000)
-      integer fksfather
-      common/cfksfather/fksfather
-      logical notagluon
-      common/cnotagluon/notagluon
-      logical isspecial
-      common/cisspecial/isspecial
-
-      include 'orders.inc'
-      logical split_type(nsplitorders) 
-      common /c_split_type/split_type
-
-      logical is_leading_cflow(max_bcol)
-      integer num_leading_cflows
-      common/c_leading_cflows/is_leading_cflow,num_leading_cflows
-c
-      if(ipartners(0).gt.nexternal-1)then
-        write(*,*)'Error #1 in check_mc_matrices',ipartners(0)
-        stop
-      endif
-c
-      
-      if (split_type(QCD_pos)) then
-      ! these tests only apply for QCD-type splittings
-        do i=1,ipartners(0)
-          ipart=ipartners(i)
-          if( ipart.eq.fksfather .or.
-     #        ipart.le.0 .or. ipart.gt.nexternal-1 .or.
-     #        ( abs(particle_type(ipart)).ne.3 .and.
-     #          particle_type(ipart).ne.8 ) )then
-            write(*,*)'Error #2 in check_mc_matrices',i,ipart,
-     #  particle_type(ipart)
-            stop
-          endif
-        enddo
-c
-        do i=1,nexternal-1
-          ithere(i)=1
-        enddo
-        do i=1,ipartners(0)
-          ipart=ipartners(i)
-          ithere(ipart)=ithere(ipart)-1
-          if(ithere(ipart).lt.0)then
-            write(*,*)'Error #3 in check_mc_matrices',i,ipart
-            stop
-          endif
-        enddo
-c
-c ntot is the total number of colour plus anticolour partners of father
-        ntot=0
-        do i=1,ipartners(0)
-          ntot=ntot+colorflow(i,0)
-c
-          if( colorflow(i,0).le.0 .or.
-     #        colorflow(i,0).gt.max_bcol )then
-            write(*,*)'Error #4 in check_mc_matrices',i,colorflow(i,0)
-            stop
-          endif
-c
-          do j=1,max_bcol
-            ithere(j)=1
-          enddo
-          do j=1,colorflow(i,0)
-            iflow=colorflow(i,j)
-            ithere(iflow)=ithere(iflow)-1
-            if(ithere(iflow).lt.0)then
-              write(*,*)'Error #5 in check_mc_matrices',i,j,iflow
-              stop
-            endif
-          enddo
-c
-        enddo
-c
-        if( (notagluon.and.ntot.ne.num_leading_cflows) .or.
-     #    ( (.not.notagluon).and.
-     #      ( (.not.isspecial).and.ntot.ne.(2*num_leading_cflows) .or.
-     #        (isspecial.and.ntot.ne.num_leading_cflows) ) ) )then
-          write(*,*)'Error #6 in check_mc_matrices',
-     #     notagluon,ntot,num_leading_cflows,max_bcol
-          stop
-        endif
-c
-        if(num_leading_cflows.gt.max_bcol)then
-          write(*,*)'Error #7 in check_mc_matrices',
-     #     num_leading_cflows,max_bcol
-          stop
-        endif
-
-      else if (split_type(QED_pos)) then
-        ! write here possible checks for QED-type splittings
-        continue
-      endif
-      return
-      end
-
-
-      subroutine find_ipartner_QED(pp, nofpartners)
-      implicit none
-      include 'nexternal.inc'
-      double precision pp(0:3, nexternal)
-      integer nofpartners
-
-      integer fks_j_from_i(nexternal,0:nexternal)
-     &     ,particle_type(nexternal),pdg_type(nexternal)
-      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
-      double precision particle_charge(nexternal)
-      common /c_charges/particle_charge
-
-      integer i_fks,j_fks,fksfather
-      common/fks_indices/i_fks,j_fks
-      double precision pmass(nexternal)
-      double precision zero
-      parameter (zero=0d0)
-
-      include 'genps.inc'
-      include "born_nhel.inc"
-      integer idup(nexternal-1,maxproc)
-      integer mothup(2,nexternal-1,maxproc)
-      integer icolup(2,nexternal-1,max_bcol)
-      include 'born_leshouche.inc'
-      include 'coupl.inc'
-      integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
-      common /MC_info/ ipartners,colorflow
-c
-c     Shower MonteCarlo
-c     
-      character*10 shower_mc
-      common /cMonteCarloType/shower_mc
-
-      logical found
-      logical same_state
-      double precision ppmin, ppnow
-      integer partner
-      integer i,j
-      double precision chargeprod
-      double precision dot
-
-      include 'pmass.inc'
-      
-      found=.false.
-      ppmin=1d99
-      fksfather=min(i_fks,j_fks)
-
-      if (shower_mc.eq.'PYTHIA8') then
-        ! this should follow what is done in TimeShower::setupQEDdip
-        ! first, look for the lowest-mass same- (opposite-)flavour pair of
-        ! particles in the opposite (same) state of the system
-        do j=1,nexternal
-          if (j.ne.fksfather.and.j.ne.i_fks) then
-            same_state = (j.gt.nincoming.and.fksfather.gt.nincoming).or.
-     $                   (j.le.nincoming.and.fksfather.le.nincoming)
-
-            if ((pdg_type(j).eq.pdg_type(fksfather).and..not.same_state).or. 
-     $          (pdg_type(j).eq.-pdg_type(fksfather).and.same_state)) then
-
-              ppnow=dot(pp(0,fksfather),pp(0,j)) - pmass(fksfather)*pmass(j)
-              if (ppnow.lt.ppmin) then
-                found=.true.
-                partner=j
-              endif
-            endif
-          endif
-        enddo
-        
-        ! if no partner has been found, then look for the
-        ! lowest-mass/chargeprod pair
-        if (.not.found) then
-          do j=1,nexternal
-            if (j.ne.fksfather.and.j.ne.i_fks) then
-              if (particle_charge(fksfather).ne.0d0.and.particle_charge(j).ne.0d0) then
-                ppnow=dot(pp(0,fksfather),pp(0,j)) - pmass(fksfather)*pmass(j) / 
-     $            (particle_charge(fksfather) * particle_charge(j))
-                if (ppnow.lt.ppmin) then
-                  found=.true.
-                  partner=j
-                endif
-              endif
-            endif
-          enddo
-        endif
-
-        ! if no partner has been found, then look for the
-        ! lowest-mass pair
-        if (.not.found) then
-          do j=1,nexternal
-            if (j.ne.fksfather.and.j.ne.i_fks) then
-              ppnow=dot(pp(0,fksfather),pp(0,j)) - pmass(fksfather)*pmass(j) 
-              if (ppnow.lt.ppmin) then
-                found=.true.
-                partner=j
-              endif
-            endif
-          enddo
-        endif
-
-      else
-        ! other showers need to be implemented
-        write(*,*) 'ERROR in find_ipartner_QED, not implemented', shower_mc
-        stop 1
-      endif
-
-      if (.not.found) then
-        write(*,*) 'ERROR in find_ipartner_QED, no parthern found'
-        stop 1
-      endif
-
-      ! now, set ipartners
-      ipartners(0) = 1
-      nofpartners = ipartners(0)
-      ipartners(ipartners(0)) = partner
-      ! all color flows have to be included here
-      colorflow(ipartners(0),0)= max_bcol
-      do i = 1, max_bcol
-        colorflow(ipartners(0),i)=i
-      enddo
-      return
-      end
-
-
-      subroutine xmcsubt_wrap(pp,xi_i_fks,y_ij_fks,wgt)
-      implicit none
-      include "nexternal.inc"
-      double precision pp(0:3,nexternal),wgt
-      double precision xi_i_fks,y_ij_fks
-      double precision xmc,xrealme,gfactsf,gfactcl,probne
-      double precision xmcxsec(nexternal),z(nexternal)
-      integer nofpartners
-      logical lzone(nexternal),flagmc
-
-      ! amp split stuff
-      include 'orders.inc'
-      integer iamp
-      double precision amp_split_mc(amp_split_size)
-      common /to_amp_split_mc/amp_split_mc
-      double precision amp_split_gfunc(amp_split_size)
-      common /to_amp_split_gfunc/amp_split_gfunc
-
-
-c True MC subtraction term
-      ! this fills the amp_split_mc
-      ! no need to set them to zero before the call, they are
-      ! reset inside xmcsubt
-      call xmcsubt(pp,xi_i_fks,y_ij_fks,gfactsf,gfactcl,probne,
-     #                   xmc,nofpartners,lzone,flagmc,z,xmcxsec)
-c G-function matrix element, to recover the real soft limit
-      call xmcsubtME(pp,xi_i_fks,y_ij_fks,gfactsf,gfactcl,xrealme)
-      
-      wgt=xmc+xrealme
-      do iamp=1, amp_split_size
-        amp_split_mc(iamp) = amp_split_mc(iamp) + amp_split_gfunc(iamp)
-      enddo
-
-      return
-      end
-
-
-
-      subroutine xmcsubtME(pp,xi_i_fks,y_ij_fks,gfactsf,gfactcl,wgt)
-      implicit none
-      include "nexternal.inc"
-      include "coupl.inc"
-      double precision pp(0:3,nexternal),gfactsf,gfactcl,wgt,wgts,wgtc,wgtsc
-      double precision xi_i_fks,y_ij_fks
-
-      double precision zero,one
-      parameter (zero=0d0)
-      parameter (one=1d0)
-
-      integer izero,ione,itwo
-      parameter (izero=0)
-      parameter (ione=1)
-      parameter (itwo=2)
-
-      double precision p1_cnt(0:3,nexternal,-2:2)
-      double precision wgt_cnt(-2:2)
-      double precision pswgt_cnt(-2:2)
-      double precision jac_cnt(-2:2)
-      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
-
-      double precision xi_i_fks_cnt(-2:2)
-      common /cxiifkscnt/xi_i_fks_cnt
-
-      integer i_fks,j_fks
-      common/fks_indices/i_fks,j_fks
-
-c Particle types (=colour) of i_fks, j_fks and fks_mother
-      integer i_type,j_type,m_type
-      double precision ch_i,ch_j,ch_m
-      common/cparticle_types/i_type,j_type,m_type,ch_i,ch_j,ch_m
-
-      double precision pmass(nexternal)
-      logical is_aorg(nexternal)
-      common /c_is_aorg/is_aorg
-
-      ! amp split stuff
-      include 'orders.inc'
-      integer iamp
-      double precision amp_split_gfunc(amp_split_size)
-      common /to_amp_split_gfunc/amp_split_gfunc
-      double precision amp_split_s(amp_split_size), 
-     $                 amp_split_c(amp_split_size), 
-     $                 amp_split_sc(amp_split_size)
-
-      include "pmass.inc"
-c
-      wgt=0d0
-      do iamp=1, amp_split_size
-        amp_split_gfunc(iamp) = 0d0
-      enddo
-      ! this contribution is needed only for i_fks being a gluon/photon
-      ! (soft limit)
-      if (is_aorg(i_fks))then
-c i_fks is gluon/photon
-         call set_cms_stuff(izero)
-         call sreal(p1_cnt(0,1,0),zero,y_ij_fks,wgts)
-         do iamp=1, amp_split_size
-           amp_split_s(iamp) = amp_split(iamp)
-         enddo
-         call set_cms_stuff(ione)
-         call sreal(p1_cnt(0,1,1),xi_i_fks,one,wgtc)
-         do iamp=1, amp_split_size
-           amp_split_c(iamp) = amp_split(iamp)
-         enddo
-         call set_cms_stuff(itwo)
-         call sreal(p1_cnt(0,1,2),zero,one,wgtsc)
-         do iamp=1, amp_split_size
-           amp_split_sc(iamp) = amp_split(iamp)
-         enddo
-         wgt=wgts+(1-gfactcl)*(wgtc-wgtsc)
-         wgt=wgt*(1-gfactsf)
-         do iamp = 1, amp_split_size
-           amp_split_gfunc(iamp) = amp_split_s(iamp)+(1-gfactcl)*(amp_split_c(iamp)-amp_split_sc(iamp))
-           amp_split_gfunc(iamp) = amp_split_gfunc(iamp)*(1-gfactsf)
-         enddo
-      elseif (abs(i_type).ne.3.and.ch_i.eq.0d0)then
-         ! we should never get here
-         write(*,*) 'FATAL ERROR #1 in xmcsubtME',i_type,i_fks
-         stop
-      endif
-c
-      return
-      end
-
-
-
-c Main routine for MC counterterms
-
-      subroutine xmcsubt(pp,xi_i_fks,y_ij_fks,gfactsf,gfactcl,probne,
-     &                   wgt,nofpartners,lzone,flagmc,z,xmcxsec)
-      implicit none
-      include "nexternal.inc"
-      include "coupl.inc"
-      include "born_nhel.inc"
-      include "fks_powers.inc"
-      include "madfks_mcatnlo.inc"
-      include "run.inc"
-      include "../../Source/MODEL/input.inc"
-      include 'nFKSconfigs.inc'
-      include 'orders.inc'
-      logical split_type(nsplitorders) 
-      common /c_split_type/split_type
-      integer fks_j_from_i(nexternal,0:nexternal)
-     &     ,particle_type(nexternal),pdg_type(nexternal)
-      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
-      double precision particle_charge(nexternal)
-      common /c_charges/particle_charge
-
-      double precision pp(0:3,nexternal),gfactsf,gfactcl,probne,wgt
-      double precision xi_i_fks,y_ij_fks,xm12,xm22
-      double precision xmcxsec(nexternal)
-      integer nofpartners
-      logical lzone(nexternal),flagmc,limit,non_limit
-
-      double precision emsca_bare,ptresc,rrnd,ref_scale,
-     & scalemin,scalemax,wgt1,qMC,emscainv,emscafun
-      double precision emscwgt(nexternal),emscav(nexternal)
-      integer jpartner,mpartner
-      logical emscasharp
-
-      double precision shattmp,dot,xkern(2),xkernazi(2),
-     & born_red(nsplitorders), born_red_tilde(nsplitorders)
-      double precision bornbars(max_bcol,nsplitorders), 
-     &                 bornbarstilde(max_bcol,nsplitorders)
-
-      integer i,j,npartner,cflows,ileg,N_p,iord
-      double precision tk,uk,q1q,q2q,E0sq(nexternal),x,yi,yj,xij,ap(2),Q(2),
-     & s,w1,w2,beta,xfact,prefact,kn,knbar,kn0,betae0,betad,betas,gfactazi,
-     & gfunction,bogus_probne_fun,
-     & z(nexternal),xi(nexternal),xjac(nexternal),ztmp,xitmp,xjactmp,
-     & zHW6,xiHW6,xjacHW6_xiztoxy,zHWPP,xiHWPP,xjacHWPP_xiztoxy,zPY6Q,
-     & xiPY6Q,xjacPY6Q_xiztoxy,zPY6PT,xiPY6PT,xjacPY6PT_xiztoxy,zPY8,
-     & xiPY8,xjacPY8_xiztoxy,wcc
-
-      common/cscaleminmax/xm12,ileg
-      double precision veckn_ev,veckbarn_ev,xp0jfks
-      common/cgenps_fks/veckn_ev,veckbarn_ev,xp0jfks
-      double precision p_born(0:3,nexternal-1)
-      common/pborn/p_born
-      integer i_fks,j_fks
-      common/fks_indices/i_fks,j_fks
-      double precision ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
-      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,
-     #                        sqrtshat,shat
-
-      integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
-      common /MC_info/ ipartners,colorflow
-      logical isspecial
-      common/cisspecial/isspecial
-
-      integer fksfather
-      common/cfksfather/fksfather
-
-      logical softtest,colltest
-      common/sctests/softtest,colltest
-
-      double precision emsca
-      common/cemsca/emsca,emsca_bare,emscasharp,scalemin,scalemax
-
-      double precision ran2,iseed
-      external ran2
-      logical extra
-
-c Stuff to be written (depending on AddInfoLHE) onto the LHE file
-      INTEGER NFKSPROCESS
-      COMMON/C_NFKSPROCESS/NFKSPROCESS
-      integer iSorH_lhe,ifks_lhe(fks_configs) ,jfks_lhe(fks_configs)
-     &     ,fksfather_lhe(fks_configs) ,ipartner_lhe(fks_configs)
-      double precision scale1_lhe(fks_configs),scale2_lhe(fks_configs)
-      common/cto_LHE1/iSorH_lhe,ifks_lhe,jfks_lhe,
-     #                fksfather_lhe,ipartner_lhe
-      common/cto_LHE2/scale1_lhe,scale2_lhe
-
-c Radiation hardness needed (pt_hardness) for the theta function
-c Should be zero if there are no jets at the Born
-      double precision shower_S_scale(fks_configs*2)
-     &     ,shower_H_scale(fks_configs*2),ref_H_scale(fks_configs*2)
-     &     ,pt_hardness
-      common /cshowerscale2/shower_S_scale,shower_H_scale,ref_H_scale
-     &     ,pt_hardness
-
-      double precision becl,delta
-c alsf and besf are the parameters that control gfunsoft
-      double precision alsf,besf
-      common/cgfunsfp/alsf,besf
-c alazi and beazi are the parameters that control gfunazi
-      double precision alazi,beazi
-      common/cgfunazi/alazi,beazi
-
-c Particle types (=color) of i_fks, j_fks and fks_mother
-      integer i_type,j_type,m_type
-      double precision ch_i,ch_j,ch_m
-      common/cparticle_types/i_type,j_type,m_type,ch_i,ch_j,ch_m
-
-C stuff to keep track of the individual orders
-      double precision amp_split_bornbars(amp_split_size,max_bcol,nsplitorders),
-     $                 amp_split_bornbarstilde(amp_split_size,max_bcol,nsplitorders)
-      common /to_amp_split_bornbars/amp_split_bornbars,
-     $                              amp_split_bornbarstilde
-      double precision amp_split_bornred(amp_split_size,nsplitorders),
-     $                 amp_split_bornredtilde(amp_split_size,nsplitorders),
-     &                 amp_split_xmcxsec(amp_split_size,nexternal)
-      common /to_amp_split_xmcxsec/amp_split_xmcxsec
-      integer iamp
-      double precision amp_split_mc(amp_split_size)
-      common /to_amp_split_mc/amp_split_mc
-
-      double precision zero,one,tiny,vtiny,ymin
-      parameter (zero=0d0)
-      parameter (one=1d0)
-      parameter (vtiny=1d-10)
-      parameter (ymin=0.9d0)
-
-      double precision pi
-      parameter(pi=3.1415926535897932384626433d0)
-
-      double precision vcf,vtf,vca
-      parameter (vcf=4d0/3d0)
-      parameter (vtf=1d0/2d0)
-      parameter (vca=3d0)
-
-      double precision pmass(nexternal)
-      include "pmass.inc"
-
-c Initialise
-
-      if (split_type(QED_pos)) then
-        ! QED partners are dynamically found
-        call find_ipartner_QED(pp,nofpartners)
-      endif
-
-      do i = 1,nexternal
-        xmcxsec(i) = 0d0
-        do iamp = 1, amp_split_size
-          amp_split_xmcxsec(iamp,i) = 0d0
-        enddo
-      enddo
-      flagmc   = .false.
-      wgt      = 0d0
-      ztmp     = 0d0
-      xitmp    = 0d0
-      xjactmp  = 0d0
-      gfactazi = 0d0
-      do i = 1,2
-        xkern(i)    = 0d0
-        xkernazi(i) = 0d0
-      enddo
-      do i = 1, nsplitorders
-        born_red(i) = 0d0
-        born_red_tilde(i) = 0d0
-      enddo
-      do iamp = 1, amp_split_size
-        amp_split_mc(iamp)=0d0
-      enddo
-      kn       = veckn_ev
-      knbar    = veckbarn_ev
-      kn0      = xp0jfks
-      nofpartners = ipartners(0)
-      tiny = 1d-6
-      if (softtest.or.colltest)tiny = 1d-12
-c Logical variables to control the IR limits:
-c one can remove any reference to xi_i_fks
-      limit = 1-y_ij_fks.lt.tiny .and. xi_i_fks.ge.tiny
-      non_limit = xi_i_fks.ge.tiny
-
-c Discard if unphysical kinematics
-      if(pp(0,1).le.0d0)return
-
-c Determine invariants, ileg, and MC hardness qMC
-      extra=dampMCsubt.or.AddInfoLHE.or.UseSudakov
-      call kinematics_driver(xi_i_fks,y_ij_fks,shat,pp,ileg,
-     &                       xm12,xm22,tk,uk,q1q,q2q,qMC,extra)
-      w1=-q1q+q2q-tk
-      w2=-q2q+q1q-uk
-      if(extra.and.qMC.lt.0d0)then
-         write(*,*)'Error in xmcsubt: qMC=',qMC
-         stop
-      endif
-
-c Check ileg, and special case for PYTHIA6PT
-      if(ileg.lt.0.or.ileg.gt.4)then
-         write(*,*)'Error in xmcsubt: ileg=',ileg
-         stop
-      endif
-      if(ileg.gt.2.and.shower_mc.eq.'PYTHIA6PT')then
-         write(*,*)'FSR not allowed when matching PY6PT'
-         stop
-      endif
-
-c New or standard MC@NLO formulation
-      probne=bogus_probne_fun(qMC)
-      if(.not.UseSudakov)probne=1.d0
-
-c Call barred Born and assign shower scale
-      call get_mbar(pp,y_ij_fks,ileg,bornbars,bornbarstilde)
-      call assign_emsca(pp,xi_i_fks,y_ij_fks)
-
-c Distinguish ISR and FSR
-      if(ileg.le.2)then
-         delta=min(1d0,deltaI)
-         yj=0d0
-         yi=y_ij_fks
-      elseif(ileg.ge.3)then
-         delta=min(1d0,deltaO)
-         yj=y_ij_fks
-         yi=0d0
-      endif
-      x=1-xi_i_fks
-      s=shat
-      xij=2*(1-xm12/s-(1-x))/(2-(1-x)*(1-yj)) 
-
-c G-function parameters 
-      gfactsf=gfunction(x,alsf,besf,2d0)
-      if(abs(i_type).eq.3)gfactsf=1d0
-      becl=-(1d0-ymin)
-      gfactcl=gfunction(y_ij_fks,alsf,becl,1d0)
-      if(alazi.lt.0d0)gfactazi=1-gfunction(y_ij_fks,-alazi,beazi,delta)
-
-c For processes that have jets at the Born level, we need to include a
-c theta-function: The radiation from the shower should always be softer
-c than the jets at the Born, hence no need to include the MC counter
-c terms when the radiation is hard.
-      if(pt_hardness.gt.shower_S_scale(nFKSprocess*2-1))then
-         emsca=2d0*sqrt(ebeam(1)*ebeam(2))
-         return
-      endif
-
-c Shower variables
-      if(shower_mc.eq.'HERWIGPP')then
-         ztmp=zHWPP(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-         xitmp=xiHWPP(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-         xjactmp=xjacHWPP_xiztoxy(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-      elseif(shower_mc.eq.'PYTHIA6Q')then
-         ztmp=zPY6Q(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-         xitmp=xiPY6Q(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-         xjactmp=xjacPY6Q_xiztoxy(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-      elseif(shower_mc.eq.'PYTHIA6PT')then
-         ztmp=zPY6PT(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-         xitmp=xiPY6PT(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-         xjactmp=xjacPY6PT_xiztoxy(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-      elseif(shower_mc.eq.'PYTHIA8')then
-         ztmp=zPY8(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-         xitmp=xiPY8(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-         xjactmp=xjacPY8_xiztoxy(ileg,xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-      endif
-
-c Main loop over colour partners
-      do npartner=1,ipartners(0)
-
-         E0sq(npartner)=dot(p_born(0,fksfather),p_born(0,ipartners(npartner)))
-         if(E0sq(npartner).lt.0d0)then
-            write(*,*)'Error in xmcsubt: negative E0sq'
-            write(*,*)E0sq(npartner),ileg,npartner
-            stop
-         endif
-
-         z(npartner)=ztmp
-         xi(npartner)=xitmp
-         xjac(npartner)=xjactmp
-         if(shower_mc.eq.'HERWIG6')then
-            z(npartner)=zHW6(ileg,E0sq(npartner),xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-            xi(npartner)=xiHW6(ileg,E0sq(npartner),xm12,xm22,shat,x,yi,yj,tk,uk,q1q,q2q)
-            xjac(npartner)=xjacHW6_xiztoxy(ileg,E0sq(npartner),xm12,xm22,shat,x,yi,yj,
-     &                                                                  tk,uk,q1q,q2q)
-         endif
-
-c Compute dead zones
-         call get_dead_zone(ileg,z(npartner),xi(npartner),s,x,yi,xm12,xm22,w1,w2,qMC,
-     &                      scalemax,ipartners(npartner),fksfather,lzone(npartner),wcc)
-
-c Compute MC subtraction terms
-         if(lzone(npartner))then
-            if(.not.flagmc)flagmc=.true.
-            if( (ileg.ge.3 .and. (m_type.eq.8.or.(m_type.eq.1.and.dabs(ch_m).lt.tiny))) .or.
-     &          (ileg.le.2 .and. (j_type.eq.8.or.(j_type.eq.1.and.dabs(ch_j).lt.tiny))) )then
-               if(i_type.eq.8)then
-c g->gg, go->gog (icode=1)
-                  if(ileg.le.2)then
-                     N_p=2
-                     if(isspecial)N_p=1
-                     if(limit)then
-                        xkern(1)=(g**2/N_p)*8*vca*(1-x*(1-x))**2/(s*x**2)
-                        xkernazi(1)=-(g**2/N_p)*16*vca*(1-x)**2/(s*x**2)
-                        xkern(2)=0d0
-                        xkernazi(2)=0d0
-                     elseif(non_limit)then
-                        xfact=(1-yi)*(1-x)/x
-                        prefact=4/(s*N_p)
-                        call AP_reduced(m_type,i_type,ch_m,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                        call Qterms_reduced_spacelike(m_type,i_type,ch_m,ch_i,one,z(npartner),Q)
-                        do i=1,2
-                          Q(i)=Q(i)/(1-z(npartner))
-                          xkernazi(i)=prefact*xfact*xjac(npartner)*Q(i)/xi(npartner)
-                        enddo
-                        if (xkern(2).ne.0d0 .or.xkernazi(2).ne.0d0) then
-                            write(*,*) 'ERROR#1, g->gg splitting QED'
-     %                        //'contributions should be 0', xkern, xkernazi
-                            stop
-                        endif
-                     endif
-c
-                  elseif(ileg.eq.3)then
-                     N_p=2
-                     if(isspecial)N_p=1
-                     if(non_limit)then
-                        xfact=(2-(1-x)*(1-(kn0/kn)*yj))/kn*knbar*(1-x)*(1-yj)
-                        prefact=2/(s*N_p)
-                        call AP_reduced_SUSY(j_type,i_type,ch_j,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                     endif
-c
-                  elseif(ileg.eq.4)then
-                     N_p=2
-                     if(isspecial)N_p=1
-                     if(limit)then
-                        xkern(1)=(g**2/N_p)*( 8*vca*
-     &                       (s**2*(1-(1-x)*x)-s*(1+x)*xm12+xm12**2)**2 )/
-     &                       ( s*(s-xm12)**2*(s*x-xm12)**2 )
-                        xkernazi(1)=-(g**2/N_p)*(16*vca*s*(1-x)**2)/((s-xm12)**2)
-                        xkern(2)=0d0
-                        xkernazi(2)=0d0
-                     elseif(non_limit)then
-                        xfact=(2-(1-x)*(1-yj))/xij*(1-xm12/s)*(1-x)*(1-yj)
-                        prefact=2/(s*N_p)
-                        call AP_reduced(j_type,i_type,ch_j,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                        call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,one,z(npartner),Q)
-                        do i=1,2
-                          Q(i)=Q(i)/(1-z(npartner))
-                          xkernazi(i)=prefact*xfact*xjac(npartner)*Q(i)/xi(npartner)
-                        enddo
-                        if (xkern(2).ne.0d0 .or.xkernazi(2).ne.0d0) then
-                            write(*,*) 'ERROR#1, g->gg splitting QED'
-     %                        //'contributions should be 0', xkern, xkernazi
-                            stop
-                        endif
-                     endif
-                  endif
-               elseif(abs(i_type).eq.3.or.(i_type.eq.1.and.dabs(ch_i).gt.tiny))then
-c g->qq, a->qq, a->ee (icode=2)
-                  if(ileg.le.2)then
-                     N_p=1
-                     if(limit)then
-                        xkern(1)=(g**2/N_p)*4*vtf*(1-x)*((1-x)**2+x**2)/(s*x)
-                        xkern(2)=xkern(1) * dble(gal(1))**2 / g**2 * 
-     &                                       ch_i**2 * abs(i_type) / vtf
-                     elseif(non_limit)then
-                        xfact=(1-yi)*(1-x)/x
-                        prefact=4/(s*N_p)
-                        call AP_reduced(m_type,i_type,ch_m,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                     endif
-c
-                  elseif(ileg.eq.4)then
-                     N_p=2
-                     if(isspecial)N_p=1
-                     if(limit)then
-                        xkern(1)=(g**2/N_p)*( 4*vtf*(1-x)*
-     &                        (s**2*(1-2*(1-x)*x)-2*s*x*xm12+xm12**2) )/
-     &                        ( (s-xm12)**2*(s*x-xm12) )
-                        xkern(2)=xkern(1) * dble(gal(1))**2 / g**2 *
-     &                                      ch_i**2 * abs(i_type) / vtf
-                        xkernazi(1)=(g**2/N_p)*(16*vtf*s*(1-x)**2)/((s-xm12)**2)
-                        xkernazi(2)=xkernazi(1) * dble(gal(1))**2 / g**2 *
-     &                              ch_i**2 * abs(i_type) / vtf
-                     elseif(non_limit)then
-                        xfact=(2-(1-x)*(1-yj))/xij*(1-xm12/s)*(1-x)*(1-yj)
-                        prefact=2/(s*N_p)
-                        call AP_reduced(j_type,i_type,ch_j,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                        call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,one,z(npartner),Q)
-                        do i=1,2
-                          Q(i)=Q(i)/(1-z(npartner))
-                          xkernazi(i)=prefact*xfact*xjac(npartner)*Q(i)/xi(npartner)
-                        enddo
-                     endif
-                  endif
-               else
-                  write(*,*)'Error 1 in xmcsubt: unknown particle type'
-                  write(*,*)i_type
-                  stop
-               endif
-            elseif( (ileg.ge.3 .and. (abs(m_type).eq.3.or.(m_type.eq.1.and.dabs(ch_m).gt.tiny))) .or.
-     &              (ileg.le.2 .and. (abs(j_type).eq.3.or.(j_type.eq.1.and.dabs(ch_j).gt.tiny))) )then
-               if(abs(i_type).eq.3.or.(i_type.eq.1.and.dabs(ch_i).gt.tiny))then
-c q->gq, q->aq, e->ae (icode=3)
-                  if(ileg.le.2)then
-                     N_p=2
-                     if(isspecial)N_p=1
-                     if(limit)then
-                        xkern(1)=(g**2/N_p)*4*vcf*(1-x)*((1-x)**2+1)/(s*x**2)
-                        xkern(2)=xkern(1) * (dble(gal(1))**2 / g**2) * 
-     &                                      (ch_i**2 / vcf)
-                        xkernazi(1)=-(g**2/N_p)*16*vcf*(1-x)**2/(s*x**2)
-                        xkernazi(2)=xkernazi(1) * (dble(gal(1))**2 / g**2) *
-     &                                      (ch_i**2 / vcf)
-                     elseif(non_limit)then
-                        xfact=(1-yi)*(1-x)/x
-                        prefact=4/(s*N_p)
-                        call AP_reduced(m_type,i_type,ch_m,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                        call Qterms_reduced_spacelike(m_type,i_type,ch_m,ch_i,one,z(npartner),Q)
-                        do i=1,2
-                          Q(i)=Q(i)/(1-z(npartner))
-                          xkernazi(i)=prefact*xfact*xjac(npartner)*Q(i)/xi(npartner)
-                        enddo
-                     endif
-c
-                  elseif(ileg.eq.3)then
-                     N_p=1
-                     if(non_limit)then
-                        xfact=(2-(1-x)*(1-(kn0/kn)*yj))/kn*knbar*(1-x)*(1-yj)
-                        prefact=2/(s*N_p)
-                        call AP_reduced(j_type,i_type,ch_j,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                     endif
-c
-                  elseif(ileg.eq.4)then
-                     N_p=1
-                     if(limit)then
-                        xkern(1)=(g**2/N_p)*
-     &                       ( 4*vcf*(1-x)*(s**2*(1-x)**2+(s-xm12)**2) )/
-     &                       ( (s-xm12)*(s*x-xm12)**2 )
-                        xkern(2)=xkern(1) * (dble(gal(1))**2 / g**2) * 
-     &                                      (ch_i**2 / vcf)
-                     elseif(non_limit)then
-                        xfact=(2-(1-x)*(1-yj))/xij*(1-xm12/s)*(1-x)*(1-yj)
-                        prefact=2/(s*N_p)
-                        call AP_reduced(j_type,i_type,ch_j,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                     endif
-                  endif
-               elseif(i_type.eq.8.or.(i_type.eq.1.and.dabs(ch_i).lt.tiny))then
-c q->qg, q->qa, sq->sqg, sq->sqa, e->ea (icode=4)
-                  if(ileg.le.2)then
-                     N_p=1
-                     if(limit)then
-                        xkern(1)=(g**2/N_p)*4*vcf*(1+x**2)/(s*x)
-                        xkern(2)=xkern(1) * (dble(gal(1))**2 / g**2) * 
-     &                                      (ch_m**2 / vcf)
-                     elseif(non_limit)then
-                        xfact=(1-yi)*(1-x)/x
-                        prefact=4/(s*N_p)
-                        call AP_reduced(m_type,i_type,ch_m,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                     endif
-c
-                  elseif(ileg.eq.3)then
-                     N_p=1
-                     if(non_limit)then
-                        xfact=(2-(1-x)*(1-(kn0/kn)*yj))/kn*knbar*(1-x)*(1-yj)
-                        prefact=2/(s*N_p)
-                        if(abs(PDG_type(j_fks)).le.6)then
-                           if(shower_mc.ne.'HERWIGPP')
-     &                     call AP_reduced(j_type,i_type,ch_j,ch_i,one,z(npartner),ap)
-                           if(shower_mc.eq.'HERWIGPP')
-     &                     call AP_reduced_massive(j_type,i_type,ch_j,ch_i,one,z(npartner),
-     &                                                           xi(npartner),xm12,ap)
-                        else
-                           call AP_reduced_SUSY(j_type,i_type,ch_j,ch_i,one,z(npartner),ap)
-                        endif
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                     endif
-c
-                  elseif(ileg.eq.4)then
-                     N_p=1
-                     if(limit)then
-                        xkern(1)=(g**2/N_p)*4*vcf*
-     &                        ( s**2*(1+x**2)-2*xm12*(s*(1+x)-xm12) )/
-     &                        ( s*(s-xm12)*(s*x-xm12) )
-                        xkern(2)=xkern(1) * (dble(gal(1))**2 / g**2) * 
-     &                                      (ch_j**2 / vcf)
-                     elseif(non_limit)then
-                        xfact=(2-(1-x)*(1-yj))/xij*(1-xm12/s)*(1-x)*(1-yj)
-                        prefact=2/(s*N_p)
-                        call AP_reduced(j_type,i_type,ch_j,ch_i,one,z(npartner),ap)
-                        do i=1,2
-                          ap(i)=ap(i)/(1-z(npartner))
-                          xkern(i)=prefact*xfact*xjac(npartner)*ap(i)/xi(npartner)
-                        enddo
-                     endif
-                  endif
-               else
-                  write(*,*)'Error 2 in xmcsubt: unknown particle type'
-                  write(*,*)i_type
-                  stop
-               endif
-            else
-               write(*,*)'Error 3 in xmcsubt: unknown particle type'
-               write(*,*)j_type,i_type
-               stop
-            endif
-
-            if(dampMCsubt)then
-               if(emscasharp)then
-                  if(qMC.le.scalemax)then
-                     emscwgt(npartner)=1d0
-                     emscav(npartner)=emsca_bare
-                  else
-                     emscwgt(npartner)=0d0
-                     emscav(npartner)=scalemax
-                  endif
-               else
-                  ptresc=(qMC-scalemin)/(scalemax-scalemin)
-                  if(ptresc.le.0d0)then
-                     emscwgt(npartner)=1d0
-                     emscav(npartner)=emsca_bare
-                  elseif(ptresc.lt.1d0)then
-                     emscwgt(npartner)=1-emscafun(ptresc,one)
-                     emscav(npartner)=emsca_bare
-                  else
-                     emscwgt(npartner)=0d0
-                     emscav(npartner)=scalemax
-                  endif
-               endif
-            endif
-
-        else
-c Dead zone
-          do i = 1,2
-            xkern(i)=0d0
-            xkernazi(i)=0d0
-          enddo   
-          if(dampMCsubt)then
-              emscav(npartner)=2d0*sqrt(ebeam(1)*ebeam(2))
-              emscwgt(npartner)=0d0
-          endif
-        endif
-
-        !
-        ! remember, we are still in the npartner = 1, ipartners(0) loop
-        !
-        do iamp = 1, amp_split_size
-          do i = 1,nsplitorders
-            amp_split_bornred(iamp,i)=0d0
-            amp_split_bornredtilde(iamp,i)=0d0
-          enddo
-        enddo
-        do i = 1,nsplitorders
-          born_red(i)=0d0
-          born_red_tilde(i)=0d0
-        enddo
-
-        do i=1,2
-           xkern(i)=xkern(i)*gfactsf*wcc
-           xkernazi(i)=xkernazi(i)*gfactazi*gfactsf*wcc
-        enddo
-        do cflows=1,colorflow(npartner,0)
-          do iord = 1, nsplitorders
-            born_red(iord)=born_red(iord)+bornbars(colorflow(npartner,cflows),iord)
-            born_red_tilde(iord)=born_red_tilde(iord)+bornbarstilde(colorflow(npartner,cflows),iord)
-            do iamp=1,amp_split_size
-              amp_split_bornred(iamp,iord)=amp_split_bornred(iamp,iord)+
-     &           amp_split_bornbars(iamp,colorflow(npartner,cflows),iord)
-              amp_split_bornredtilde(iamp,iord)=amp_split_bornredtilde(iamp,iord)+
-     &           amp_split_bornbarstilde(iamp,colorflow(npartner,cflows),iord)
-            enddo
-          enddo
-        enddo
-
-c check the lines below
-        do iord = 1, nsplitorders
-           if (.not.split_type(iord).or.(iord.ne.qed_pos.and.iord.ne.qcd_pos)) cycle
-           if (iord.eq.qcd_pos) then
-             xmcxsec(npartner) =
-     &          xkern(1)*born_red(iord)+xkernazi(1)*born_red_tilde(iord)
-             do iamp=1,amp_split_size
-               amp_split_xmcxsec(iamp,npartner) =
-     &            xkern(1)*amp_split_bornred(iamp,iord) + 
-     &            xkernazi(1)*amp_split_bornredtilde(iamp,iord)
-             enddo
-           else if(iord.eq.qed_pos) then
-             xmcxsec(npartner) =
-     &          xkern(2)*born_red(iord)+xkernazi(2)*born_red_tilde(iord)
-             do iamp=1,amp_split_size
-               amp_split_xmcxsec(iamp,npartner) =
-     &            xkern(2)*amp_split_bornred(iamp,iord) + 
-     &            xkernazi(2)*amp_split_bornredtilde(iamp,iord)
-             enddo
-           endif
-        enddo
-        if (dampMCsubt) then 
-          xmcxsec(npartner)=xmcxsec(npartner)*emscwgt(npartner)
-          do iamp=1,amp_split_size
-            amp_split_xmcxsec(iamp,npartner) = amp_split_xmcxsec(iamp,npartner)*emscwgt(npartner)
-          enddo
-        endif
-
-        wgt=wgt+xmcxsec(npartner)
-        do iamp=1,amp_split_size
-          amp_split_mc(iamp)=amp_split_mc(iamp)+amp_split_xmcxsec(iamp,npartner)
-        enddo
-
-        if(xmcxsec(npartner).lt.0d0)then
-           write(*,*)'Fatal error in xmcsubt'
-           write(*,*)npartner,xmcxsec(npartner)
-           stop
-        endif
-c End of loop over colour partners
-      enddo
-
-c Assign emsca on statistical basis
-      if(dampMCsubt.and.wgt.gt.1d-30)then
-        rrnd=ran2()
-        wgt1=0d0
-        jpartner=0
-        do npartner=1,ipartners(0)
-           if(lzone(npartner).and.jpartner.eq.0)then
-              wgt1=wgt1+xmcxsec(npartner)
-              if(wgt1.ge.rrnd*wgt)then
-                 jpartner=ipartners(npartner)
-                 mpartner=npartner
-              endif
-           endif
-        enddo
-        if(jpartner.eq.0)then
-           write(*,*)'Error in xmcsubt: emsca unweighting failed'
-           stop
-        else
-           emsca=emscav(mpartner)
-        endif
-      endif
-      if(dampMCsubt.and.wgt.lt.1d-30)emsca=scalemax
-
-c Additional information for LHE
-      if(AddInfoLHE)then
-         fksfather_lhe(nFKSprocess)=fksfather
-         if(jpartner.ne.0)then
-            ipartner_lhe(nFKSprocess)=jpartner
-         else
-c min() avoids troubles if ran2()=1
-            ipartner_lhe(nFKSprocess)=min( int(ran2()*ipartners(0))+1,ipartners(0) )
-            ipartner_lhe(nFKSprocess)=ipartners(ipartner_lhe(nFKSprocess))
-         endif
-         scale1_lhe(nFKSprocess)=qMC
-      endif
-
-      if(dampMCsubt)then
-         if(emsca.lt.scalemin)then
-            write(*,*)'Error in xmcsubt: emsca too small'
-            write(*,*)emsca,jpartner,lzone(jpartner)
-            stop
-         endif
-      endif
-
-      do i=1,nexternal
-         if (i.le.ipartners(0)) then 
-           xmcxsec(i)=xmcxsec(i)*probne
-           do iamp = 1, amp_split_size
-             amp_split_xmcxsec(iamp,i)=amp_split_xmcxsec(iamp,i)*probne
-           enddo
-         else
-           xmcxsec(i)=0d0
-           do iamp = 1, amp_split_size
-             amp_split_xmcxsec(iamp,i)=0d0
-           enddo
-         endif
-      enddo
-
-      return
-      end
-
-
-
-      subroutine get_mbar(p,y_ij_fks,ileg,bornbars,bornbarstilde)
-c Computes barred amplitudes (bornbars) squared according
-c to Odagiri's prescription (hep-ph/9806531).
-c Computes barred azimuthal amplitudes (bornbarstilde) with
-c the same method 
-      implicit none
-
-      include "genps.inc"
-      include "nexternal.inc"
-      include "born_nhel.inc"
-      include "orders.inc"
-
-      double precision p(0:3,nexternal)
-      double precision y_ij_fks,bornbars(max_bcol,nsplitorders),
-     &                          bornbarstilde(max_bcol,nsplitorders)
-
-      double precision zero
-      parameter (zero=0.d0)
-      double complex czero
-      parameter (czero=dcmplx(0d0,0d0))
-      double precision p_born_rot(0:3,nexternal-1)
-
-      integer imother_fks,ileg
-
-      double precision p_born(0:3,nexternal-1)
-      common/pborn/p_born
-
-      double Precision amp2(ngraphs), jamp2(0:ncolor)
-      common/to_amps/  amp2,       jamp2
-
-      integer i_fks,j_fks
-      common/fks_indices/i_fks,j_fks
-
-      double precision wgt_born
-      double complex W1(6),W2(6),W3(6),W4(6),Wij_angle,Wij_recta
-      double complex azifact
-
-      double complex xij_aor
-      common/cxij_aor/xij_aor
-
-      double precision sumborn
-      integer i
-
-      double precision vtiny,pi(0:3),pj(0:3),cphi_mother,sphi_mother
-      parameter (vtiny=1d-12)
-      double complex ximag
-      parameter (ximag=(0.d0,1.d0))
-
-      double precision xi_i_fks_ev,y_ij_fks_ev,t
-      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
-
-      double precision cthbe,sthbe,cphibe,sphibe
-      common/cbeangles/cthbe,sthbe,cphibe,sphibe
-
-      logical calculatedBorn
-      common/ccalculatedBorn/calculatedBorn
-      double precision iden_comp
-      common /c_iden_comp/iden_comp
-
-c Particle types (=color) of i_fks, j_fks and fks_mother
-      integer i_type,j_type,m_type
-      double precision ch_i,ch_j,ch_m
-      common/cparticle_types/i_type,j_type,m_type,ch_i,ch_j,ch_m
-
-      double precision born(nsplitorders)
-      double complex borntilde(nsplitorders)
-      logical split_type(nsplitorders) 
-      common /c_split_type/split_type
-      complex*16 ans_cnt(2, nsplitorders), wgt1(2)
-      common /c_born_cnt/ ans_cnt
-      double complex ans_extra_cnt(2,nsplitorders)
-      integer iord, iextra_cnt, isplitorder_born, isplitorder_cnt
-      common /c_extra_cnt/iextra_cnt, isplitorder_born, isplitorder_cnt
-
-      integer iamp
-      double precision amp_split_born(amp_split_size,nsplitorders) 
-      double complex amp_split_borntilde(amp_split_size,nsplitorders)
-      double precision amp_split_bornbars(amp_split_size,max_bcol,nsplitorders),
-     $                 amp_split_bornbarstilde(amp_split_size,max_bcol,nsplitorders)
-      common /to_amp_split_bornbars/amp_split_bornbars,
-     $                              amp_split_bornbarstilde
-
-
-c
-      logical is_leading_cflow(max_bcol)
-      integer num_leading_cflows
-      common/c_leading_cflows/is_leading_cflow,num_leading_cflows
-      
-c
-C BORN/BORNTILDE
-
-C check if momenta have to be rotated
-
-      if ((ileg.eq.1.or.ileg.eq.2).and.(j_fks.eq.2 .and. nexternal-1.ne.3)) then
-c Rotation according to innerpin.m. Use rotate_invar() if a more 
-c general rotation is needed.
-c Exclude 2->1 (at the Born level) processes: matrix elements are
-c independent of the PS point, but non-zero helicity configurations
-c might flip when rotating the momenta.
-        do i=1,nexternal-1
-           p_born_rot(0,i)=p_born(0,i)
-           p_born_rot(1,i)=-p_born(1,i)
-           p_born_rot(2,i)=p_born(2,i)
-           p_born_rot(3,i)=-p_born(3,i)
-        enddo
-        calculatedBorn=.false.
-        call sborn(p_born_rot,wgt_born)
-        if (iextra_cnt.gt.0) call extra_cnt(p_born_rot, iextra_cnt, ans_extra_cnt)
-        calculatedBorn=.false.
-      else
-        call sborn(p_born,wgt_born)
-        if (iextra_cnt.gt.0) call extra_cnt(p_born, iextra_cnt, ans_extra_cnt)
-      endif
-
-      do iord = 1, nsplitorders
-        if (.not.split_type(iord).or.(iord.ne.qed_pos.and.iord.ne.qcd_pos)) cycle
-C check if any extra_cnt is needed
-        if (iextra_cnt.gt.0) then
-            write(*,*) 'FIXEXTRACNTMC'
-            stop
-            if (iord.eq.isplitorder_born) then
-            ! this is the contribution from the born ME
-               wgt1(1) = ans_cnt(1,iord)
-               wgt1(2) = ans_cnt(2,iord)
-            else if (iord.eq.isplitorder_cnt) then
-            ! this is the contribution from the extra cnt
-               wgt1(1) = ans_extra_cnt(1,iord)
-               wgt1(2) = ans_extra_cnt(2,iord)
-            else
-               write(*,*) 'ERROR in sborncol_isr', iord
-               stop
-            endif
-        else
-           wgt1(1) = ans_cnt(1,iord)
-           wgt1(2) = ans_cnt(2,iord)
-        endif
-        if (abs(m_type).eq.3.or.dabs(ch_m).gt.0d0) wgt1(2) = czero
-        born(iord) = dble(wgt1(1))
-        borntilde(iord) = wgt1(2)
-        do iamp=1, amp_split_size
-          amp_split_born(iamp,iord) = dble(amp_split_cnt(iamp,1,iord))
-          if (abs(m_type).eq.3.or.dabs(ch_m).gt.0d0) then
-            amp_split_borntilde(iamp,iord) = czero
-          else
-            amp_split_borntilde(iamp,iord) = amp_split_cnt(iamp,2,iord)
-          endif
-        enddo
-      enddo
-
-c BORN TILDE
-      if(ileg.eq.1.or.ileg.eq.2)then
-c Insert <ij>/[ij] which is not included by sborn()
-        if (1d0-y_ij_fks.lt.vtiny)then
-           azifact=xij_aor
-        else
-           do i=0,3
-              pi(i)=p_i_fks_ev(i)
-              pj(i)=p(i,j_fks)
-           enddo
-           if(j_fks.eq.2)then
-c Rotation according to innerpin.m. Use rotate_invar() if a more 
-c general rotation is needed
-              pi(1)=-pi(1)
-              pi(3)=-pi(3)
-              pj(1)=-pj(1)
-              pj(3)=-pj(3)
-           endif
-           CALL IXXXSO(pi ,ZERO ,+1,+1,W1)        
-           CALL OXXXSO(pj ,ZERO ,-1,+1,W2)        
-           CALL IXXXSO(pi ,ZERO ,-1,+1,W3)        
-           CALL OXXXSO(pj ,ZERO ,+1,+1,W4)        
-           Wij_angle=(0d0,0d0)
-           Wij_recta=(0d0,0d0)
-           do i=1,4
-              Wij_angle = Wij_angle + W1(i)*W2(i)
-              Wij_recta = Wij_recta + W3(i)*W4(i)
-           enddo
-           azifact=Wij_angle/Wij_recta
-        endif
-c Insert the extra factor due to Madgraph convention for polarization vectors
-        if(j_fks.eq.2)then
-           cphi_mother=-1.d0
-           sphi_mother=0.d0
-        else
-           cphi_mother=1.d0
-           sphi_mother=0.d0
-        endif
-        do iord=1, nsplitorders
-           borntilde(iord) = -(cphi_mother+ximag*sphi_mother)**2 *
-     #                borntilde(iord) * dconjg(azifact)
-           do iamp=1, amp_split_size
-             amp_split_borntilde(iamp,iord) = -(cphi_mother+ximag*sphi_mother)**2 *
-     #                amp_split_borntilde(iamp,iord) * dconjg(azifact)
-           enddo
-        enddo
-      elseif(ileg.eq.3.or.ileg.eq.4)then
-         if((abs(j_type).eq.3.or.ch_j.ne.0d0).and.
-     &     (i_type.eq.8.or.i_type.eq.1).and.
-     &     ch_i.eq.0d0)then
-            do iord=1, nsplitorders
-               borntilde(iord)=czero
-               do iamp=1, amp_split_size
-                 amp_split_borntilde(iamp,iord) = czero
-               enddo
-            enddo
-         elseif((m_type.eq.8.or.m_type.eq.1).and.ch_m.eq.0d0)then
-c Insert <ij>/[ij] which is not included by sborn()
-            if(1.d0-y_ij_fks.lt.vtiny)then
-               azifact=xij_aor
-            else
-               do i=0,3
-                  pi(i)=p_i_fks_ev(i)
-                  pj(i)=p(i,j_fks)
-               enddo
-               CALL IXXXSO(pi ,ZERO ,+1,+1,W1)        
-               CALL OXXXSO(pj ,ZERO ,-1,+1,W2)        
-               CALL IXXXSO(pi ,ZERO ,-1,+1,W3)        
-               CALL OXXXSO(pj ,ZERO ,+1,+1,W4)        
-               Wij_angle=(0d0,0d0)
-               Wij_recta=(0d0,0d0)
-               do i=1,4
-                  Wij_angle = Wij_angle + W1(i)*W2(i)
-                  Wij_recta = Wij_recta + W3(i)*W4(i)
-               enddo
-               azifact=Wij_angle/Wij_recta
-            endif
-c Insert the extra factor due to Madgraph convention for polarization vectors
-            imother_fks=min(i_fks,j_fks)
-            call getaziangles(p_born(0,imother_fks),
-     #                           cphi_mother,sphi_mother)
-            do iord=1, nsplitorders
-               borntilde(iord) = -(cphi_mother-ximag*sphi_mother)**2 *
-     #                  borntilde(iord) * azifact
-               do iamp=1, amp_split_size
-                 amp_split_borntilde(iamp,iord) = -(cphi_mother-ximag*sphi_mother)**2 *
-     #                amp_split_borntilde(iamp,iord) * azifact
-               enddo
-            enddo
-         else
-            write(*,*)'FATAL ERROR in get_mbar',
-     #           i_type,j_type,i_fks,j_fks
-            stop
-         endif
-      else
-         write(*,*)'unknown ileg in get_mbar'
-         stop
-      endif
-
-CMZ! this has to be all changed according to the correct jamps
-
-c born is the total born amplitude squared
-      sumborn=0.d0
-      do i=1,max_bcol
-         if(is_leading_cflow(i))sumborn=sumborn+jamp2(i)
-c sumborn is the sum of the leading-color amplitudes squared
-      enddo
-      
-
-c BARRED AMPLITUDES
-      do i=1,max_bcol
-        do iord=1,nsplitorders
-          if (sumborn.ne.0d0.and.is_leading_cflow(i)) then
-            bornbars(i,iord)=jamp2(i)/sumborn * born(iord) *iden_comp
-            do iamp=1,amp_split_size
-              amp_split_bornbars(iamp,i,iord)=jamp2(i)/sumborn * 
-     &                              amp_split_born(iamp,iord) *iden_comp
-            enddo
-          elseif (born(iord).eq.0d0 .or. jamp2(i).eq.0d0
-     &           .or..not.is_leading_cflow(i)) then
-            bornbars(i,iord)=0d0
-            do iamp=1,amp_split_size
-              amp_split_bornbars(iamp,i,iord)=0d0
-            enddo
-          else
-            write (*,*) 'ERROR #1, dividing by zero'
-            stop
-          endif
-          if (sumborn.ne.0d0.and.is_leading_cflow(i)) then
-            bornbarstilde(i,iord)=jamp2(i)/sumborn * dble(borntilde(iord)) *iden_comp
-            do iamp=1,amp_split_size
-              amp_split_bornbarstilde(iamp,i,iord)=jamp2(i)/sumborn * 
-     &                      dble(amp_split_borntilde(iamp,iord)) *iden_comp
-            enddo
-          elseif (borntilde(iord).eq.0d0 .or. jamp2(i).eq.0d0
-     &           .or..not.is_leading_cflow(i)) then
-            bornbarstilde(i,iord)=0d0
-            do iamp=1,amp_split_size
-              amp_split_bornbarstilde(iamp,i,iord)=0d0 
-            enddo
-          else
-            write (*,*) 'ERROR #2, dividing by zero'
-            stop
-          endif      
-c bornbars(i) is the i-th leading-color amplitude squared re-weighted
-c in such a way that the sum of bornbars(i) is born rather than sumborn.
-c the same holds for bornbarstilde(i).
-        enddo
-      enddo
-
-      return
-      end
-
-
-
-      function gfunction(w,alpha,beta,delta)
-c Gets smoothly to 0 as w goes to 1
-      implicit none
-      double precision gfunction,alpha,beta,delta,w,wmin,wg,tt,tmp
-
-      if(beta.lt.0d0)then
-         wmin=0d0
-      else
-         wmin=max(0d0,1d0-delta)
-      endif
-      wg=min(1d0-(1-wmin)*abs(beta),0.99d0)
-      tt=(abs(w)-wg)/(1d0-wg)
-      if(tt.gt.1d0)then
-         write(*,*)'Fatal error in gfunction',tt
-         stop
-      endif
-      tmp=1d0
-      if(alpha.gt.0d0)then
-         if(tt.gt.0d0.and.abs(w).lt.0.99d0)
-     &   tmp=(1-tt)**(2*alpha)/(tt**(2*alpha)+(1-tt)**(2*alpha))
-         if(abs(w).ge.0.99d0)tmp=0d0
-      endif
-      gfunction=tmp
-      return
-      end
-
-
-
-      subroutine kinematics_driver(xi_i_fks,y_ij_fks,sh,pp,ileg,
-     &                             xm12,xm22,xtk,xuk,xq1q,xq2q,qMC,extra)
-c Determines Mandelstam invariants and assigns ileg and shower-damping
-c variable qMC
-      implicit none
-      include "nexternal.inc"
-      include "coupl.inc"
-      include "run.inc"
-      double precision pp(0:3,nexternal),pp_rec(0:3)
-      double precision xi_i_fks,y_ij_fks,xij
-
-      integer ileg,j,i,nfinal
-      double precision xp1(0:3),xp2(0:3),xk1(0:3),xk2(0:3),xk3(0:3)
-      common/cpkmomenta/xp1,xp2,xk1,xk2,xk3
-      double precision sh,xtk,xuk,w1,w2,xq1q,xq2q,xm12,xm22
-      double precision qMC,zPY8,zeta1,zeta2,get_zeta,z,qMCarg,dot
-      logical extra
-      double precision p_born(0:3,nexternal-1)
-      common/pborn/p_born
-      integer fksfather
-      common/cfksfather/fksfather
-      integer i_fks,j_fks
-      common/fks_indices/i_fks,j_fks
-      double precision tiny
-      parameter(tiny=1d-5)
-      double precision zero
-      parameter(zero=0d0)
-
-      integer isqrtneg
-      save isqrtneg
-
-      double precision pmass(nexternal)
-      include "pmass.inc"
-
-c Initialise
-      do i=0,3
-         pp_rec(i)=0d0
-         xp1(i)=0d0
-         xp2(i)=0d0
-         xk1(i)=0d0
-         xk2(i)=0d0
-         xk3(i)=0d0
-      enddo
-      nfinal=nexternal-2
-      xm12=0d0
-      xm22=0d0
-      xq1q=0d0
-      xq2q=0d0
-      qMC=-1d0
-
-c Discard if unphysical FKS variables
-      if(xi_i_fks.lt.0d0.or.xi_i_fks.gt.1d0.or.
-     &   abs(y_ij_fks).gt.1d0)then
-         write(*,*)'Error 0 in kinematics_driver: fks variables'
-         write(*,*)xi_i_fks,y_ij_fks
-         stop
-      endif
-
-c Determine ileg
-c ileg = 1 ==> emission from left     incoming parton
-c ileg = 2 ==> emission from right    incoming parton
-c ileg = 3 ==> emission from massive  outgoing parton
-c ileg = 4 ==> emission from massless outgoing parton
-c Instead of pmass(j_fks), one should use pmass(fksfather), but the
-c kernels where pmass(fksfather) != pmass(j_fks) are non-singular
-      if(fksfather.le.2)then
-        ileg=fksfather
-      elseif(pmass(j_fks).ne.0d0)then
-        ileg=3
-      elseif(pmass(j_fks).eq.0d0)then
-        ileg=4
-      else
-        write(*,*)'Error 1 in kinematics_driver: unknown ileg'
-        write(*,*)ileg,fksfather,pmass(j_fks)
-        stop
-      endif
-
-c Determine and assign momenta:
-c xp1 = incoming left parton  (emitter (recoiler) if ileg = 1 (2))
-c xp2 = incoming right parton (emitter (recoiler) if ileg = 2 (1))
-c xk1 = outgoing parton       (emitter (recoiler) if ileg = 3 (4))
-c xk2 = outgoing parton       (emitter (recoiler) if ileg = 4 (3))
-c xk3 = extra parton          (FKS parton)
-      do j=0,3
-c xk1 and xk2 are never used for ISR
-         xp1(j)=pp(j,1)
-         xp2(j)=pp(j,2)
-         xk3(j)=pp(j,i_fks)
-         if(ileg.gt.2)pp_rec(j)=pp(j,1)+pp(j,2)-pp(j,i_fks)-pp(j,j_fks)
-         if(ileg.eq.3)then
-            xk1(j)=pp(j,j_fks)
-            xk2(j)=pp_rec(j)
-         elseif(ileg.eq.4)then
-            xk1(j)=pp_rec(j)
-            xk2(j)=pp(j,j_fks)
-         endif
-      enddo
-
-c Determine the Mandelstam invariants needed in the MC functions in terms
-c of FKS variables: the argument of MC functions are (p+k)^2, NOT 2 p.k
-c
-c Definitions of invariants in terms of momenta
-c
-c xm12 =     xk1 . xk1
-c xm22 =     xk2 . xk2
-c xtk  = - 2 xp1 . xk3
-c xuk  = - 2 xp2 . xk3
-c xq1q = - 2 xp1 . xk1 + xm12
-c xq2q = - 2 xp2 . xk2 + xm22
-c w1   = + 2 xk1 . xk3        = - xq1q + xq2q - xtk
-c w2   = + 2 xk2 . xk3        = - xq2q + xq1q - xuk
-c xq1c = - 2 xp1 . xk2        = - s - xtk - xq1q + xm12
-c xq2c = - 2 xp2 . xk1        = - s - xuk - xq2q + xm22
-c
-c Parametrisation of invariants in terms of FKS variables
-c
-c ileg = 1
-c xp1  =  sqrt(s)/2 * ( 1 , 0 , 0 , 1 )
-c xp2  =  sqrt(s)/2 * ( 1 , 0 , 0 , -1 )
-c xk3  =  B * ( 1 , 0 , sqrt(1-yi**2) , yi )
-c xk1  =  irrelevant
-c xk2  =  irrelevant
-c yi = y_ij_fks
-c x = 1 - xi_i_fks
-c B = sqrt(s)/2*(1-x)
-c
-c ileg = 2
-c xp1  =  sqrt(s)/2 * ( 1 , 0 , 0 , 1 )
-c xp2  =  sqrt(s)/2 * ( 1 , 0 , 0 , -1 )
-c xk3  =  B * ( 1 , 0 , sqrt(1-yi**2) , -yi )
-c xk1  =  irrelevant
-c xk2  =  irrelevant
-c yi = y_ij_fks
-c x = 1 - xi_i_fks
-c B = sqrt(s)/2*(1-x)
-c
-c ileg = 3
-c xp1  =  sqrt(s)/2 * ( 1 , 0 , sqrt(1-yi**2) , yi )
-c xp2  =  sqrt(s)/2 * ( 1 , 0 , -sqrt(1-yi**2) , -yi )
-c xk1  =  ( sqrt(veckn_ev**2+xm12) , 0 , 0 , veckn_ev )
-c xk2  =  xp1 + xp2 - xk1 - xk3
-c xk3  =  B * ( 1 , 0 , sqrt(1-yj**2) , yj )
-c yj = y_ij_fks
-c yi = irrelevant
-c x = 1 - xi_i_fks
-c veckn_ev is such that xk2**2 = xm22
-c B = sqrt(s)/2*(1-x)
-c azimuth = irrelevant (hence set = 0)
-c
-c ileg = 4
-c xp1  =  sqrt(s)/2 * ( 1 , 0 , sqrt(1-yi**2) , yi )
-c xp2  =  sqrt(s)/2 * ( 1 , 0 , -sqrt(1-yi**2) , -yi )
-c xk1  =  xp1 + xp2 - xk2 - xk3
-c xk2  =  A * ( 1 , 0 , 0 , 1 )
-c xk3  =  B * ( 1 , 0 , sqrt(1-yj**2) , yj )
-c yj = y_ij_fks
-c yi = irrelevant
-c x = 1 - xi_i_fks
-c A = (s*x-xm12)/(sqrt(s)*(2-(1-x)*(1-yj)))
-c B = sqrt(s)/2*(1-x)
-c azimuth = irrelevant (hence set = 0)
-
-      if(ileg.eq.1)then
-         xtk=-sh*xi_i_fks*(1-y_ij_fks)/2
-         xuk=-sh*xi_i_fks*(1+y_ij_fks)/2
-         if(extra)then
-            if(shower_mc.eq.'HERWIG6'  .or.
-     &         shower_mc.eq.'HERWIGPP' )qMC=xi_i_fks/2*sqrt(sh*(1-y_ij_fks**2))
-            if(shower_mc.eq.'PYTHIA6Q' )qMC=sqrt(-xtk)
-            if(shower_mc.eq.'PYTHIA6PT'.or.
-     &         shower_mc.eq.'PYTHIA8'  )qMC=sqrt(-xtk*xi_i_fks)
-         endif
-      elseif(ileg.eq.2)then
-         xtk=-sh*xi_i_fks*(1+y_ij_fks)/2
-         xuk=-sh*xi_i_fks*(1-y_ij_fks)/2
-         if(extra)then
-            if(shower_mc.eq.'HERWIG6'  .or.
-     &         shower_mc.eq.'HERWIGPP' )qMC=xi_i_fks/2*sqrt(sh*(1-y_ij_fks**2))
-            if(shower_mc.eq.'PYTHIA6Q' )qMC=sqrt(-xuk)
-            if(shower_mc.eq.'PYTHIA6PT'.or.
-     &         shower_mc.eq.'PYTHIA8'  )qMC=sqrt(-xuk*xi_i_fks)
-         endif
-      elseif(ileg.eq.3)then
-         xm12=pmass(j_fks)**2
-         xm22=dot(pp_rec,pp_rec)
-         xtk=-2*dot(xp1,xk3)
-         xuk=-2*dot(xp2,xk3)
-         xq1q=-2*dot(xp1,xk1)+xm12
-         xq2q=-2*dot(xp2,xk2)+xm22
-         w1=-xq1q+xq2q-xtk
-         w2=-xq2q+xq1q-xuk
-         if(extra)then
-            if(shower_mc.eq.'HERWIG6'.or.
-     &         shower_mc.eq.'HERWIGPP')then
-               zeta1=get_zeta(sh,w1,w2,xm12,xm22)
-               qMCarg=zeta1*((1-zeta1)*w1-zeta1*xm12)
-               if(qMCarg.lt.0d0.and.qMCarg.ge.-tiny)qMCarg=0d0
-               if(qMCarg.lt.-tiny) then
-                  isqrtneg=isqrtneg+1
-                  write(*,*)'Error 2 in kinematics_driver: negtive sqrt'
-                  write(*,*)qMCarg,isqrtneg
-                  if(isqrtneg.ge.100)stop
-               endif
-               qMC=sqrt(qMCarg)
-            elseif(shower_mc.eq.'PYTHIA6Q')then
-               qMC=sqrt(w1+xm12)
-            elseif(shower_mc.eq.'PYTHIA6PT')then
-               write(*,*)'PYTHIA6PT not available for FSR'
-               stop
-            elseif(shower_mc.eq.'PYTHIA8')then
-               z=zPY8(ileg,xm12,xm22,sh,1d0-xi_i_fks,0d0,y_ij_fks,xtk,xuk,xq1q,xq2q)
-               qMC=sqrt(z*(1-z)*w1)
-            endif
-         endif
-      elseif(ileg.eq.4)then
-         xm12=dot(pp_rec,pp_rec)
-         xm22=0d0
-         xtk=-2*dot(xp1,xk3)
-         xuk=-2*dot(xp2,xk3)
-         xij=2*(1-xm12/sh-xi_i_fks)/(2-xi_i_fks*(1-y_ij_fks))
-         w2=sh*xi_i_fks*xij*(1-y_ij_fks)/2d0
-         xq2q=-sh*xij*(2-dot(xp1,xk2)*4d0/(sh*xij))/2d0
-         xq1q=xuk+xq2q+w2
-         w1=-xq1q+xq2q-xtk
-         if(extra)then
-            if(shower_mc.eq.'HERWIG6'.or.
-     &         shower_mc.eq.'HERWIGPP')then
-               zeta2=get_zeta(sh,w2,w1,xm22,xm12)
-               qMCarg=zeta2*(1-zeta2)*w2
-               if(qMCarg.lt.0d0.and.qMCarg.ge.-tiny)qMCarg=0d0
-               if(qMCarg.lt.-tiny)then
-                  isqrtneg=isqrtneg+1
-                  write(*,*)'Error 3 in kinematics_driver: negtive sqrt'
-                  write(*,*)qMCarg,isqrtneg
-                  if(isqrtneg.ge.100)stop
-               endif
-               qMC=sqrt(qMCarg)
-            elseif(shower_mc.eq.'PYTHIA6Q')then
-               qMC=sqrt(w2)
-            elseif(shower_mc.eq.'PYTHIA6PT')then
-               write(*,*)'PYTHIA6PT not available for FSR'
-               stop
-            elseif(shower_mc.eq.'PYTHIA8')then
-               z=zPY8(ileg,xm12,xm22,sh,1d0-xi_i_fks,0d0,y_ij_fks,xtk,xuk,xq1q,xq2q)
-               qMC=sqrt(z*(1-z)*w2)
-            endif
-         endif
-      else
-         write(*,*)'Error 4 in kinematics_driver: assigned wrong ileg',ileg
-         stop
-      endif
-
-c Checks on invariants
-      call check_invariants(ileg,sh,xtk,xuk,w1,w2,xq1q,xq2q,xm12,xm22)
-
-      return
-      end
-
-
-
-      subroutine check_invariants(ileg,sh,xtk,xuk,w1,w2,xq1q,xq2q,xm12,xm22)
-      implicit none
-      integer ileg
-
-      double precision xp1(0:3),xp2(0:3),xk1(0:3),xk2(0:3),xk3(0:3)
-      common/cpkmomenta/xp1,xp2,xk1,xk2,xk3
-      double precision sh,xtk,xuk,w1,w2,xq1q,xq2q,xm12,xm22
-
-      double precision tiny,dot
-      parameter(tiny=1d-5)
-
-      if(ileg.le.2)then
-         if((abs(xtk+2*dot(xp1,xk3))/sh.ge.tiny).or.
-     &      (abs(xuk+2*dot(xp2,xk3))/sh.ge.tiny))then
-            write(*,*)'Imprecision 1 in check_invariants'
-            write(*,*)abs(xtk+2*dot(xp1,xk3))/sh,
-     &                abs(xuk+2*dot(xp2,xk3))/sh
-            stop
-         endif
-c
-      elseif(ileg.eq.3)then
-         if(sqrt(w1+xm12).ge.sqrt(sh)-sqrt(xm22))then
-            write(*,*)'Imprecision 2a in check_invariants'
-            write(*,*)sqrt(w1),sqrt(sh),xm22
-            stop
-         endif
-         if(((abs(w1-2*dot(xk1,xk3))/sh.ge.tiny)).or.
-     &      ((abs(w2-2*dot(xk2,xk3))/sh.ge.tiny)))then
-            write(*,*)'Imprecision 2b in check_invariants'
-            write(*,*)abs(w1-2*dot(xk1,xk3))/sh,
-     &                abs(w2-2*dot(xk2,xk3))/sh
-            stop
-         endif
-         if(xm12.eq.0d0)then
-            write(*,*)'Error 2c in check_invariants'
-            stop
-         endif
-c
-      elseif(ileg.eq.4)then
-         if(sqrt(w2).ge.sqrt(sh)-sqrt(xm12))then
-            write(*,*)'Imprecision 3a in check_invariants'
-            write(*,*)sqrt(w2),sqrt(sh),xm12
-            stop
-         endif
-         if(((abs(w2-2*dot(xk2,xk3))/sh.ge.tiny)).or.
-     &      ((abs(xq2q+2*dot(xp2,xk2))/sh.ge.tiny)).or.
-     &      ((abs(xq1q+2*dot(xp1,xk1)-xm12)/sh.ge.tiny)))then
-            write(*,*)'Imprecision 3b in check_invariants'
-            write(*,*)abs(w2-2*dot(xk2,xk3))/sh,
-     &                abs(xq2q+2*dot(xp2,xk2))/sh,
-     &                abs(xq1q+2*dot(xp1,xk1)-xm12)/sh
-            stop
-         endif
-         if(xm22.ne.0d0)then
-            write(*,*)'Error 3c in check_invariants'
-            stop
-         endif
-
-      endif
-
-      return
-      end
-
-
-
-c Monte Carlo functions
-c
-c The invariants given in input to these routines follow FNR conventions
-c (i.e., are defined as (p+k)^2, NOT 2 p.k). 
-c The invariants used inside these routines follow MNR conventions
-c (i.e., are defined as -2p.k, NOT (p+k)^2)
-
-c Herwig6
-
-      function zHW6(ileg,e0sq,xm12,xm22,s,x,yi,yj,xtk,xuk,xq1q,xq2q)
-c Shower energy variable
-      implicit none
-      integer ileg
-      double precision zHW6,e0sq,xm12,xm22,s,x,yi,yj,xtk,xuk,xq1q,xq2q,tiny,
-     &ss,w1,w2,tbeta1,zeta1,tbeta2,zeta2,get_zeta,beta,betae0,betad,betas
-      parameter (tiny=1d-5)
-
-      w1=-xq1q+xq2q-xtk
-      w2=-xq2q+xq1q-xuk
-c
-      if(ileg.eq.1)then
-         if(1-x.lt.tiny)then
-            zHW6=1-(1-x)*(s*(1-yi)+4*e0sq*(1+yi))/(8*e0sq)
-         elseif(1-yi.lt.tiny)then
-            zHW6=x-(1-yi)*(1-x)*(s*x**2-4*e0sq)/(8*e0sq)
-         else
-            ss=1-(1+xuk/s)/(e0sq/xtk)
-            if(ss.lt.0d0)goto 999
-            zHW6=2*(e0sq/xtk)*(1-sqrt(ss))
-         endif
-c
-      elseif(ileg.eq.2)then
-         if(1-x.lt.tiny)then
-            zHW6=1-(1-x)*(s*(1-yi)+4*e0sq*(1+yi))/(8*e0sq)
-         elseif(1-yi.lt.tiny)then
-            zHW6=x-(1-yi)*(1-x)*(s*x**2-4*e0sq)/(8*e0sq)
-         else
-            ss=1-(1+xtk/s)/(e0sq/xuk)
-            if(ss.lt.0d0)goto 999
-            zHW6=2*(e0sq/xuk)*(1-sqrt(ss))
-         endif
-c
-      elseif(ileg.eq.3)then
-         if(e0sq.le.(w1+xm12))goto 999
-         if(1-x.lt.tiny)then
-            beta=1-xm12/s
-            betae0=sqrt(1-xm12/e0sq)
-            betad=sqrt((1-(xm12-xm22)/s)**2-(4*xm22/s))
-            betas=1+(xm12-xm22)/s
-            zHW6=1+(1-x)*( s*(yj*betad-betas)/(4*e0sq*(1+betae0))-
-     &                     betae0*(xm12-xm22+s*(1+(1+yj)*betad-betas))/
-     &                     (betad*(xm12-xm22+s*(1+betad))) )
-         else
-            tbeta1=sqrt(1-(w1+xm12)/e0sq)
-            zeta1=get_zeta(s,w1,w2,xm12,xm22)
-            zHW6=1-tbeta1*zeta1-w1/(2*(1+tbeta1)*e0sq)
-         endif
-c
-      elseif(ileg.eq.4)then
-         if(e0sq.le.w2)goto 999
-         if(1-x.lt.tiny)then
-            zHW6=1-(1-x)*( (s-xm12)*(1-yj)/(8*e0sq)+
-     &                     s*(1+yj)/(2*(s-xm12)) )
-         elseif(1-yj.lt.tiny)then
-            zHW6=(s*x-xm12)/(s-xm12)+(1-yj)*(1-x)*(s*x-xm12)*
-     &           ( (s-xm12)**2*(s*(1-2*x)+xm12)+
-     &             4*e0sq*s*(s*x-xm12*(2-x)) )/
-     &           ( 8*e0sq*(s-xm12)**3 )
-         else
-            tbeta2=sqrt(1-w2/e0sq)
-            zeta2=get_zeta(s,w2,w1,xm22,xm12)
-            zHW6=1-tbeta2*zeta2-w2/(2*(1+tbeta2)*e0sq)
-         endif
-c
-      else
-         write(*,*)'zHW6: unknown ileg'
-         stop
-      endif
-
-      if(zHW6.lt.0d0.or.zHW6.gt.1d0)goto 999
-
-      return
- 999  continue
-      zHW6=-1d0
-      return
-      end
-
-
-
-      function xiHW6(ileg,e0sq,xm12,xm22,s,x,yi,yj,xtk,xuk,xq1q,xq2q)
-c Shower evolution variable
-      implicit none
-      integer ileg
-      double precision xiHW6,e0sq,xm12,xm22,s,x,yi,yj,xtk,xuk,xq1q,xq2q,tiny,z,
-     &zHW6,w1,w2,beta,betae0,betad,betas
-      parameter (tiny=1d-5)
-
-      z=zHW6(ileg,e0sq,xm12,xm22,s,x,yi,yj,xtk,xuk,xq1q,xq2q)
-      if(z.lt.0d0)goto 999
-      w1=-xq1q+xq2q-xtk
-      w2=-xq2q+xq1q-xuk
-c
-      if(ileg.eq.1)then
-         if(1-x.lt.tiny)then
-            xiHW6=2*s*(1-yi)/(s*(1-yi)+4*e0sq*(1+yi))
-         elseif(1-yi.lt.tiny)then
-            xiHW6=(1-yi)*s*x**2/(4*e0sq)
-         else
-            xiHW6=2*(1+xuk/(s*(1-z)))
-         endif
-c
-      elseif(ileg.eq.2)then
-         if(1-x.lt.tiny)then
-            xiHW6=2*s*(1-yi)/(s*(1-yi)+4*e0sq*(1+yi))
-         elseif(1-yi.lt.tiny)then
-            xiHW6=(1-yi)*s*x**2/(4*e0sq)
-         else
-            xiHW6=2*(1+xtk/(s*(1-z)))
-         endif
-c
-      elseif(ileg.eq.3)then
-         if(e0sq.le.(w1+xm12))goto 999
-         if(1-x.lt.tiny)then
-            beta=1-xm12/s
-            betae0=sqrt(1-xm12/e0sq)
-            betad=sqrt((1-(xm12-xm22)/s)**2-(4*xm22/s))
-            betas=1+(xm12-xm22)/s
-            xiHW6=( s*(1+betae0)*betad*(xm12-xm22+s*(1+betad))*(yj*betad-betas) )/
-     &            ( -4*e0sq*betae0*(1+betae0)*(xm12-xm22+s*(1+(1+yj)*betad-betas))+
-     &              (s*betad*(xm12-xm22+s*(1+betad))*(yj*betad-betas)) )
-         else
-            xiHW6=w1/(2*z*(1-z)*e0sq)
-         endif
-c
-      elseif(ileg.eq.4)then
-         if(e0sq.le.w2)goto 999
-         if(1-x.lt.tiny)then
-            xiHW6=2*(s-xm12)**2*(1-yj)/( (s-xm12)**2*(1-yj)+4*e0sq*s*(1+yj) )
-         elseif(1-yj.lt.tiny)then
-            xiHW6=(s-xm12)**2*(1-yj)/(4*e0sq*s)
-         else
-            xiHW6=w2/(2*z*(1-z)*e0sq)
-         endif
-c
-      else
-         write(*,*)'xiHW6: unknown ileg'
-         stop
-      endif
-
-      if(xiHW6.lt.0d0)goto 999
-
-      return
- 999  continue
-      xiHW6=-1d0
-      return
-      end
-
-
-
-      function xjacHW6_xiztoxy(ileg,e0sq,xm12,xm22,s,x,yi,yj,xtk,xuk,xq1q,xq2q)
-c Returns the jacobian d(z,xi)/d(x,y), where z and xi are the shower 
-c variables, and x and y are FKS variables
-      implicit none
-      integer ileg
-      double precision xjacHW6_xiztoxy,e0sq,xm12,xm22,s,x,yi,yj,xtk,xuk,xq1q,
-     &xq2q,tiny,tmp,z,zHW6,xi,xiHW6,w1,w2,tbeta1,zeta1,dw1dx,dw2dx,dw1dy,dw2dy,
-     &tbeta2,get_zeta,beta,betae0,betad,betas,eps1,eps2,beta1,beta2,zmo
-      parameter (tiny=1d-5)
-
-      tmp=0d0
-      z=zHW6(ileg,e0sq,xm12,xm22,s,x,yi,yj,xtk,xuk,xq1q,xq2q)
-      xi=xiHW6(ileg,e0sq,xm12,xm22,s,x,yi,yj,xtk,xuk,xq1q,xq2q)
-      if(z.lt.0d0.or.xi.lt.0d0)goto 999
-      w1=-xq1q+xq2q-xtk
-      w2=-xq2q+xq1q-xuk
-c
-      if(ileg.eq.1)then
-         if(1-x.lt.tiny)then
-            tmp=-2*s/(s*(1-yi)+4*(1+yi)*e0sq)
-         elseif(1-yi.lt.tiny)then
-            tmp=-s*x**2/(4*e0sq)
-         else
-            tmp=-s*(1-x)*z**3/(4*e0sq*(1-z)*(xi*(1-z)+z))