Merge lp:~maddevelopers/mg5amcnlo/2.1.2-lhapdf6 into lp:~maddevelopers/mg5amcnlo/2.1.3

[marco zaro]

enable lhapdf6

[Valentin Hirschi (valentin-hirschi)]

Shouldn't this be merged to 2.2.0 instead?

[marco zaro (marco-zaro)]

i don’t know, which one will be the next release?

On 20 Aug 2014, at 11:53, Valentin Hirschi <email address hidden> wrote:

> Shouldn't this be merged to 2.2.0 instead?
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.1.2-lhapdf6/+merge/231518
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/2.1.3.

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

2.1.3 is for bug fixing only. If we need to rush a release for political reason.
So this should be for 2.2.0

Note that I can not test for lhapdf6 since I can't install lhapdf6 (and do not want to change the configuration of the machime for that).

I have test however the lhapdf5 one, and it doesn't work ...
Note sure what the problem is but he code didn't recognize correctly the lhapdf mode.

If I do the following it works nicely:

generate p p > e+ e-
launch
set pdlabel lhapdf
set lhaid 10772

I'm however worry about some misterious print statement which not always appear multiple times.
This might indicate a side effect and/or that the pdf change during the run...

Here is the log where the see the line twice.

INFO: Using LHAPDF v5.9.1 interface for PDFs
write compile file for card: /Users/omatt/Documents/eclipse/2.1.2-lhapdf6/PROC_sm_0/Cards/param_card.dat
run_card missed argument xmtcentral. Takes default: 0.0
run_card missed argument d. Takes default: 1.0
run_card missed argument gridrun. Takes default: False
run_card missed argument fixed_couplings. Takes default: True
Using random number seed offset = 21
INFO: Running Survey
Creating Jobs
Working on SubProcesses
    P0_qq_ll
INFO: Idle: 0, Running: 1, Completed: 0 [ current time: 12h15 ]
DEBUG: Found too many jobs. Recovering
INFO: Idle: 0, Running: 0, Completed: 1 [ 5s ]
INFO: End survey
refine 10000
INFO: compile directory
INFO: Using LHAPDF v5.9.1 interface for PDFs <------ WHY HERE?

I guess the main question is why is there also twice the line:
INFO: compile directory

Now they are for sure a side effect since, If I run the following:

Then the two first run go trough but the latest returns 0 cross-section:
The log indicate the following line:
 Unimplemented distribution= lhapdf
 Implemented are:
So clearly the switch between the mode are not working in both direction.

For lhapdf6 in itself, I will abstain. My understanding was that Paolo was going to check that part, and since I didn't hear about him that's why I didn't include that part in 2.1.2.

Cheers,

Olivier

review: Needs Fixing

[marco zaro (marco-zaro)]

Ciao Olivier,
but are these problems related to the lhapdf6 branch or they are also there in 2.1.2?
The changes i did in lhapdf6 are very minimal…
I’ll have a look anyway.
Can oyu please check that if one does one run with internal pdfs and the next one with lhapdf, the libpdf and libgeneric are correctly recompiled with the lhapdf related wrappers?
It seems not to be the case for madevent runs…
I’ll have a look later.

Cheers,

Marco

On 21 Aug 2014, at 13:49, Olivier Mattelaer <email address hidden> wrote:

> Review: Needs Fixing
>
> Hi,
>
> 2.1.3 is for bug fixing only. If we need to rush a release for political reason.
> So this should be for 2.2.0
>
> Note that I can not test for lhapdf6 since I can't install lhapdf6 (and do not want to change the configuration of the machime for that).
>
> I have test however the lhapdf5 one, and it doesn't work ...
> Note sure what the problem is but he code didn't recognize correctly the lhapdf mode.
>
> If I do the following it works nicely:
>
> generate p p > e+ e-
> launch
> set pdlabel lhapdf
> set lhaid 10772
>
> I'm however worry about some misterious print statement which not always appear multiple times.
> This might indicate a side effect and/or that the pdf change during the run...
>
> Here is the log where the see the line twice.
>
> INFO: Using LHAPDF v5.9.1 interface for PDFs
> write compile file for card: /Users/omatt/Documents/eclipse/2.1.2-lhapdf6/PROC_sm_0/Cards/param_card.dat
> run_card missed argument xmtcentral. Takes default: 0.0
> run_card missed argument d. Takes default: 1.0
> run_card missed argument gridrun. Takes default: False
> run_card missed argument fixed_couplings. Takes default: True
> Using random number seed offset = 21
> INFO: Running Survey
> Creating Jobs
> Working on SubProcesses
> P0_qq_ll
> INFO: Idle: 0, Running: 1, Completed: 0 [ current time: 12h15 ]
> DEBUG: Found too many jobs. Recovering
> INFO: Idle: 0, Running: 0, Completed: 1 [ 5s ]
> INFO: End survey
> refine 10000
> INFO: compile directory
> INFO: Using LHAPDF v5.9.1 interface for PDFs <------ WHY HERE?
>
> I guess the main question is why is there also twice the line:
> INFO: compile directory
>
>
> Now they are for sure a side effect since, If I run the following:
>
> Then the two first run go trough but the latest returns 0 cross-section:
> The log indicate the following line:
> Unimplemented distribution= lhapdf
> Implemented are:
> So clearly the switch between the mode are not working in both direction.
>
>
> For lhapdf6 in itself, I will abstain. My understanding was that Paolo was going to check that part, and since I didn't hear about him that's why I didn't include that part in 2.1.2.
>
> Cheers,
>
> Olivier
>
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.1.2-lhapdf6/+merge/231518
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/2.1.3.

[marco zaro (marco-zaro)]

Ciao Olivier,
I just managed to install lhapdf6 (v 6.1.4) on my mac, with the —disable-python flag when you configure it…
It is working…
I will now update the lhapdf6 branch with the latest 2.2.0 and check out the problems you were mentioning…
Cheers,

Marco

On 21 Aug 2014, at 14:01, marco zaro <email address hidden> wrote:

> Ciao Olivier,
> but are these problems related to the lhapdf6 branch or they are also there in 2.1.2?
> The changes i did in lhapdf6 are very minimal…
> I’ll have a look anyway.
> Can oyu please check that if one does one run with internal pdfs and the next one with lhapdf, the libpdf and libgeneric are correctly recompiled with the lhapdf related wrappers?
> It seems not to be the case for madevent runs…
> I’ll have a look later.
>
> Cheers,
>
> Marco
>
>
>
> On 21 Aug 2014, at 13:49, Olivier Mattelaer <email address hidden> wrote:
>
>> Review: Needs Fixing
>>
>> Hi,
>>
>> 2.1.3 is for bug fixing only. If we need to rush a release for political reason.
>> So this should be for 2.2.0
>>
>> Note that I can not test for lhapdf6 since I can't install lhapdf6 (and do not want to change the configuration of the machime for that).
>>
>> I have test however the lhapdf5 one, and it doesn't work ...
>> Note sure what the problem is but he code didn't recognize correctly the lhapdf mode.
>>
>> If I do the following it works nicely:
>>
>> generate p p > e+ e-
>> launch
>> set pdlabel lhapdf
>> set lhaid 10772
>>
>> I'm however worry about some misterious print statement which not always appear multiple times.
>> This might indicate a side effect and/or that the pdf change during the run...
>>
>> Here is the log where the see the line twice.
>>
>> INFO: Using LHAPDF v5.9.1 interface for PDFs
>> write compile file for card: /Users/omatt/Documents/eclipse/2.1.2-lhapdf6/PROC_sm_0/Cards/param_card.dat
>> run_card missed argument xmtcentral. Takes default: 0.0
>> run_card missed argument d. Takes default: 1.0
>> run_card missed argument gridrun. Takes default: False
>> run_card missed argument fixed_couplings. Takes default: True
>> Using random number seed offset = 21
>> INFO: Running Survey
>> Creating Jobs
>> Working on SubProcesses
>> P0_qq_ll
>> INFO: Idle: 0, Running: 1, Completed: 0 [ current time: 12h15 ]
>> DEBUG: Found too many jobs. Recovering
>> INFO: Idle: 0, Running: 0, Completed: 1 [ 5s ]
>> INFO: End survey
>> refine 10000
>> INFO: compile directory
>> INFO: Using LHAPDF v5.9.1 interface for PDFs <------ WHY HERE?
>>
>> I guess the main question is why is there also twice the line:
>> INFO: compile directory
>>
>>
>> Now they are for sure a side effect since, If I run the following:
>>
>> Then the two first run go trough but the latest returns 0 cross-section:
>> The log indicate the following line:
>> Unimplemented distribution= lhapdf
>> Implemented are:
>> So clearly the switch between the mode are not working in both direction.
>>
>>
>> For lhapdf6 in itself, I will abstain. My understanding was that Paolo was going to check that part, and since I didn't hear about him that's why I didn't include that part in 2.1.2.
>>
>> Cheers,
>>
>> Olivier
>>
>>
>>
>> --
>...

[lp:~maddevelopers/mg5amcnlo/2.1.2-lhapdf6]

271. By marco zaro on 2014-08-21

merged with 2.2.0 r303

[marco zaro (marco-zaro)]

Hi Olivier,
btw, it runs smoothly for madevent, with this output (the lhapdf print statement occurs only once)
Cheers,

Marco

INFO: load configuration from /Users/marcozaro/Physics/MadGraph/2.1.2-lhapdf6/PROC_sm_0/Cards/me5_configuration.txt
INFO: load configuration from /Users/marcozaro/Physics/MadGraph/2.1.2-lhapdf6/input/mg5_configuration.txt
INFO: load configuration from /Users/marcozaro/Physics/MadGraph/2.1.2-lhapdf6/PROC_sm_0/Cards/me5_configuration.txt
generate_events run_01
The following switches determine which programs are run:
 1 Run the pythia shower/hadronization: pythia=NOT INSTALLED
 2 Run PGS as detector simulator: pgs=NOT INSTALLED
 3 Run Delphes as detector simulator: delphes=NOT INSTALLED
 4 Decay particles with the MadSpin module: madspin=OFF
 5 Add weight to events based on coupling parameters: reweight=OFF
  Either type the switch number (1 to 5) to change its default setting,
  or set any switch explicitly (e.g. type 'madspin=ON' at the prompt)
  Type '0', 'auto', 'done' or just press enter when you are done.
 [0, 4, 5, auto, done, madspin=ON, madspin=OFF, madspin, reweight=ON, ... ][60s to answer]
>
Do you want to edit a card (press enter to bypass editing)?
  1 / param : param_card.dat
  2 / run : run_card.dat
 you can also
   - enter the path to a valid card or banner.
   - use the 'set' command to modify a parameter directly.
     The set option works only for param_card and run_card.
     Type 'help set' for more information on this command.
   - call an external program (ASperGE/MadWidth/...).
     Type 'help' for the list of available command
 [0, done, 1, param, 2, run, enter path][60s to answer]
>set pdlabel lhapdf
INFO: modify parameter pdlabel of the run_card.dat to lhapdf
Do you want to edit a card (press enter to bypass editing)?
  1 / param : param_card.dat
  2 / run : run_card.dat
 you can also
   - enter the path to a valid card or banner.
   - use the 'set' command to modify a parameter directly.
     The set option works only for param_card and run_card.
     Type 'help set' for more information on this command.
   - call an external program (ASperGE/MadWidth/...).
     Type 'help' for the list of available command
 [0, done, 1, param, 2, run, enter path]
>
Generating 10000 events with run name run_01
survey run_01
INFO: compile directory
INFO: Using LHAPDF v6.1.4 interface for PDFs
INFO: Trying to download NNPDF23_nlo_as_0119
NNPDF23_nlo_as_0119.tar.gz: 20.6 MB [100.0%]
INFO: NNPDF23_nlo_as_0119 successfully downloaded and stored in /Users/marcozaro/Physics/LHAPDF-6/share/LHAPDF
write compile file for card: /Users/marcozaro/Physics/MadGraph/2.1.2-lhapdf6/PROC_sm_0/Cards/param_card.dat
run_card missed argument xmtcentral. Takes default: 0.0
run_card missed argument d. Takes default: 1.0
run_card missed argument gridrun. Takes default: False
run_card missed argument fixed_couplings. Takes default: True
Using random number seed offset = 21
INFO: Running Survey
Creating Jobs
Working on SubProcesses
    P0_qq_lvl
INFO: Idle: 0, Running: 1, Completed: 0 [ current time: 15h27 ]
...

[Olivier Mattelaer (olivier-mattelaer)]

For the print statement, it was happening only once (the very first time). No idea why...
I guess that we can forget that proble,

I realize that I forget to put the command that I use for the second bug (the critical one).

generate p p > e+ e-
launch
set pdlabel lhapdf
set lhaid 70251
launch
set pdlabel cteq6l1
launch
set pdlabel lhapdf

actually the command crash as well:
generate p p > e+ e-
launch
set pdlabel cteq6l1
launch
set pdlabel lhapd
set lhaid 70251

Cheers,

Olivier

[marco zaro (marco-zaro)]

Ciao Olivier,
this crash is there in 2.1.2 as well
it is related to the fact that if you run first with internal pdfs, then a
second time with the lhapdf ones, some files in Source need to be
recompiled in order to have LHAPDF working.
The error

 Unimplemented distribution= lhapdf
 Implemented are:
mrs02nl,mrs02nn,mrs0119,mrs0117,mrs0121,mrs01_j,mrs99_1,mrs99_2,mrs99_3,mrs99_4,mrs99_5,mrs99_6,mrs99_7,mrs99_8,mrs99_9,mrs9910,mrs9911,mrs9912,mrs98z1,mrs98z2,mrs98z3,mrs98z4,mrs98z5,mrs98ht,mrs98l1,mrs98l2,mrs98l3,mrs98l4,mrs98l5,cteq3_m,cteq3_l,cteq3_d,cteq4_m,cteq4_d,cteq4_l,cteq4a1,cteq4a2,cteq4a3,cteq4a4,cteq4a5,cteq4hj,cteq4lq,cteq5_m,cteq5_d,cteq5_l,cteq5hj,cteq5hq,cteq5f3,cteq5f4,cteq5m1,ctq5hq1,cteq5l1,cteq6_m,cteq6_d,cteq6_l,cteq6l1,nn23lo,nn23lo1,nn23nlo,

comes from PDF/pdfwrap.f, which should not be compiled and linked if one is
using lhapdf.

in the NLO Source makefile, there is 'clean4pdf' which removes the existing
libraries

clean4pdf:

        rm -f ../lib/libpdf.a

        rm -f ../lib/libgeneric.a

and it is not there in the LO one.
It is called at every run, when Source is compiled.

This should fix the madevent behaviour.

Cheers,

Marco

Marco Zaro

On Thu, Aug 21, 2014 at 4:17 PM, Olivier Mattelaer <
<email address hidden>> wrote:

>
> For the print statement, it was happening only once (the very first time).
> No idea why...
> I guess that we can forget that proble,
>
> I realize that I forget to put the command that I use for the second bug
> (the critical one).
>
> generate p p > e+ e-
> launch
> set pdlabel lhapdf
> set lhaid 70251
> launch
> set pdlabel cteq6l1
> launch
> set pdlabel lhapdf
>
> actually the command crash as well:
> generate p p > e+ e-
> launch
> set pdlabel cteq6l1
> launch
> set pdlabel lhapd
> set lhaid 70251
>
>
> Cheers,
>
> Olivier
> --
>
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.1.2-lhapdf6/+merge/231518
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/2.1.3.
>

[lp:~maddevelopers/mg5amcnlo/2.1.2-lhapdf6]

272. By marco zaro on 2014-08-26

the pdf set is installed locally only for v5 (for v6 an error is raised)

[Olivier Mattelaer (olivier-mattelaer)]

Hi Marco,

>it is related to the fact that if you run first with internal pdfs, then a
second time with the lhapdf ones, some files in Source need to be
recompiled in order to have LHAPDF working.

Did you fix this? Tell me when it's done and I will retry.

Cheers,

Olivier

[Olivier Mattelaer (olivier-mattelaer)]

Ok I put this review in rejected, since they are a second proposal directly in 2.2.0

Cheers,

Olivier

Unmerged revisions

275. By Olivier Mattelaer on 2014-09-05

add a nice error message if the pdlabel is not valid (doesn't need to dig into the fortran log file anymore).

274. By marco zaro on 2014-09-05

merged with 2.2.0 r307

273. By marco zaro on 2014-09-05

fixed bug related to madevent

the commands
    generate p p > e+ e-
    launch
    set pdlabel cteq6l1
    launch
    set pdlabel lhapdf
    set lhaid 21100
now work

272. By marco zaro on 2014-08-26

the pdf set is installed locally only for v5 (for v6 an error is raised)

271. By marco zaro on 2014-08-21

merged with 2.2.0 r303

270. By Marco Zaro <email address hidden> on 2014-06-19

fix in the madevent interface, lhapdf6 should be working for LO as well

269. By Marco Zaro <email address hidden> on 2014-06-19

imerged with 2.1.2 r304

------------- th 2.1.2This line and the following will be ignored --------------

modified:
  Template/LO/Source/PDF/pdf_lhapdf.f
  Template/NLO/MCatNLO/srcCommon/myfastjetfortran.cc
  Template/NLO/Source/PDF/pdf_lhapdf.f
  Template/NLO/SubProcesses/cuts.f
  Template/NLO/SubProcesses/fastjet_wrapper.f
  Template/NLO/SubProcesses/fastjetfortran_madfks_core.cc
  Template/NLO/SubProcesses/fastjetfortran_madfks_full.cc
  madgraph/various/cluster.py
unknown:
  --auto
  ME5_debug@
  PROCNLO_loop_sm_0/
  PROCNLO_loop_sm_1/
  PROC_sm_0/
  additional_command
  py.py
  input/mg5_configuration.txt
  models/loop_sm/model.pkl
  models/sm/model.pkl
  vendor/CutTools/includects/
  vendor/StdHEP/log.mcfio.5938
  vendor/StdHEP/log.stdhep.6496
pending merges:
  Marco Zaro 2014-06-19 fix to push 300
    Marco Zaro 2014-06-19 [merge] merged with 2.1.2-fastjet-eta
    Marco Zaro 2014-06-18 compilation errors fixed in
    Marco Zaro 2014-06-17 the fastjet functions now can directly return tje jets within etamax.
    Rikkert Frederix 2014-06-18 [merge] To be sure to have a correct termination when running on a cluster,
    Rikkert Frederix 2014-06-17 To determine if a cluster run finished correctly, also check that the

268. By Marco Zaro <email address hidden> on 2014-06-17

merged with 2.1.2 r301

267. By Marco Zaro <email address hidden> on 2014-06-12

fixed case in srcHerwig/jimmy.f

266. By Marco Zaro <email address hidden> on 2014-06-11

fixed lhapdf related commands not to copy a pdf set if it exists

=== modified file 'Template/LO/Source/make_opts'
--- Template/LO/Source/make_opts	2014-01-07 10:50:42 +0000
+++ Template/LO/Source/make_opts	2014-08-26 07:56:12 +0000
@@ -1,14 +1,15 @@
 # Environment variables
 
 ifeq ($(origin FFLAGS),undefined)
-  FFLAGS= -O -w -fbounds-check
+  FFLAGS= -O -w -fbounds-check -fPIC
   #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
 endif
 FFLAGS+= -ffixed-line-length-132
 
 # Set FC unless it's defined by an environment variable
 ifeq ($(origin FC),default)
-  FC=gfortran
+#  FC=gfortran
+   FC=gfortran
 endif
 
 # Options: dynamic, lhapdf
@@ -17,6 +18,12 @@
 
 # Option dynamic
 
+ifeq ($(UNAME), Darwin)
+  dylibext=dylib
+else
+  dylibext=so  
+endif
+
 ifdef dynamic
   ifeq ($(UNAME), Darwin)
     libext=dylib

=== modified file 'Template/NLO/Cards/shower_card.dat'
--- Template/NLO/Cards/shower_card.dat	2014-07-01 13:57:28 +0000
+++ Template/NLO/Cards/shower_card.dat	2014-08-26 07:56:12 +0000
@@ -1,81 +1,111 @@
-################################################################################
-#                                                                               
-# This file contains the settings for partons showers to be used in aMC@NLO     
-#   mind the format:   variable = name           #         comment              
-#                                                                               
-################################################################################
-#
-# SHOWER SETTINGS:
-#
-nevents    = -1     # number of events to shower. If negative, shower all events
-nsplit_jobs= 1      # number of jobs to be run on parallel or on a cluster (<100)
-combine_td = T      # combine or not the td files if splitting
-ue_enabled = F      # underlying event on (T)/ off (F) (is MSTP(81) 0/1 for PY6)
-pdfcode    = 0      # pdf code: 0 = internal pdf, 1 = same as NLO, other = lhaglue
-hadronize  = T      # hadronization on/off (is MSTP(111) 0/1 for PY6) 
-maxprint   = 2      # maximum number of events to be printed in the log
-maxerrs    = 0.1    # maximum tolerated fraction of errors
-b_stable   = F      # b hadrons are stable
-pi_stable  = T      # pi0's are stable
-wp_stable  = F      # w+'s are stable
-wm_stable  = F      # w-'s are stable
-z_stable   = F      # z0's are stable
-h_stable   = F      # Higgs bosons are stable
-tap_stable = F      # tau+'s are stable
-tam_stable = F      # tau-'s are stable
-mup_stable = F      # mu+'s are stable
-mum_stable = F      # mu-'s are stable
-rnd_seed   = 0      # random seed (0 is default)
-rnd_seed2  = 0      # 2nd random seed (only for HW6, 0 is default)
-lambda_5   = -1     # lambda_5 value, -1 for default
-b_mass     = -1     # b mass, -1 for default
-is_4lep    = F      # true if it is 4 lepton production (only for PY6)
-is_bbar    = F      # true if it is b-b~ production (only for HW6)
-modbos_1   = 5      # dec mode for boson 1 (only for HW6, overridden by DM_ below)
-modbos_2   = 5      # dec mode for boson 2 (only for HW6, overridden by DM_ below)
-################################################################################
-# DECAY CHANNELS
-# Write down the decay channels for the resonances, to be performed by the shower.
-# The syntax (for a two-body decay) is
-# DM_I = M > D1 D2 @ BR @ ME
-# where I = 1, ..., 99, M is the decaying resonance, D1, D2 are the decay products
-# (up to D5 if  such a decay is supported by the shower), BR is the branching ratio
-# (only used by the HERWIG6 shower, ignored otherwise) and ME is the type of matrix
-# element to be used in the decay (only used by HERWIG6, ignored otherwise).
-# BR's are correctly understood by HERWIG6 only if they add up to one and only if
-# no more than three modes are required for a given resonance.
-# ME corresponds to the third entry of subroutine HWMODK, see the relevant manual.
-# Examples of syntax:
-# Z -> e+ e- or mu+ mu- with BR = 0.5 each
-# DM_1 = 23 > -11 11 @ 0.5d0 @ 100
-# DM_2 = 23 > -13 13 @ 0.5d0 @ 100
-# H -> tau+ tau- with BR = 1
-# DM_3 = 25 > -15 15 @ 1.0d0 @ 0
-# t -> nu_e e+ b @ 1d0
-# DM_4 = 6 > 12 -11 5 @ 1d0 @ 100
-# WARNING: for HERWIG6 the order of decay products in >2-body decays IS RELEVANT.
-# WARNING: 1 -> n decays (with n > 2) are handled by PYTHIA6 and PYTHIA8 through
-# a sequence of 1 -> 2 decays.
-################################################################################
+#***********************************************************************
+#                        MadGraph5_aMC@NLO                             *
+#                                                                      *
+#                      shower_card.dat aMC@NLO                         *
+#                                                                      *
+#  This file is used to set the parameters for the shower.             *
+#                                                                      *
+#  Some notation/conventions:                                          *
+#                                                                      *
+#   Lines starting with a hash (#) are info or comments                *
+#                                                                      *
+#   mind the format:   variable    = value     # comment               *
+#***********************************************************************
+#
+#****************
+# Shower settings
+#****************
+#
+#***********************************************************************
+# Number of events, jobs, errors, and random seeds                     *
+#***********************************************************************
+nevents      = -1     # N evts to shower (< 0 = all)
+nsplit_jobs  = 1      # N jobs to run in parallel (< 100!!)
+combine_td   = T      # combine the topdrawer files if nsplit_jobs > 1
+maxprint     = 2      # N evts to print in the log
+maxerrs      = 0.1    # max fraction of errors
+rnd_seed     = 0      # 1st random seed (0 = default)
+rnd_seed2    = 0      # 2nd random seed (0 = default) !ONLY FOR HWERIG6!
+#***********************************************************************
+# PDFs and non-perturbative modelling                                  *
+#***********************************************************************
+pdfcode      = 0      # 0 = internal, 1 = same as NLO, other = lhaglue
+ue_enabled   = F      # underlying event
+hadronize    = T      # hadronisation on/off        !IGNORED BY HERWIG6!
+lambda_5     = -1     # Lambda_5 (< 0 = default)    !IGNORED BY PYTHIA8!
+#***********************************************************************
+# Stable or unstable particles                                         *
+#***********************************************************************
+b_stable     = F      # set B hadrons stable
+pi_stable    = T      # set pi0's stable
+wp_stable    = F      # set w+'s stable
+wm_stable    = F      # set w-'s stable
+z_stable     = F      # set z0's stable
+h_stable     = F      # set Higgs' stable
+tap_stable   = F      # set tau+'s stable
+tam_stable   = F      # set tau-'s stable
+mup_stable   = F      # set mu+'s stable
+mum_stable   = F      # set mu-'s stable
+#***********************************************************************
+# Mass of the b quark                                                  *
+#***********************************************************************
+b_mass       = -1     # if < 0 = read from SubProcesses/MCmasses_*.inc
+#***********************************************************************
+# Special settings                                                     *
+#***********************************************************************
+is_4lep      = F    # T if 4-lepton production        !ONLY FOR PYTHIA6!
+is_bbar      = F    # T if bb~ production             !ONLY FOR HERWIG6!
+#***********************************************************************
+# Decay channels                                                       *
+# Write down the decay channels for the resonances, to be performed by *
+# the shower.                                                          *
+# The syntax (for a two-body decay) is                                 *
+# DM_I = M > D1 D2 @ BR @ ME                                           *
+# where I < 100, M is the decaying resonance, D1, D2 are the decay     *
+# products (up to D5 if such a decay is supported by the shower), BR   *
+# is the branching ratio (only used by the HERWIG6 shower, ignored     *
+# otherwise) and ME is the type of matrix element to be used in the    *
+# decay (only used by HERWIG6, ignored otherwise).                     *
+# BR's are correctly understood by HERWIG6 only if they add up to 1    *
+# and only if no more than three modes are required for a given        *
+# resonance.                                                           *
+# ME corresponds to the third entry of subroutine HWMODK, see the      *
+# relevant manual.                                                     *
+#                                                                      *
+# WARNING: in HERWIG6, the order of decay products in > 2-body decays  *
+# IS RELEVANT.                                                         *
+# WARNING: in PYTHIA6, turning hadronisation off disables top decays   *
+# WARNING: in PYTHIA6 and PYTHIA8, 1 -> n decays (with n > 2) are      *
+# handled through a sequence of 1 -> 2 decays.                         *
+#                                                                      *
+# Examples of syntax:                                                  *
+# Z -> e+ e- or mu+ mu- with BR = 0.5 each                             *
+# DM_1 = 23 > -11 11 @ 0.5d0 @ 100                                     *
+# DM_2 = 23 > -13 13 @ 0.5d0 @ 100                                     *
+# H -> tau+ tau- with BR = 1                                           *
+# DM_3 = 25 > -15 15 @ 1.0d0 @ 0                                       *
+# t -> nu_e e+ b with BR = 1 (HERWIG)                                  *
+# DM_4 = 6 > 12 -11 5 @ 1d0 @ 100                                      *
+# t -> nu_e e+ b with BR = 1 (PYTHIA)                                  *
+# DM_5 = 6 > 24 5 @ 1d0 @ 100                                          *
+# DM_6 = 24 > 12 -11 @ 1d0 @ 100                                       *
+#***********************************************************************
 
-################################################################################
-#
-# EXTRA LIBRARIES/ANALYSES
-# The following lines need to be changed if the user does not want to create the
-# StdHEP/HEPMC file, but to directly run his/her own analyse.
-# Please note that this works only for HW6 and PY6, and that the analysis should
-# be in the HWAnalyzer/ (or PYAnalyzer/) folder. 
-# Please use files in those folders as examples.
-# "None" and an empty value are equivalent.
-#
-################################################################################
-EXTRALIBS   = stdhep Fmcfio         # Needed extra-libraries (not LHAPDF). 
-                                    #   Default: "stdhep Fmcfio"
-EXTRAPATHS  = ../lib                # Path to the extra-libraries. 
-                                    #   Default: "../lib"
-INCLUDEPATHS=                       # Path to the dirs containing header files neede by C++.
-                                    # Directory names are separated by white spaces
-ANALYSE     =                       # User's analysis and histogramming routines 
-                                    # (please use .o as extension and use spaces to separate files)
+#***********************************************************************
+# Extra Libraries/analyses                                             *
+# The following lines need to be changed if the user does not want to  *
+# create a StdHEP/HepMC file, but to directly run an own analysis (to  *
+# be placed in HWAnalyzer or analogous MCatNLO subfolders).            *
+# Please use files in those folders as examples.                       *
+#***********************************************************************
+EXTRALIBS    = stdhep Fmcfio     # Extra-libraries (not LHAPDF) 
+                                 # Default: "stdhep Fmcfio"
+EXTRAPATHS   = ../lib            # Path to the extra-libraries
+                                 # Default: "../lib"
+INCLUDEPATHS =                   # Path to header files needed by c++
+                                 # Dir names separated by white spaces
+ANALYSE      =                   # User's analysis and histogramming
+                                 # routines (please use .o as extension
+                                 # and use spaces to separate files)
 
 

=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIG6.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIG6.Script	2014-06-19 07:02:07 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_HERWIG6.Script	2014-08-26 07:56:12 +0000
@@ -316,7 +316,6 @@
  '$PART1'  '$PART2'               ! hadron types
   $beammom1 $beammom2               ! beam momenta
   $LHSOFT                         ! Underlying event on (.TRUE.) or off (.FALSE.)
-  $MODBOS_1 $MODBOS_2             ! Gauge boson decay modes
 EOF
 if [ $pdftype = 1 ]
 then

=== modified file 'Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f'
--- Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2014-06-06 09:51:20 +0000
+++ Template/NLO/MCatNLO/srcHerwig/madfks_hwdriver.f	2014-08-26 07:56:12 +0000
@@ -69,9 +69,6 @@
       IF(.NOT.LHSOFT)WRITE(*,*)'Underlying event WILL NOT be generated'
       WRITE(*,*)
 C
-      WRITE(*,*)'Enter decay modes for gauge bosons'
-      READ(*,*)MODBOS(1),MODBOS(2)
-C
       IF(LHSOFT.AND.IPDF.NE.1)THEN
          WRITE(*,*)' '
          WRITE(*,*)'Underlying event requires external PDF sets'

=== removed file 'Template/NLO/Source/make_opts'
--- Template/NLO/Source/make_opts	2014-06-24 08:40:18 +0000
+++ Template/NLO/Source/make_opts	1970-01-01 00:00:00 +0000
@@ -1,86 +0,0 @@
-# Environment variables
-
-ifeq ($(origin FFLAGS),undefined)
-  #FFLAGS= -O -w
-  FFLAGS = -O -fno-automatic
-  #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
-endif
-
-ifeq ($(origin CXXFLAGS),undefined)
-  #CXXFLAGS= -g 
-  CXXFLAGS = -O 
-endif
-
-ifeq ($(origin CFLAGS),undefined)
-  #CFLAGS= -g 
-  CFLAGS = -O 
-endif
-
-  
-
-FFLAGS+= -ffixed-line-length-132
-
-# Set FC unless it's defined by an environment variable
-ifeq ($(origin FC),default)
-#  FC=g77
-  FC=gfortran
-endif
-
-# Set CXX unless it's defined by an environment variable
-ifeq ($(origin CXX),default)
-  CXX=g++
-endif
-
-# Options: dynamic, lhapdf
-
-UNAME := $(shell uname)
-
-# Option dynamic
-
-ifdef dynamic
-  ifeq ($(UNAME), Darwin)
-    libext=dylib
-    FFLAGS+= -fno-common
-    LDFLAGS += -bundle
-    define CREATELIB
-      $(FC) -dynamiclib -undefined dynamic_lookup -o $(1) $(2)
-    endef
-  else
-    libext=so
-    FFLAGS+= -fPIC
-    LDFLAGS += -shared
-    define CREATELIB
-      $(FC) $(FFLAGS) $(LDFLAGS) -o $(1) $(2)
-    endef
-  endif
-else
-  libext=a
-  define CREATELIB
-    $(AR) cru $(1) $(2)
-    ranlib $(1)
-  endef
-endif
-
-# Option lhapdf
-
-ifdef lhapdf
-  alfas_functions=alfas_functions_lhapdf
-  lhapdf=-lLHAPDF -lstdc++
-  reweight_xsec_events_pdf_dummy=
-else
-  alfas_functions=alfas_functions
-  lhapdf=
-  reweight_xsec_events_pdf_dummy=reweight_xsec_events_pdf_dummy.o
-endif
-
-# Madloop
-ifeq ($(has_loops),true)
-  libmadloop=libMadLoop.a
-  libcuttools=-lcts
-else
-  libmadloop=
-  libcuttools=
-endif
-
-# OLP (virtuals) necessary link
-libOLP=

=== added file 'Template/NLO/Source/make_opts.inc'
--- Template/NLO/Source/make_opts.inc	1970-01-01 00:00:00 +0000
+++ Template/NLO/Source/make_opts.inc	2014-08-26 07:56:12 +0000
@@ -0,0 +1,86 @@
+# Environment variables
+
+ifeq ($(origin FFLAGS),undefined)
+  #FFLAGS= -O -w
+  FFLAGS = -O -fno-automatic
+  #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
+endif
+
+ifeq ($(origin CXXFLAGS),undefined)
+  #CXXFLAGS= -g 
+  CXXFLAGS = -O 
+endif
+
+ifeq ($(origin CFLAGS),undefined)
+  #CFLAGS= -g 
+  CFLAGS = -O 
+endif
+
+  
+
+FFLAGS+= -ffixed-line-length-132
+
+# Set FC unless it's defined by an environment variable
+ifeq ($(origin FC),default)
+#  FC=g77
+  FC=gfortran
+endif
+
+# Set CXX unless it's defined by an environment variable
+ifeq ($(origin CXX),default)
+  CXX=g++
+endif
+
+# Options: dynamic, lhapdf
+
+UNAME := $(shell uname)
+
+# Option dynamic
+
+ifdef dynamic
+  ifeq ($(UNAME), Darwin)
+    libext=dylib
+    FFLAGS+= -fno-common
+    LDFLAGS += -bundle
+    define CREATELIB
+      $(FC) -dynamiclib -undefined dynamic_lookup -o $(1) $(2)
+    endef
+  else
+    libext=so
+    FFLAGS+= -fPIC
+    LDFLAGS += -shared
+    define CREATELIB
+      $(FC) $(FFLAGS) $(LDFLAGS) -o $(1) $(2)
+    endef
+  endif
+else
+  libext=a
+  define CREATELIB
+    $(AR) cru $(1) $(2)
+    ranlib $(1)
+  endef
+endif
+
+# Option lhapdf
+
+ifdef lhapdf
+  alfas_functions=alfas_functions_lhapdf
+  lhapdf=-lLHAPDF -lstdc++
+  reweight_xsec_events_pdf_dummy=
+else
+  alfas_functions=alfas_functions
+  lhapdf=
+  reweight_xsec_events_pdf_dummy=reweight_xsec_events_pdf_dummy.o
+endif
+
+# Madloop
+ifdef madloop
+  libmadloop=libMadLoop.a
+  libcuttools=-lcts %(link_tir_libs)s
+else
+  libmadloop=
+  libcuttools=
+endif
+
+# OLP (virtuals) necessary link
+libOLP=

=== modified file 'Template/NLO/Source/makefile'
--- Template/NLO/Source/makefile	2014-06-25 23:19:58 +0000
+++ Template/NLO/Source/makefile	2014-08-26 07:56:12 +0000
@@ -5,6 +5,7 @@
 BINDIR= ../bin/
 PDFDIR= ./PDF/
 CUTTOOLSDIR= ./CutTools/
+IREGIDIR= ./IREGI/src/
 STDHEPDIR= ./StdHEP/
 
 PROCESS= hfill.o matrix.o myamp.o 
@@ -80,10 +81,16 @@
 
 CutTools: $(LIBDIR)libcts.a
 
+IREGI: $(LIBDIR)libiregi.a
+
 $(LIBDIR)libcts.a: $(CUTTOOLSDIR)
 	cd $(CUTTOOLSDIR); make
-	ln -sf ../Source/CutTools/includects/libcts.a $(LIBDIR)libcts.a
-	ln -sf ../Source/CutTools/includects/mpmodule.mod $(LIBDIR)mpmodule.mod
+	ln -sf ../Source/$(CUTTOOLSDIR)includects/libcts.a $(LIBDIR)libcts.a
+	ln -sf ../Source/$(CUTTOOLSDIR)includects/mpmodule.mod $(LIBDIR)mpmodule.mod
+
+$(LIBDIR)libiregi.a: $(IREGIDIR)
+	cd $(IREGIDIR); make
+	ln -sf ../Source/$(IREGIDIR)libiregi.a $(LIBDIR)libiregi.a
 
 StdHEP: $(MCATNLOLIBDIR)libstdhep.a  
 
@@ -99,6 +106,9 @@
 cleanCT:
 	cd $(CUTTOOLSDIR); make clean; cd ..
 
+cleanIR:
+	cd $(IREGIDIR); make clean; cd ..
+
 cleanStdHEP:
 	cd $(STDHEPDIR); make clean; cd ..
 

=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
--- Template/NLO/SubProcesses/BinothLHA.f	2013-12-18 07:47:00 +0000
+++ Template/NLO/SubProcesses/BinothLHA.f	2014-08-26 07:56:12 +0000
@@ -28,8 +28,10 @@
       parameter (pi=3.1415926535897932385d0)
       parameter (zero=0d0)
       double precision p(0:3,nexternal-1)
-      double precision virt_wgt,born_wgt,double,single,virt_wgts(3)
-     $     ,virt_wgts_hel(3),born_wgt_recomputed,born_wgt_recomp_direct
+      double precision virt_wgt,born_wgt,double,single
+     $     ,born_wgt_recomputed,born_wgt_recomp_direct
+      double precision, allocatable :: virt_wgts(:,:)
+      double precision, allocatable :: virt_wgts_hel(:,:)
       double precision mu,ao2pi,conversion,alpha_S
       save conversion
       double precision fkssymmetryfactor,fkssymmetryfactorBorn,
@@ -39,6 +41,8 @@
      &                         fkssymmetryfactorDeg,ngluons,nquarks
       logical firsttime,firsttime_conversion
       data firsttime,firsttime_conversion /.true.,.true./
+      logical firsttime_run
+      data firsttime_run /.true./
       double precision qes2
       common /coupl_es/ qes2
       integer nvtozero
@@ -52,16 +56,18 @@
       parameter (fksprefact=.true.)
       integer ret_code
       double precision madfks_single, madfks_double
-      double precision tolerance,prec_found
+      double precision tolerance, prec_found
+      double precision, allocatable :: accuracies(:)
       integer i,j
       integer nbad, nbadmax
       double precision target,ran2,accum
       external ran2
       double precision wgt_hel(max_bhel)
       common/c_born_hel/wgt_hel
+      integer nsqso, MLResArrayDim
 c statistics for MadLoop
       double precision avgPoleRes(2),PoleDiff(2)
-      integer ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1
+      integer ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1(0:9)
       common/ups_stats/ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1
       parameter (nbadmax = 5)
       double precision pmass(nexternal)
@@ -86,30 +92,39 @@
       single  = 0d0
       double  = 0d0
       prec_found = 1.0d0
-      if (firsttime) then
+      if (firsttime_run) then
+         call get_nsqso_loop(nsqso)
+         call get_answer_dimension(MLResArrayDim)
+         allocate(accuracies(0:nsqso))
+         allocate(virt_wgts(0:3,0:MLResArrayDim))
+         allocate(virt_wgts_hel(0:3,0:MLResArrayDim))
 c Make sure that whenever in the initialisation phase, MadLoop calls
 c itself again to perform stability check to make sure no unstable EPS
 c splips unnoticed.
          CALL FORCE_STABILITY_CHECK(.TRUE.)
-
+         firsttime_run = .false.
+      endif
+      if (firsttime) then
          write(*,*) "alpha_s value used for the virtuals"/
      &        /" is (for the first PS point): ", alpha_S
          tolerance=IRPoleCheckThreshold/10d0 ! for the pole check below
-         call sloopmatrix_thres(p, virt_wgts, tolerance, prec_found,
+         call sloopmatrix_thres(p, virt_wgts, tolerance, accuracies,
      $        ret_code)
-         virt_wgt= virt_wgts(1)/dble(ngluons)
-         single  = virt_wgts(2)/dble(ngluons)
-         double  = virt_wgts(3)/dble(ngluons)
+         prec_found = accuracies(0)
+         virt_wgt= virt_wgts(1,0)/dble(ngluons)
+         single  = virt_wgts(2,0)/dble(ngluons)
+         double  = virt_wgts(3,0)/dble(ngluons)
       else
          tolerance=PrecisionVirtualAtRunTime
 c Just set the accuracy found to a positive value as it is not specified
 c once the initial pole check is performed.
          if (mc_hel.eq.0) then
-            call sloopmatrix_thres(p,virt_wgts,tolerance,prec_found
+            call sloopmatrix_thres(p,virt_wgts,tolerance,accuracies
      $           ,ret_code)
-            virt_wgt= virt_wgts(1)/dble(ngluons)
-            single  = virt_wgts(2)/dble(ngluons)
-            double  = virt_wgts(3)/dble(ngluons)
+            prec_found = accuracies(0)            
+            virt_wgt= virt_wgts(1,0)/dble(ngluons)
+            single  = virt_wgts(2,0)/dble(ngluons)
+            double  = virt_wgts(3,0)/dble(ngluons)
          elseif (mc_hel.eq.1) then
 c Use the Born helicity amplitudes to sample the helicities of the
 c virtual as flat as possible
@@ -142,12 +157,13 @@
 c$$$            fillh=.true.
             fillh=.false.
             call sloopmatrixhel_thres(p,hel(ihel),virt_wgts_hel
-     $           ,tolerance,prec_found,ret_code)
-            virt_wgt = virt_wgt + virt_wgts_hel(1)*dble(goodhel(ihel))
-     $           /volh/4d0/dble(ngluons)
-            single   = single   + virt_wgts_hel(2)*dble(goodhel(ihel))
-     $           /volh/4d0/dble(ngluons)
-            double   = double   + virt_wgts_hel(3)*dble(goodhel(ihel))
+     $           ,tolerance,accuracies,ret_code)
+            prec_found = accuracies(0)
+            virt_wgt = virt_wgt + virt_wgts_hel(1,0)*dble(goodhel(ihel))
+     $           /volh/4d0/dble(ngluons)
+            single   = single   + virt_wgts_hel(2,0)*dble(goodhel(ihel))
+     $           /volh/4d0/dble(ngluons)
+            double   = double   + virt_wgts_hel(3,0)*dble(goodhel(ihel))
      $           /volh/4d0/dble(ngluons)
 c Average over initial state helicities (and take the ngluon factor into
 c account)
@@ -200,7 +216,12 @@
             if (mc_hel.ne.0) then
 c Set-up the MC over helicities. This assumes that the 'HelFilter.dat'
 c exists, which should be the case when firsttime is false.
-               open (unit=67,file='HelFilter.dat',status='old',err=201)
+               if (NHelForMCoverHels.lt.0) then
+                   mc_hel=0
+                   goto 203
+               endif
+               open (unit=67,file='../MadLoop5_resources/HelFilter.dat',
+     $   status='old',err=201)
                hel(0)=0
                j=0
 c optimized loop output
@@ -214,8 +235,18 @@
                  endif
                enddo
                goto 203
+201            continue
+               write (*,*) 'Cannot do MC over hel:'/
+     &     /' "HelFilter.dat" does not exist'/
+     &     /' or does not have the correct format.'/
+     $     /' Change NHelForMCoverHels in FKS_params.dat '/
+     &     /'to explicitly summ over them instead.'
+               stop
+c               mc_hel=0
+c               goto 203
+202            continue
 c non optimized loop output
- 202           rewind(67)
+               rewind(67)
                read(67,*,err=201) (include_hel(i),i=1,max_bhel)
                do i=1,max_bhel
                   if (include_hel(i).eq.'T') then
@@ -225,7 +256,7 @@
                      hel(j)=i
                   endif
                enddo
- 203           continue
+203            continue
 c Only do MC over helicities if there are NHelForMCoverHels
 c or more non-zero (independent) helicities
                if (NHelForMCoverHels.eq.-1) then
@@ -290,9 +321,7 @@
       else
          n10=n10+1              ! no known ret_code (10)
       endif
-      if (mod(ret_code,10).ne.0) then
-         n1=n1+1                ! no known ret_code (1)
-      endif
+         n1(mod(ret_code,10))=n1(mod(ret_code,10))+1                ! unit ret code distribution
 
 c Write out the unstable, non-rescued phase-space points (MadLoop return
 c code is in the four hundreds) or the ones that are found by the pole
@@ -342,10 +371,6 @@
      $     .and. ret_code/100.eq.1)virt_wgt=0d0
 
       return
- 201  write (*,*) 'Cannot do MC over hel:'/
-     &     /' "HelFilter.dat" does not exist'/
-     &     /' or does not have the correct format'
-      stop
       end
 
       subroutine BinothLHAInit(filename)

=== modified file 'Template/NLO/SubProcesses/add_write_info.f'
--- Template/NLO/SubProcesses/add_write_info.f	2014-02-19 09:46:05 +0000
+++ Template/NLO/SubProcesses/add_write_info.f	2014-08-26 07:56:12 +0000
@@ -228,7 +228,7 @@
       enddo
 c Assume helicity summed
       do i=1,nexternal
-         jpart(7,i)=0
+         jpart(7,i)=9
       enddo
       if (firsttime2 .and. isum_hel.ne.0) then
          write (*,*) 'WARNING: for writing the events, no helicity '//

=== modified file 'Template/NLO/SubProcesses/ajob_template'
--- Template/NLO/SubProcesses/ajob_template	2014-06-12 08:52:03 +0000
+++ Template/NLO/SubProcesses/ajob_template	2014-08-26 07:56:12 +0000
@@ -7,8 +7,9 @@
 fi
 }
 
-
-
+if [[ -e MadLoop5_resources.tar && ! -e MadLoop5_resources ]]; then
+    tar -xf MadLoop5_resources.tar
+fi
 
 TAGTAGTAGTAGTAGTAGTAG for i in 1 ; do
     
@@ -72,22 +73,17 @@
     fi
     if  [[ ! -e symfact.dat ]] ; then
     	if [[  -e ../symfact.dat ]] ; then
-	    ln -sf ../symfact.dat .
-	else
-	    ln -sf ../../symfact.dat .    
-	fi
+			ln -sf ../symfact.dat .
+		else
+			ln -sf ../../symfact.dat .    
+		fi
     fi
 
-    link1up MadLoopParams.dat
-    link1up HelFilter.dat 
-    link1up LoopFilter.dat
-    link1up ColorDenomFactors.dat
-    link1up HelConfigs.dat
-    link1up ColorNumFactors.dat
     link1up FKS_params.dat
+#   Not necessary anymore
+#    link1up MadLoop5_resources
     link1up OLE_order.olc
     link1up param_card.dat
-    link1up libLHAPDF.so
     
 #   check where is the basic file for the creation of input_app.txt
 #       

=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2013-12-13 08:38:35 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2014-08-26 07:56:12 +0000
@@ -9,10 +9,13 @@
       double precision pi, zero
       parameter (pi=3.1415926535897932385d0)
       parameter (zero = 0d0)
-      integer npoints, npointsChecked,ret_code
+      integer npoints, npointsChecked
       integer i, j, k
-      double precision tolerance
+      integer return_code
+      double precision tolerance, tolerance_default
+      double precision, allocatable :: accuracies(:)
       double precision accuracy
+      parameter (tolerance_default = 1d-5)
       double precision ren_scale, energy
       parameter (ren_scale = 1d2)
       parameter (energy = 1d3)
@@ -23,7 +26,8 @@
       double precision p_born(0:3,nexternal-1)
       common/pborn/p_born
       double precision pswgt
-      double precision virt_wgts(3), fks_double, fks_single
+      double precision fks_double, fks_single
+      double precision, allocatable :: virt_wgts(:,:)
       double precision double, single, finite
       double complex born(2)
       logical calculatedborn
@@ -46,19 +50,21 @@
       include 'run.inc'
       include 'coupl.inc'
       include 'q_es.inc'
+      integer nsqso      
       double precision pmass(nexternal), pmass_rambo(nexternal)
       integer nfail
+      logical first_time
+      data first_time/.TRUE./
       
-c general MadFKS parameters
-      include "FKSParams.inc"
-
 C-----
 C  BEGIN CODE
 C-----  
-c
-c     Read general MadFKS parameters
-c
-      call FKSParamReader(paramFileName,.TRUE.,.FALSE.)
+      if (first_time) then
+          call get_nsqso_loop(nsqso)          
+          allocate(virt_wgts(0:3,0:nsqso))
+          allocate(accuracies(0:nsqso))
+          first_time = .false.
+      endif
 
       call setrun                !Sets up run parameters
       call setpara('param_card.dat')   !Sets up couplings and masses
@@ -71,9 +77,9 @@
       read(*,*) npoints
       write(*,*)'Insert the relative tolerance'
       write(*,*)' A negative number will mean use the default one: ',
-     1 IRPoleCheckThreshold 
+     1 tolerance_default 
       read(*,*) tolerance
-      if (tolerance .le. zero) tolerance = IRPoleCheckThreshold
+      if (tolerance .le. zero) tolerance = tolerance_default
 
       mu_r = ren_scale
       qes2 = ren_scale**2
@@ -160,12 +166,13 @@
           enddo
 
           call sborn(p_born, born)
-          call sloopmatrix_thres(p_born,virt_wgts,tolerance,accuracy
-     $         ,ret_code) 
+          call sloopmatrix_thres(p_born,virt_wgts,tolerance,
+     1 accuracies,return_code)
+          accuracy=accuracies(0)
 
-          finite = virt_wgts(1)/dble(ngluons)
-          single = virt_wgts(2)/dble(ngluons)
-          double = virt_wgts(3)/dble(ngluons)
+          finite = virt_wgts(1,0)/dble(ngluons)
+          single = virt_wgts(2,0)/dble(ngluons)
+          double = virt_wgts(3,0)/dble(ngluons)
 
 C         If MadLoop was still in initialization mode, then skip this
 C         point for the checks

=== removed symlink 'Template/NLO/SubProcesses/check_sa_loop.f'
=== target was u'../../loop_material/StandAlone/SubProcesses/check_sa.f'
=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2014-07-02 13:45:42 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2014-08-26 07:56:12 +0000
@@ -107,7 +107,7 @@
       common/SHevents/Hevents
       character*10 dum
 c statistics for MadLoop      
-      integer ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1
+      integer ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1(0:9)
       common/ups_stats/ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1
 
       double precision virtual_over_born
@@ -148,7 +148,9 @@
       nqdp=0
       nini=0
       n10=0
-      n1=0
+      do i=0,9
+        n1(i)=0
+      enddo
       
       open (unit=lun+1,file='../dname.mg',status='unknown',err=11)
       read (lun+1,'(a130)',err=11,end=11) buf
@@ -320,7 +322,12 @@
          write(*,*)
      &        "  Unknown return code (10):                        ",n10
          write(*,*)
-     &        "  Unknown return code (1):                         ",n1
+     &        "  Unit return code distribution (1):               "
+         do j=0,9
+           if (n1(j).ne.0) then
+              write(*,*) "#Unit ",j," = ",n1(j)
+           endif
+         enddo
       endif
 
       call cpu_time(tAfter)

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2014-07-02 13:45:42 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2014-08-26 07:56:12 +0000
@@ -128,7 +128,7 @@
       common/SHevents/Hevents
       character*10 dum
 c statistics for MadLoop      
-      integer ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1
+      integer ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1(0:9)
       common/ups_stats/ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1
 
 c timing statistics
@@ -161,7 +161,9 @@
       nqdp=0
       nini=0
       n10=0
-      n1=0
+      do i=0,9
+        n1(i)=0
+      enddo
       
       open (unit=lun+1,file='../dname.mg',status='unknown',err=11)
       read (lun+1,'(a130)',err=11,end=11) buf
@@ -513,7 +515,12 @@
          write(*,*)
      &        "  Unknown return code (10):                        ",n10
          write(*,*)
-     &        "  Unknown return code (1):                         ",n1
+     &        "  Unit return code distribution (1):               "
+         do j=0,9
+           if (n1(j).ne.0) then
+              write(*,*) "#Unit ",j," = ",n1(j)
+           endif
+         enddo
       endif
 
       call cpu_time(tAfter)

=== modified file 'Template/NLO/SubProcesses/driver_vegas.f'
--- Template/NLO/SubProcesses/driver_vegas.f	2013-11-28 17:05:38 +0000
+++ Template/NLO/SubProcesses/driver_vegas.f	2014-08-26 07:56:12 +0000
@@ -98,7 +98,7 @@
 
       integer n_mp, n_disc
 c statistics for MadLoop      
-      integer ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1
+      integer ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1(0:9)
       common/ups_stats/ntot,nsun,nsps,nups,neps,n100,nddp,nqdp,nini,n10,n1
 
 c general MadFKS parameters
@@ -126,7 +126,9 @@
       nqdp=0
       nini=0
       n10=0
-      n1=0
+      do i=0,9
+        n1(i)=0
+      enddo
 
       open (unit=lun+1,file='../dname.mg',status='unknown',err=11)
       read (lun+1,'(a130)',err=11,end=11) buf
@@ -311,7 +313,12 @@
          write(*,*)
      &        "  Unknown return code (10):                        ",n10
          write(*,*)
-     &        "  Unknown return code (1):                         ",n1
+     &        "  Unit return code distribution (1):               "
+         do j=0,9
+           if (n1(j).ne.0) then
+              write(*,*) "#Unit ",j," = ",n1(j)
+           endif
+         enddo
       endif
 
       if(savegrid)call initplot_dummy

=== removed file 'Template/NLO/SubProcesses/makefile_loop'
--- Template/NLO/SubProcesses/makefile_loop	2012-11-15 11:00:21 +0000
+++ Template/NLO/SubProcesses/makefile_loop	1970-01-01 00:00:00 +0000
@@ -1,33 +0,0 @@
-include ../../../Source/make_opts
-
-LIBDIR = ../../../lib/
-LOOPLIB= libMadLoop.a
-LINKLIBS =  -L$(LIBDIR) -lcts -ldhelas -lmodel
-LIBS = $(LIBDIR)libcts.$(libext) $(LIBDIR)libdhelas.$(libext)	\
-	$(LIBDIR)libmodel.$(libext)
-PROCESS= loop_matrix.o improve_ps.o born_matrix.o loop_num.o CT_interface.o		\
-		 $(patsubst %.f,%.o,$(wildcard polynomial.f)) \
-		 $(patsubst %.f,%.o,$(wildcard MadLoopParamReader.f)) \
-		 $(patsubst %.f,%.o,$(wildcard helas_calls*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard jamp?_calls_*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard mp_born_amps_and_wfs.f)) \
-		 $(patsubst %.f,%.o,$(wildcard mp_compute_loop_coefs.f)) \
-		 $(patsubst %.f,%.o,$(wildcard mp_helas_calls*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard coef_construction_*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard loop_CT_calls_*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard mp_coef_construction_*.f))
-
-%.o : %.f
-	$(FC) $(FFLAGS) -c $<
-
-all: cleanlib $(LOOPLIB) 
-
-$(LOOPLIB): $(PROCESS)   
-	@ar cr $(LOOPLIB) $(PROCESS)
-	mv $(LOOPLIB) ../
-
-clean: cleanlib
-		@rm -f *.o 
-
-cleanlib:
-	rm -f ../$(LOOPLIB) $(PROG)

=== added file 'Template/NLO/SubProcesses/makefile_loop.inc'
--- Template/NLO/SubProcesses/makefile_loop.inc	1970-01-01 00:00:00 +0000
+++ Template/NLO/SubProcesses/makefile_loop.inc	2014-08-26 07:56:12 +0000
@@ -0,0 +1,36 @@
+include ../../../Source/make_opts
+
+LIBDIR = ../../../lib/
+LOOPLIB= libMadLoop.a
+
+LINKLIBS =  -L$(LIBDIR) -lcts -ldhelas -lmodel %(link_tir_libs)s
+LIBS = $(LIBDIR)libcts.$(libext) $(LIBDIR)libdhelas.$(libext)	\
+	$(LIBDIR)libmodel.$(libext) %(tir_libs)s
+PROCESS= loop_matrix.o improve_ps.o born_matrix.o loop_num.o CT_interface.o	MadLoopCommons.o \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard polynomial.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard MadLoopParamReader.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard helas_calls*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard jamp?_calls_*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_born_amps_and_wfs.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_compute_loop_coefs.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_helas_calls*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard coef_construction_*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard loop_CT_calls_*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_coef_construction_*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard TIR_interface.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard GOLEM_interface.f))		 
+
+%(doto)s : %(dotf)s
+	$(FC) $(FFLAGS) -c $< %(tir_include)s
+
+all: cleanlib $(LOOPLIB) 
+
+$(LOOPLIB): $(PROCESS)   
+	@ar cr $(LOOPLIB) $(PROCESS)
+	mv $(LOOPLIB) ../
+
+clean: cleanlib
+		@rm -f *.o 
+
+cleanlib:
+	rm -f ../$(LOOPLIB) $(PROG)

=== modified file 'Template/NLO/SubProcesses/montecarlocounter.f'
--- Template/NLO/SubProcesses/montecarlocounter.f	2014-05-27 16:58:59 +0000
+++ Template/NLO/SubProcesses/montecarlocounter.f	2014-08-26 07:56:12 +0000
@@ -349,6 +349,7 @@
       include "fks_powers.inc"
       include "madfks_mcatnlo.inc"
       include "run.inc"
+      include "../../Source/MODEL/input.inc"
       include 'nFKSconfigs.inc'
       integer fks_j_from_i(nexternal,0:nexternal)
      &     ,particle_type(nexternal),pdg_type(nexternal)
@@ -366,7 +367,7 @@
       integer jpartner,mpartner
       logical emscasharp
 
-      double precision shattmp,dot,xkern,xkernazi,born_red,
+      double precision shattmp,dot,xkern(2),xkernazi(2),born_red,
      & born_red_tilde
       double precision bornbars(max_bcol), bornbarstilde(max_bcol)
 
@@ -452,6 +453,7 @@
       parameter (vtf=1d0/2d0)
       parameter (vca=3d0)
 
+      double precision g_ew,charge,qi2,qj2
       double precision pmass(nexternal)
       include "pmass.inc"
 
@@ -470,6 +472,9 @@
       knbar    = veckbarn_ev
       kn0      = xp0jfks
       nofpartners = ipartners(0)
+      g_ew=sqrt(4.d0*pi/aewm1)
+      qi2=charge(pdg_type(i_fks))**2
+      qj2=charge(pdg_type(j_fks))**2
       tiny = 1d-4
       if (softtest.or.colltest)tiny = 1d-6
 c Logical variables to control the IR limits:
@@ -585,8 +590,8 @@
 c Compute MC subtraction terms
          if(lzone(npartner))then
             if(.not.flagmc)flagmc=.true.
-            if( (ileg.ge.3 .and. m_type.eq.8) .or.
-     &          (ileg.le.2 .and. j_type.eq.8) )then
+            if( (ileg.ge.3 .and. (m_type.eq.8.or.m_type.eq.0)) .or.
+     &          (ileg.le.2 .and. (j_type.eq.8.or.j_type.eq.0)) )then
                if(i_type.eq.8)then
 c g  --> g  g ( icode = 1 )
 c go --> go g
@@ -594,17 +599,17 @@
                      N_p=2
                      if(isspecial)N_p=1
                      if(limit)then
-                        xkern=(g**2/N_p)*8*vca*(1-x*(1-x))**2/(s*x**2)
-                        xkernazi=-(g**2/N_p)*16*vca*(1-x)**2/(s*x**2)
+                        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)
                      elseif(non_limit)then
                         xfact=(1-yi)*(1-x)/x
                         prefact=4/(s*N_p)
                         call AP_reduced(m_type,i_type,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
                         call Qterms_reduced_spacelike(m_type,i_type,one,z(npartner),Q)
                         Q=Q/(1-z(npartner))
-                        xkernazi=prefact*xfact*xjac(npartner)*Q/xi(npartner)
+                        xkernazi(1)=prefact*xfact*xjac(npartner)*Q/xi(npartner)
                      endif
 c
                   elseif(ileg.eq.3)then
@@ -615,59 +620,66 @@
                         prefact=2/(s*N_p)
                         call AP_reduced_SUSY(j_type,i_type,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
                      endif
 c
                   elseif(ileg.eq.4)then
                      N_p=2
                      if(isspecial)N_p=1
                      if(limit)then
-                        xkern=(g**2/N_p)*( 8*vca*
+                        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=-(g**2/N_p)*(16*vca*s*(1-x)**2)/((s-xm12)**2)
+                        xkernazi(1)=-(g**2/N_p)*(16*vca*s*(1-x)**2)/((s-xm12)**2)
                      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,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
                         call Qterms_reduced_timelike(j_type,i_type,one,z(npartner),Q)
                         Q=Q/(1-z(npartner))
-                        xkernazi=prefact*xfact*xjac(npartner)*Q/xi(npartner)
+                        xkernazi(1)=prefact*xfact*xjac(npartner)*Q/xi(npartner)
                      endif
                   endif
                elseif(abs(i_type).eq.3)then
 c g --> q q~ ( icode = 2 )
+c a --> q q~
                   if(ileg.le.2)then
                      N_p=1
                      if(limit)then
-                        xkern=(g**2/N_p)*4*vtf*(1-x)*((1-x)**2+x**2)/(s*x)
+                        xkern(1)=(g**2/N_p)*4*vtf*(1-x)*((1-x)**2+x**2)/(s*x)
+                        xkern(2)=xkern(1)*(g_ew**2/g**2)*(qi2*vca/vtf)
                      elseif(non_limit)then
                         xfact=(1-yi)*(1-x)/x
                         prefact=4/(s*N_p)
                         call AP_reduced(m_type,i_type,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(2)=xkern(1)*(g_ew**2/g**2)*(qi2*vca/vtf)
                      endif
 c
                   elseif(ileg.eq.4)then
                      N_p=2
                      if(isspecial)N_p=1
                      if(limit)then
-                        xkern=(g**2/N_p)*( 4*vtf*(1-x)*
+                        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) )
-                        xkernazi=(g**2/N_p)*(16*vtf*s*(1-x)**2)/((s-xm12)**2)
+                        xkern(2)=xkern(1)*(g_ew**2/g**2)*(qi2*vca/vtf)
+                        xkernazi(1)=(g**2/N_p)*(16*vtf*s*(1-x)**2)/((s-xm12)**2)
+                        xkernazi(2)=xkernazi(1)*(g_ew**2/g**2)*(qi2*vca/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,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(2)=xkern(1)*(g_ew**2/g**2)*(qi2*vca/vtf)
                         call Qterms_reduced_timelike(j_type,i_type,one,z(npartner),Q)
                         Q=Q/(1-z(npartner))
-                        xkernazi=prefact*xfact*xjac(npartner)*Q/xi(npartner)
+                        xkernazi(1)=prefact*xfact*xjac(npartner)*Q/xi(npartner)
+                        xkernazi(2)=xkernazi(1)*(g_ew**2/g**2)*(qi2*vca/vtf)
                      endif
                   endif
                else
@@ -679,21 +691,26 @@
      &              (ileg.le.2 .and. abs(j_type).eq.3) )then
                if(abs(i_type).eq.3)then
 c q --> g q ( icode = 3 )
+c a --> a q
                   if(ileg.le.2)then
                      N_p=2
                      if(isspecial)N_p=1
                      if(limit)then
-                        xkern=(g**2/N_p)*4*vcf*(1-x)*((1-x)**2+1)/(s*x**2)
-                        xkernazi=-(g**2/N_p)*16*vcf*(1-x)**2/(s*x**2)
+                        xkern(1)=(g**2/N_p)*4*vcf*(1-x)*((1-x)**2+1)/(s*x**2)
+                        xkern(2)=xkern(1)*(g_ew**2/g**2)*(qi2/vcf)
+                        xkernazi(1)=-(g**2/N_p)*16*vcf*(1-x)**2/(s*x**2)
+                        xkernazi(2)=xkernazi(1)*(g_ew**2/g**2)*(qi2/vcf)
                      elseif(non_limit)then
                         xfact=(1-yi)*(1-x)/x
                         prefact=4/(s*N_p)
                         call AP_reduced(m_type,i_type,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(2)=xkern(1)*(g_ew**2/g**2)*(qi2/vcf)
                         call Qterms_reduced_spacelike(m_type,i_type,one,z(npartner),Q)
                         Q=Q/(1-z(npartner))
-                        xkernazi=prefact*xfact*xjac(npartner)*Q/xi(npartner)
+                        xkernazi(1)=prefact*xfact*xjac(npartner)*Q/xi(npartner)
+                        xkernazi(2)=xkernazi(1)*(g_ew**2/g**2)*(qi2/vcf)
                      endif
 c
                   elseif(ileg.eq.3)then
@@ -703,21 +720,24 @@
                         prefact=2/(s*N_p)
                         call AP_reduced(j_type,i_type,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(2)=xkern(1)*(g_ew**2/g**2)*(qi2/vcf)
                      endif
 c
                   elseif(ileg.eq.4)then
                      N_p=1
                      if(limit)then
-                        xkern=(g**2/N_p)*
-     &                        ( 4*vcf*(1-x)*(s**2*(1-x)**2+(s-xm12)**2) )/
-     &                        ( (s-xm12)*(s*x-xm12)**2 )
+                        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)*(g_ew**2/g**2)*(qi2/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,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(2)=xkern(1)*(g_ew**2/g**2)*(qi2/vcf)
                      endif
                   endif
                elseif(i_type.eq.8)then
@@ -726,13 +746,13 @@
                   if(ileg.le.2)then
                      N_p=1
                      if(limit)then
-                        xkern=(g**2/N_p)*4*vcf*(1+x**2)/(s*x)
+                        xkern(1)=(g**2/N_p)*4*vcf*(1+x**2)/(s*x)
                      elseif(non_limit)then
                         xfact=(1-yi)*(1-x)/x
                         prefact=4/(s*N_p)
                         call AP_reduced(m_type,i_type,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
                      endif
 c
                   elseif(ileg.eq.3)then
@@ -750,13 +770,13 @@
                            call AP_reduced_SUSY(j_type,i_type,one,z(npartner),ap)
                         endif
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
                      endif
 c
                   elseif(ileg.eq.4)then
                      N_p=1
                      if(limit)then
-                        xkern=(g**2/N_p)*4*vcf*
+                        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) )
                      elseif(non_limit)then
@@ -764,7 +784,57 @@
                         prefact=2/(s*N_p)
                         call AP_reduced(j_type,i_type,one,z(npartner),ap)
                         ap=ap/(1-z(npartner))
-                        xkern=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(1)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                     endif
+                  endif
+               elseif(i_type.eq.0)then
+c q  --> q  a ( icode = 4 )
+c sq --> sq a
+                  if(ileg.le.2)then
+                     N_p=1
+                     if(limit)then
+                        xkern(2)=(g_ew**2/N_p)*4*qj2*(1+x**2)/(s*x)
+                     else
+                        xfact=(1-yi)*(1-x)/x
+                        prefact=4/(s*N_p)
+                        call AP_reduced(m_type,i_type,one,z(npartner),ap)
+                        ap=ap/(1-z(npartner))
+                        xkern(2)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(2)=xkern(2)*(g_ew**2/g**2)*(qj2/vcf)
+                     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,one,z(npartner),ap)
+                           if(shower_mc.eq.'HERWIGPP')
+     &                     call AP_reduced_massive(j_type,i_type,one,z(npartner),
+     &                                                 xi(npartner),xm12,ap)
+                        else
+                           call AP_reduced_SUSY(j_type,i_type,one,z(npartner),ap)
+                        endif
+                        ap=ap/(1-z(npartner))
+                        xkern(2)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(2)=xkern(2)*(g_ew**2/g**2)*(qj2/vcf)
+                     endif
+c
+                  elseif(ileg.eq.4)then
+                     N_p=1
+                     if(limit)then
+                        xkern(2)=(g_ew**2/N_p)*4*qj2*
+     &                       ( s**2*(1+x**2)-2*xm12*(s*(1+x)-xm12) )/
+     &                       ( s*(s-xm12)*(s*x-xm12) )
+                     else
+                        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,one,z(npartner),ap)
+                        ap=ap/(1-z(npartner))
+                        xkern(2)=prefact*xfact*xjac(npartner)*ap/xi(npartner)
+                        xkern(2)=xkern(2)*(g_ew**2/g**2)*(qj2/vcf)
                      endif
                   endif
                else
@@ -814,13 +884,16 @@
 
         born_red=0d0
         born_red_tilde=0d0
-        xkern=xkern*gfactsf*wcc
-        xkernazi=xkernazi*gfactazi*gfactsf*wcc
+        do i=1,2
+           xkern(i)=xkern(i)*gfactsf*wcc
+           xkernazi(i)=xkernazi(i)*gfactazi*gfactsf*wcc
+        enddo
         do cflows=1,colorflow(npartner,0)
            born_red=born_red+bornbars(colorflow(npartner,cflows))
            born_red_tilde=born_red_tilde+bornbarstilde(colorflow(npartner,cflows))
         enddo
-        xmcxsec(npartner)=xkern*born_red+xkernazi*born_red_tilde
+c Change here, to include also xkern(2)!
+        xmcxsec(npartner)=xkern(1)*born_red+xkernazi(1)*born_red_tilde
         if(dampMCsubt)xmcxsec(npartner)=xmcxsec(npartner)*emscwgt(npartner)
         wgt=wgt+xmcxsec(npartner)
 
@@ -2650,6 +2723,7 @@
       zm3=(1-sqrt(1-4*xi/mdip_g**2))/2
 
 c Dead zones
+c IMPLEMENT QED DZ's!
       if(shower_mc.eq.'HERWIG6')then
          lzone=.false.
          if(ileg.le.2.and.z**2.ge.xi)lzone=.true.
@@ -2713,3 +2787,38 @@
       lzone=.false.
       return
       end
+
+
+
+      function charge(ipdg)
+c computes the electric charge given the pdg code
+      implicit none
+      integer ipdg
+      double precision charge,tmp,dipdg
+
+      dipdg=dble(ipdg)
+c quarks
+      if(abs(dipdg).eq.1) tmp=-1d0/3d0*sign(1d0,dipdg)
+      if(abs(dipdg).eq.2) tmp= 2d0/3d0*sign(1d0,dipdg)
+      if(abs(dipdg).eq.3) tmp=-1d0/3d0*sign(1d0,dipdg)
+      if(abs(dipdg).eq.4) tmp= 2d0/3d0*sign(1d0,dipdg)
+      if(abs(dipdg).eq.5) tmp=-1d0/3d0*sign(1d0,dipdg)
+      if(abs(dipdg).eq.6) tmp= 2d0/3d0*sign(1d0,dipdg)
+c leptons
+      if(abs(dipdg).eq.11)tmp=-1d0*sign(1d0,dipdg)
+      if(abs(dipdg).eq.12)tmp= 0d0
+      if(abs(dipdg).eq.13)tmp=-1d0*sign(1d0,dipdg)
+      if(abs(dipdg).eq.14)tmp= 0d0
+      if(abs(dipdg).eq.15)tmp=-1d0*sign(1d0,dipdg)
+      if(abs(dipdg).eq.16)tmp= 0d0
+c bosons
+      if(dipdg.eq.21)     tmp= 0d0
+      if(dipdg.eq.22)     tmp= 0d0
+      if(dipdg.eq.23)     tmp= 0d0
+      if(abs(dipdg).eq.24)tmp= 1d0*sign(1d0,dipdg)
+      if(dipdg.eq.25)     tmp= 0d0
+c
+      charge=tmp
+
+      return
+      end

=== modified file 'Template/NLO/SubProcesses/setcuts.f'
--- Template/NLO/SubProcesses/setcuts.f	2014-04-17 12:06:34 +0000
+++ Template/NLO/SubProcesses/setcuts.f	2014-08-26 07:56:12 +0000
@@ -409,7 +409,7 @@
                if ( itree(1,i).eq.1 .or. itree(1,i).eq.2 ) t_channel=i
                if (t_channel.ne.0) exit ! only s-channels
                mass_min(i)=mass_min(itree(1,i))+mass_min(itree(2,i))
-               if (xm(i).lt.mass_min(i)) then
+               if (xm(i).lt.mass_min(i)-vtiny) then
                   write (*,*)
      $                 'ERROR in the determination of conflicting BW',i
      $                 ,xm(i),mass_min(i)

=== modified file 'Template/NLO/SubProcesses/sumres.py'
--- Template/NLO/SubProcesses/sumres.py	2013-12-12 13:33:17 +0000
+++ Template/NLO/SubProcesses/sumres.py	2014-08-26 07:56:12 +0000
@@ -65,8 +65,9 @@
             tot+= proc['result']
             err+= math.pow(proc['error'],2)
 if dirs:
-    print "%d jobs did not terminate correctly " % len(dirs)
+    print "%d jobs did not terminate correctly: " % len(dirs)
     print '\n'.join(dirs)
+    print "The results are probably not correct. Please check the relevant log files corresponding to the above jobs for more information."
 
 processes.sort(key = lambda proc: -proc['errorABS'])
 

=== modified file 'Template/loop_material/Checks/StabilityCheckDriver.f'
--- Template/loop_material/Checks/StabilityCheckDriver.f	2012-11-09 19:55:19 +0000
+++ Template/loop_material/Checks/StabilityCheckDriver.f	2014-08-26 07:56:12 +0000
@@ -24,11 +24,23 @@
 !     
       INTEGER I,J,K
       REAL*8 P(0:3,NEXTERNAL)   ! four momenta. Energy is the zeroth component.
-      REAL*8 SQRTS,MATELEM(3),BORNELEM,AO2PI           ! sqrt(s)= center of mass energy 
+
+      INTEGER MATELEM_ARRAY_DIM
+      REAL*8 , allocatable :: MATELEM(:,:)
+
+      REAL*8 SQRTS,BORNELEM,AO2PI           ! sqrt(s)= center of mass energy 
       REAL*8 PIN(0:3), POUT(0:3)
       CHARACTER*120 BUFF(NEXTERNAL)
       CHARACTER*1 EX 
       INTEGER HELCHOICE
+      INTEGER SOCHOICE
+
+!
+!     SAVED VARIABLES
+!
+      LOGICAL INIT
+      DATA INIT/.TRUE./
+      COMMON/INITCHECKSA/INIT
 
 !     
 !     EXTERNAL
@@ -41,15 +53,22 @@
 !-----
 !     
 !---  INITIALIZATION CALLS
-!     
+!    
+      IF (INIT) THEN
+         INIT=.FALSE.
+         CALL %(proc_prefix)sGET_ANSWER_DIMENSION(MATELEM_ARRAY_DIM)
+         ALLOCATE(MATELEM(0:3,0:MATELEM_ARRAY_DIM))
+
+
 !---  Call to initialize the values of the couplings, masses and widths 
 !     used in the evaluation of the matrix element. The primary parameters of the
 !     models are read from Cards/param_card.dat. The secondary parameters are calculated
 !     in Source/MODEL/couplings.f. The values are stored in common blocks that are listed
 !     in coupl.inc .
 
-      call setpara('param_card.dat')  !first call to setup the paramaters
+          call setpara('param_card.dat')  !first call to setup the paramaters
 
+      ENDIF
       AO2PI=G**2/(8.D0*(3.14159265358979323846d0**2))
 
 
@@ -79,28 +98,35 @@
         if (MU_R.lt.0.0d0) then
           MU_R=SQRTS            
         endif
-        write(*,*) "Enter Helicity tag, -1 = summed. For loops only."
+        write(*,*) "Enter Helicity tag, -1 = summed."
         read(*,*) HELCHOICE
+        write(*,*) "Enter split_orders choice, -1 = all."
+        read(*,*) SOCHOICE
+        IF (SOCHOICE.NE.-1) THEN
+          CALL %(proc_prefix)sSET_COUPLINGORDERS_TARGET(SOCHOICE)
+        ENDIF
 !---  Update the couplings with the new MU_R
         CALL UPDATE_AS_PARAM()
 !     
 !     Now we can call the matrix element!
 !
-        CALL SMATRIX(P,BORNELEM) 
         IF (HELCHOICE.EQ.-1) THEN
-          CALL SLOOPMATRIX(P,MATELEM)
+          CALL %(proc_prefix)sSLOOPMATRIX(P,MATELEM)
         ELSE
-          CALL SLOOPMATRIXHEL(P,HELCHOICE,MATELEM)
+          CALL %(proc_prefix)sSLOOPMATRIXHEL(P,HELCHOICE,MATELEM)
         ENDIF
         write(*,*) '##TAG#RESULT_START#TAG##'
         do i=1,nexternal      
           write (*,'(a2,1x,5e25.15)') 'PS',P(0,i),P(1,i),P(2,i),P(3,i)
         enddo
         write (*,'(a3,1x,i2)') 'EXP',-(2*nexternal-8)
-        write (*,'(a4,1x,1e25.15)') 'BORN',BORNELEM
-        write (*,'(a3,1x,1e25.15)') 'FIN',MATELEM(1)/BORNELEM/AO2PI
-        write (*,'(a4,1x,1e25.15)') '1EPS',MATELEM(2)/BORNELEM/AO2PI
-        write (*,'(a4,1x,1e25.15)') '2EPS',MATELEM(3)/BORNELEM/AO2PI
+        write (*,'(a4,1x,1e25.15)') 'BORN',MATELEM(0,0)
+        write (*,'(a3,1x,1e25.15)') 'FIN',
+     &MATELEM(1,0)/MATELEM(0,0)/AO2PI
+        write (*,'(a4,1x,1e25.15)') '1EPS',
+     &MATELEM(2,0)/MATELEM(0,0)/AO2PI
+        write (*,'(a4,1x,1e25.15)') '2EPS',
+     &MATELEM(3,0)/MATELEM(0,0)/AO2PI
         write (*,*) 'Export_Format Default'
         write(*,*) '##TAG#RESULT_STOP#TAG##'      
       enddo

=== modified file 'Template/loop_material/Checks/StabilityCheckDriver_loop_induced.f'
--- Template/loop_material/Checks/StabilityCheckDriver_loop_induced.f	2012-11-09 19:55:19 +0000
+++ Template/loop_material/Checks/StabilityCheckDriver_loop_induced.f	2014-08-26 07:56:12 +0000
@@ -24,11 +24,22 @@
 !     
       INTEGER I,J,K
       REAL*8 P(0:3,NEXTERNAL)   ! four momenta. Energy is the zeroth component.
-      REAL*8 SQRTS,MATELEM(3),BORNELEM,AO2PI           ! sqrt(s)= center of mass energy 
+      
+      INTEGER MATELEM_ARRAY_DIM
+      REAL*8 , allocatable :: MATELEM(:,:)
+
+      REAL*8 SQRTS,BORNELEM,AO2PI           ! sqrt(s)= center of mass energy 
       REAL*8 PIN(0:3), POUT(0:3)
       CHARACTER*120 BUFF(NEXTERNAL)
       CHARACTER*1 EX 
-      INTEGER HELCHOICE
+      INTEGER HELCHOICE, SOCHOICE
+
+!
+!     SAVED VARIABLES
+!
+      LOGICAL INIT
+      DATA INIT/.TRUE./
+      COMMON/INITCHECKSA/INIT
 
 !     
 !     EXTERNAL
@@ -41,14 +52,22 @@
 !-----
 !     
 !---  INITIALIZATION CALLS
-!     
+!    
+
+      IF (INIT) THEN
+         INIT=.FALSE.
+         CALL %(proc_prefix)sGET_ANSWER_DIMENSION(MATELEM_ARRAY_DIM)
+         ALLOCATE(MATELEM(0:3,0:MATELEM_ARRAY_DIM))
+
 !---  Call to initialize the values of the couplings, masses and widths 
 !     used in the evaluation of the matrix element. The primary parameters of the
 !     models are read from Cards/param_card.dat. The secondary parameters are calculated
 !     in Source/MODEL/couplings.f. The values are stored in common blocks that are listed
 !     in coupl.inc .
 
-      call setpara('param_card.dat')  !first call to setup the paramaters
+         call setpara('param_card.dat')  !first call to setup the paramaters
+
+      ENDIF
 
       AO2PI=G**2/(8.D0*(3.14159265358979323846d0**2))
 
@@ -81,15 +100,20 @@
         endif
         write(*,*) "Enter Helicity tag, -1 = summed. For loops only."
         read(*,*) HELCHOICE
+        write(*,*) "Enter split_orders choice, -1 = all."
+        read(*,*) SOCHOICE
+        IF (SOCHOICE.NE.-1) THEN
+          CALL %(proc_prefix)sSET_COUPLINGORDERS_TARGET(SOCHOICE)
+        ENDIF
 !---  Update the couplings with the new MU_R
         CALL UPDATE_AS_PARAM()
 !     
 !     Now we can call the matrix element!
 !
         IF (HELCHOICE.EQ.-1) THEN
-          CALL SLOOPMATRIX(P,MATELEM)
+          CALL %(proc_prefix)sSLOOPMATRIX(P,MATELEM)
         ELSE
-          CALL SLOOPMATRIXHEL(P,HELCHOICE,MATELEM)
+          CALL %(proc_prefix)sSLOOPMATRIXHEL(P,HELCHOICE,MATELEM)
         ENDIF
         write(*,*) '##TAG#RESULT_START#TAG##'
         do i=1,nexternal      
@@ -97,7 +121,7 @@
         enddo
         write (*,'(a3,1x,i2)') 'EXP',-(2*nexternal-8)
         write (*,'(a4,1x,1e25.15)') 'BORN',0.0d0
-        write (*,'(a3,1x,1e25.15)') 'FIN',MATELEM(1)
+        write (*,'(a3,1x,1e25.15)') 'FIN',MATELEM(1,0)
         write (*,'(a4,1x,1e25.15)') '1EPS',0.0d0
         write (*,'(a4,1x,1e25.15)') '2EPS',0.0d0
         write (*,*) 'Export_Format LoopInduced'

=== added directory 'Template/loop_material/StandAlone/Cards'
=== removed directory 'Template/loop_material/StandAlone/Cards'
=== renamed file 'Template/loop_material/StandAlone/Cards/MadLoopParams.dat' => 'Template/loop_material/StandAlone/Cards/MadLoopParams.dat'
--- Template/loop_material/StandAlone/Cards/MadLoopParams.dat	2013-05-01 16:26:01 +0000
+++ Template/loop_material/StandAlone/Cards/MadLoopParams.dat	2014-08-26 07:56:12 +0000
@@ -2,6 +2,37 @@
 ! The name of the variable to define must start with the '#' sign and then
 ! the value should be put immediately on the next line.
 
+!     
+#MLReductionLib
+1|4|3|2
+! Default :: 1|4|3|2
+!     The tensor integral reduction library.The current choices are:
+!     1 | CutTools
+!     2 | PJFry++
+!     3 | IREGI
+!     4 | Golem95
+!     One can use the combinations to reduce integral,e.g.
+!     1|2|3 means first use CutTools, if it is not stable, use PJFry++, 
+!     if it is still unstable, use IREGI. If it failed, use QP of CutTools.
+!  
+
+!     =================================================================================
+!     The parameters below set the parameters for IREGI
+!     =================================================================================
+
+#IREGIMODE
+2
+! Default :: 2
+!     IMODE=0, IBP reduction
+!     IMODE=1, PaVe reduction
+!     IMODE=2, PaVe reduction with stablility improved by IBP reduction
+
+#IREGIRECY
+.TRUE.
+! Default :: .TRUE.
+!     Use RECYCLING OR NOT IN IREGI
+!
+
 !     =================================================================================
 !     The parameters below set the stability checks of MadLoop at run time
 !     =================================================================================
@@ -78,9 +109,10 @@
 !     The parameters below set two CutTools internal parameters accessible to the user. 
 !     =================================================================================
 
-!     Choose here what library to chose for CutTools to compute the scalar loops of the
+!     Choose here what library to chose for CutTools/TIR to compute the scalar loops of the
 !     master integral basis. The choices are as follows:
-!     2 | AVH
+!     (Does not apply for Golem95, where OneLOop is always used)
+!     2 | OneLOop
 !     3 | QCDLoop
 #CTLoopLibrary
 2
@@ -131,6 +163,13 @@
 1.0d-9
 ! Default :: 1.0d-9
 
+!     Setting the on-shell threshold for deciding whether the invariant variables
+! of external momenta are on-shell or not. It will only be used in constructing
+! s-matrix in Golem95.
+#OSThres
+1.0d-8
+! Default :: 1.0d-8
+
 !     The setting below is recommended to be on as it allows to systematically used the 
 !     first PS point thrown at ML5 to be used for making sure that the helicity filter
 !     read from HelFilter.dat is consistent as it might be no longer up to date with
@@ -139,6 +178,14 @@
 .TRUE.
 ! Default :: .TRUE.
 
+!     This decides whether to write out the helicity and loop filters to the files 
+!     HelFilters.dat and LoopFilters.dat to save them for future runs. It usually
+!     preferable but sometimes not desired because of the need of threadlocks in the
+!     context of mpi parallelization. So it can be turned off here in such cases.
+#WriteOutFilters
+.TRUE.
+! Default :: .TRUE.
+
 !     Some loop contributions may be zero for some helicities which are however
 !     contributing. In order to save their computing time, you can chose here to try
 !     to filter them out. The gain is typically minimal, so it is turned off by default.
@@ -146,7 +193,7 @@
 .FALSE.
 ! Default :: .FALSE.
 
-!     This decides wether consecutive consistency for the loop filtering setup is also 
+!     This decides whether consecutive consistency for the loop filtering setup is also 
 !     required.
 #LoopInitStartOver
 .FALSE.

=== added directory 'Template/loop_material/StandAlone/SubProcesses/MadLoop5_resources'
=== added file 'Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.f'
--- Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.f	1970-01-01 00:00:00 +0000
+++ Template/loop_material/StandAlone/SubProcesses/MadLoopCommons.f	2014-08-26 07:56:12 +0000
@@ -0,0 +1,341 @@
+      subroutine joinPath(str1,str2,path)
+
+      character*(*) str1
+      character*(*) str2
+      character*(*) path
+
+      integer i,j,k
+
+      i =1
+      do while (i.le.LEN(str1))
+      if(str1(i:i).eq.' ') goto 800     
+      path(i:i) = str1(i:i)
+      i=i+1
+      enddo
+800   continue
+      j=1
+      do while (j.le.LEN(str2))
+      if(str2(j:j).eq.' ') goto 801
+      path(i-1+j:i-1+j) = str2(j:j)
+      j=j+1
+      enddo
+801   continue
+      k=i+j-1
+      do while (k.le.LEN(path))
+      path(k:k) = ' '
+      k=k+1
+      enddo
+
+      return
+
+      end
+
+      subroutine setMadLoopPath(path)
+
+      character(512) path
+      character(512) dummy      
+
+      character(512) prefix,fpath
+      character(17) nameToCheck
+      parameter (nameToCheck='MadLoopParams.dat')
+
+      LOGICAL ML_INIT
+      DATA ML_INIT/.TRUE./
+      common/ML_INIT/ML_INIT
+
+      LOGICAL CTINIT,TIRINIT,GOLEMINIT
+      DATA CTINIT,TIRINIT,GOLEMINIT/.TRUE.,.TRUE.,.TRUE./
+      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT, GOLEMINIT
+
+      character(512) MLPath
+      data MLPath/'[[NA]]'/      
+      common/MLPATH/MLPath
+
+      integer i
+
+C     Just a dummy call for LD to pick up this function
+C     when creating the BLHA2 dynamic library
+      dummy = ' '
+      CALL SETPARA2(dummy)
+
+      if (LEN(path).ge.4 .and. path(1:4).eq.'auto') then
+          if (MLPath(1:6).eq.'[[NA]]') then
+C     Try to automatically find the path
+          prefix='./'
+          call joinPath(prefix,nameToCheck,fpath)
+          OPEN(1, FILE=fpath, ERR=1, STATUS='OLD',      
+     $    ACTION='READ')
+          MLPath=prefix
+          goto 10
+1         continue
+          close(1)
+          prefix='./MadLoop5_resources/'
+          call joinPath(prefix,nameToCheck,fpath)
+          OPEN(1, FILE=fpath, ERR=2, STATUS='OLD',      
+     $    ACTION='READ')
+          MLPath=prefix
+          goto 10
+2         continue
+          close(1)
+          prefix='../MadLoop5_resources/'
+          call joinPath(prefix,nameToCheck,fpath)
+          OPEN(1, FILE=fpath, ERR=66, STATUS='OLD',      
+     $    ACTION='READ')
+          MLPath=prefix
+          goto 10
+66        continue
+          close(1)
+c     We could not automatically find the auxiliary files
+          write(*,*) '==='
+          write(*,*) 'ERROR: MadLoop5 could not automatically find',
+     $    ' the file MadLoopParams.dat.'
+          write(*,*) '==='
+          write(*,*) '(Try using <CALL setMadLoopPath(/my/path)>',
+     $    ' (before your first call to MadLoop) in order to',
+     $    ' set the directory where this file is located as well as',
+     $    ' other auxiliary files, such as <xxx>_ColorNumFactors.dat',
+     $    ', <xxx>_ColorDenomFactors.dat, etc..)'
+          stop
+10        continue
+          close(1)
+          return
+          endif
+      else
+C     Use the one specified by the user
+C     Make sure there is a separator added
+      i =1
+      do while (i.le.LEN(path) .and. path(i:i).ne.' ')
+      i=i+1
+      enddo
+      if (path(i-1:i-1).ne.'/') then
+          path(i:i) = '/'
+      endif
+      MLpath=path          
+      endif
+
+C     Check that the FilePath set is correct
+      call joinPath(MLPath,nameToCheck,fpath)
+      OPEN(1, FILE=fpath, ERR=3, STATUS='OLD',      
+     $ACTION='READ')
+      goto 11
+3     continue
+      close(1)
+      write(*,*) '==='
+      write(*,*) 'ERROR: The MadLoop5 auxiliary files could not',
+     $' be found in ',MLPath
+      write(*,*) '==='
+      stop
+11    continue
+      close(1)
+
+      end
+
+      INTEGER FUNCTION SET_RET_CODE_U(MLRed,DOING_QP,STABLE)
+C
+C This functions returns the value of U
+C
+C
+C     U == 0
+C         Not stable.
+C     U == 1
+C         Stable with CutTools in double precision.
+C     U == 2
+C         Stable with PJFry++.
+C     U == 3
+C         Stable with IREGI.
+C     U == 4
+C         Stable with Golem95
+C     U == 9
+C         Stable with CutTools in quadruple precision.
+C
+      IMPLICIT NONE
+C     
+C CONSTANTS
+C
+C
+C ARGUMENTS
+C
+      INTEGER MLRed
+      LOGICAL DOING_QP,STABLE
+C
+C LOCAL VARIABLES
+C
+C
+C FUNCTION
+C
+C
+C BEGIN CODE
+C
+      IF(.NOT.STABLE)THEN
+         SET_RET_CODE_U=0
+         RETURN
+      ENDIF
+      IF(DOING_QP)THEN
+         IF(MLRed.EQ.1)THEN
+            SET_RET_CODE_U=9
+            RETURN
+         ELSE
+            STOP 'Only CutTools can use quardruple precision'
+         ENDIF
+      ENDIF
+      IF(MLRed.GE.1.AND.MLRed.LE.4)THEN
+         SET_RET_CODE_U=MLRed
+      ELSE
+         STOP 'Only CutTools,PJFry++,IREGI,Golem95 are available'
+      ENDIF
+      END
+
+      SUBROUTINE DETECT_LOOPLIB(LIBNUM,NLOOPLINE,RANK,complex_mass,
+     $     LPASS)
+C
+C DETECT WHICH LOOP LIB PASSED
+C
+      IMPLICIT NONE
+C
+C CONSTANTS
+C
+C
+C ARGUMENTS
+C
+      INTEGER LIBNUM,NLOOPLINE,RANK
+      LOGICAL complex_mass,LPASS
+C     
+C LOCAL VARIABLES
+C
+C
+C GLOBAL VARIABLES
+C
+C ----------
+C BEGIN CODE
+C ----------
+      IF(LIBNUM.EQ.1)THEN
+C CutTools
+         CALL DETECT_CUTTOOLS(NLOOPLINE,RANK,complex_mass,LPASS)
+      ELSEIF(LIBNUM.EQ.2)THEN
+C PJFry++
+         CALL DETECT_PJFRY(NLOOPLINE,RANK,complex_mass,LPASS)
+      ELSEIF(LIBNUM.EQ.3)THEN
+C IREGI
+         CALL DETECT_IREGI(NLOOPLINE,RANK,complex_mass,LPASS)
+      ELSEIF(LIBNUM.EQ.4)THEN
+C Golem95
+         CALL DETECT_GOLEM(NLOOPLINE,RANK,complex_mass,LPASS)
+      ELSE
+         STOP 'ONLY CUTTOOLS,PJFry++,IREGI,Golem95 are provided'
+      ENDIF
+      RETURN
+      END
+
+      SUBROUTINE DETECT_CUTTOOLS(NLOOPLINE,RANK,complex_mass,LPASS)
+C
+C DETECT THE CUTTOOLS CAN BE USED OR NOT
+C
+      IMPLICIT NONE
+C
+C CONSTANTS
+C
+C
+C ARGUMENTS
+C
+      INTEGER NLOOPLINE,RANK
+      LOGICAL complex_mass,LPASS
+C
+C LOCAL VARIABLES
+C
+C
+C GLOBAL VARIABLES
+C
+C ----------
+C BEGIN CODE
+C ----------
+      LPASS=.TRUE.
+      IF(NLOOPLINE+1.LT.RANK)LPASS=.FALSE.
+      RETURN
+      END
+
+      SUBROUTINE DETECT_PJFRY(NLOOPLINE,RANK,complex_mass,LPASS)
+C
+C DETECT THE PJFRY++ CAN BE USED OR NOT
+C
+      IMPLICIT NONE
+C
+C CONSTANTS
+C
+C
+C ARGUMENTS
+C
+      INTEGER NLOOPLINE,RANK
+      LOGICAL complex_mass,LPASS
+C
+C LOCAL VARIABLES
+C
+C
+C GLOBAL VARIABLES
+C
+C ----------
+C BEGIN CODE
+C ----------
+      LPASS=.TRUE.
+      IF(NLOOPLINE.LT.RANK.OR.RANK.GT.5
+     $.OR.NLOOPLINE.GT.5.OR.complex_mass.OR.NLOOPLINE.eq.1) THEN
+        LPASS=.FALSE.
+      ENDIF
+      RETURN
+      END
+
+      SUBROUTINE DETECT_IREGI(NLOOPLINE,RANK,complex_mass,LPASS)
+C     
+C DETECT THE IREGI CAN BE USED OR NOT
+C
+      IMPLICIT NONE
+C
+C CONSTANTS
+C
+C
+C ARGUMENTS
+C
+      INTEGER NLOOPLINE,RANK
+      LOGICAL complex_mass,LPASS
+C     
+C LOCAL VARIABLES
+C
+C
+C GLOBAL VARIABLES
+C
+C ----------
+C BEGIN CODE
+C ----------
+      LPASS=.TRUE.
+      IF(NLOOPLINE.GE.7.OR.RANK.GE.7)LPASS=.FALSE.
+      RETURN
+      END
+
+      SUBROUTINE DETECT_GOLEM(NLOOPLINE,RANK,complex_mass,LPASS)
+C
+C DETECT THE Golem95 CAN BE USED OR NOT
+C
+            IMPLICIT NONE
+C
+C CONSTANTS
+C
+C
+C ARGUMENTS
+C
+            INTEGER NLOOPLINE,RANK
+            LOGICAL complex_mass,LPASS
+C
+C LOCAL VARIABLES
+C
+C
+C GLOBAL VARIABLES
+C
+C ----------
+C BEGIN CODE
+C ----------
+
+      LPASS=.TRUE.
+      IF(NLOOPLINE.GE.7.OR.RANK.GE.7.OR.NLOOPLINE.LE.1)LPASS=.FALSE.
+      IF(NLOOPLINE.LE.5.AND.RANK.GT.NLOOPLINE+1)LPASS=.FALSE.
+      IF(NLOOPLINE.EQ.6.AND.RANK.GT.NLOOPLINE)LPASS=.FALSE.
+      RETURN
+      END

=== modified file 'Template/loop_material/StandAlone/SubProcesses/MadLoopParamReader.f'
--- Template/loop_material/StandAlone/SubProcesses/MadLoopParamReader.f	2013-05-01 16:26:01 +0000
+++ Template/loop_material/StandAlone/SubProcesses/MadLoopParamReader.f	2014-08-26 07:56:12 +0000
@@ -2,10 +2,14 @@
 
       implicit none
 
-      CHARACTER*64 fileName, buff, buff2, mode
+      CHARACTER(512) fileName, buff, buff2, mode
+      CHARACTER*10 MLReductionLib_str,MLReductionLib_str_save
+      CHARACTER*2  MLReductionLib_char
+      INTEGER MLRed,i,j,k
+
       include "MadLoopParams.inc"
 
-      logical printParam, couldRead, paramPrinted
+      logical printParam, couldRead, paramPrinted, find
       data paramPrinted/.FALSE./
       couldRead=.False.
 !     Default parameters
@@ -68,6 +72,12 @@
                  stop 'ZeroThres must be > 0'
              endif
 
+           else if (buff .eq. '#OSThres') then
+              read(666,*,end=999) OSThres
+              if (OSThres.le.0.0d0) then
+                 stop 'OSThres must be > 0'
+              endif
+
            else if (buff .eq. '#CheckCycle') then
              read(666,*,end=999) CheckCycle
              if (CheckCycle.lt.1) then
@@ -92,13 +102,65 @@
           else if (buff .eq. '#HelInitStartOver') then
              read(666,*,end=999) HelInitStartOver
 
+          else if (buff .eq. '#WriteOutFilters') then
+             read(666,*,end=999) WriteOutFilters
+
           else if (buff .eq. '#ImprovePSPoint') then
              read(666,*,end=999) ImprovePSPoint
              if (ImprovePSPoint .lt. -1 .or.
      &           ImprovePSPoint .gt. 2 ) then
                stop 'ImprovePSPoint must be >= -1 and <=2.'
              endif
-
+          else if (buff .eq. '#MLReductionLib') then
+             read(666,*,end=999) MLReductionLib_str
+             MLReductionLib(1:5)=0
+             MLReductionLib_str_save=MLReductionLib_str
+             j=0
+             DO
+                i=index(MLReductionLib_str,'|')
+                IF(i.EQ.0)THEN
+                   MLReductionLib_char=MLReductionLib_str
+                ELSE
+                   MLReductionLib_char=MLReductionLib_str(:i-1)
+                ENDIF
+                IF(MLReductionLib_char.EQ.'1 ')THEN
+                   MLRed=1
+                ELSEIF(MLReductionLib_char.EQ.'2 ')THEN
+                   MLRed=2
+                ELSEIF(MLReductionLib_char.EQ.'3 ')THEN
+                   MLRed=3
+                ELSEIF(MLReductionLib_char.EQ.'4 ')THEN
+                   MLRed=4
+                ELSE
+                   PRINT *, 'MLReductionLib is wrong: '//
+     $                  TRIM(MLReductionLib_str_save)
+                   STOP
+                ENDIF
+                find=.FALSE.
+                DO k=1,j
+                   IF(MLReductionLib(k).EQ.MLRed)THEN
+                      find=.TRUE.
+                      EXIT
+                   ENDIF
+                ENDDO
+                IF(.NOT.find)THEN
+                   j=j+1
+                   MLReductionLib(j)=MLRed
+                ENDIF
+                IF(i.EQ.0)THEN
+                   EXIT
+                ELSE
+                   MLReductionLib_str=MLReductionLib_str(i+1:)
+                ENDIF
+             ENDDO
+          else if (buff .eq. '#IREGIRECY') then
+             read(666,*,end=999) IREGIRECY
+          else if (buff .eq. '#IREGIMODE') then
+             read(666,*,end=999) IREGIMODE
+             if (IREGIMODE .lt. 0 .or.
+     &            IREGIMODE .gt.2) then
+                stop 'IREGIMODE must be >=0 and <=2.'
+             endif
           else
              write(*,*) 'The parameter name ',buff(2:),
      &' is not reckognized.'
@@ -136,6 +198,10 @@
       endif
       write(*,*)
      & '==============================================================='
+      write(*,*) ' > MLReductionLib            = '
+     $     //TRIM(MLReductionLib_str_save)
+      write(*,*) ' > IREGIMODE                 = ',IREGIMODE
+      write(*,*) ' > IREGIRECY                 = ',IREGIRECY
       write(*,*) ' > CTModeRun                 = ',CTModeRun
       write(*,*) ' > MLStabThres               = ',MLStabThres
       write(*,*) ' > NRotations_DP             = ',NRotations_DP
@@ -152,6 +218,8 @@
       write(*,*) ' > LoopInitStartOver         = ',LoopInitStartOver
       write(*,*) ' > HelInitStartOver          = ',HelInitStartOver
       write(*,*) ' > ZeroThres                 = ',ZeroThres
+      write(*,*) ' > OSThres                   = ',OSThres
+      write(*,*) ' > WriteOutFilters           = ',WriteOutFilters      
       write(*,*)
      & '==============================================================='
       paramPrinted=.TRUE.
@@ -167,6 +235,10 @@
       
       include "MadLoopParams.inc"
 
+      MLReductionLib(1)=1
+      MLReductionLib(2:5)=0
+      IREGIMODE=2
+      IREGIRECY=.TRUE.
       CTModeInit=0
       CTModeRun=-1
       NRotations_DP=1
@@ -180,7 +252,9 @@
       DoubleCheckHelicityFilter=.True.
       LoopInitStartOver=.False.
       HelInitStartOver=.False.
+      WriteOutFilters=.True.
       ZeroThres=1.0d-9
+      OSThres=1.0d-13
       ImprovePSPoint=2
 
       end

=== modified file 'Template/loop_material/StandAlone/SubProcesses/MadLoopParams.inc'
--- Template/loop_material/StandAlone/SubProcesses/MadLoopParams.inc	2013-05-01 16:26:01 +0000
+++ Template/loop_material/StandAlone/SubProcesses/MadLoopParams.inc	2014-08-26 07:56:12 +0000
@@ -6,14 +6,16 @@
 !====================================================================
 !
 	  integer CTModeInit,CTModeRun,CheckCycle,MaxAttempts,
-     &CTLoopLibrary,NRotations_DP,NRotations_QP,ImprovePSPoint
+     &CTLoopLibrary,NRotations_DP,NRotations_QP,ImprovePSPoint,
+     &MLReductionLib(5),IREGIMODE
 
-      real*8 MLStabThres,CTStabThres,ZeroThres
+      real*8 MLStabThres,CTStabThres,ZeroThres,OSThres
 
       logical UseLoopFilter,LoopInitStartOver,DoubleCheckHelicityFilter
-	  logical HelInitStartOver
+	  logical HelInitStartOver,IREGIRECY,WriteOutFilters
 
       common /MADLOOP/CTModeInit,CTModeRun,NRotations_DP,NRotations_QP,
      &ImprovePSPoint,CheckCycle, MaxAttempts,UseLoopFilter,MLStabThres,
      &CTStabThres,CTLoopLibrary,LoopInitStartOver,
-     &DoubleCheckHelicityFilter,ZeroThres,HelInitStartOver
+     &DoubleCheckHelicityFilter,ZeroThres,OSThres,HelInitStartOver,
+     &MLReductionLib,IREGIMODE,IREGIRECY,WriteOutFilters

=== removed file 'Template/loop_material/StandAlone/SubProcesses/check_sa.f'
--- Template/loop_material/StandAlone/SubProcesses/check_sa.f	2013-12-13 07:04:16 +0000
+++ Template/loop_material/StandAlone/SubProcesses/check_sa.f	1970-01-01 00:00:00 +0000
@@ -1,583 +0,0 @@
-      PROGRAM DRIVER
-!**************************************************************************
-!     THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT.
-!     IT USES A SIMPLE PHASE SPACE GENERATOR
-!     Fabio Maltoni - 3rd Febraury 2007
-!**************************************************************************
-      IMPLICIT NONE
-!     
-!     CONSTANTS  
-!     
-      REAL*8 ZERO
-      PARAMETER (ZERO=0D0)
-
-      LOGICAL READPS
-      PARAMETER (READPS = .FALSE.)
-
-      INTEGER NPSPOINTS
-      PARAMETER (NPSPOINTS = 4)
-
-!     
-!     INCLUDE FILES
-!     
-!---  the include file with the values of the parameters and masses      
-      INCLUDE "coupl.inc"
-!---  integer nexternal ! number particles (incoming+outgoing) in the me 
-      INCLUDE "nexternal.inc" 
-!---  particle masses
-      REAL*8 PMASS(NEXTERNAL)      
-!---  integer    n_max_cg
-      INCLUDE "ngraphs.inc"     !how many diagrams (could be useful to know...)
-
-!     
-!     LOCAL
-!     
-      INTEGER I,J,K
-      REAL*8 P(0:3,NEXTERNAL)   ! four momenta. Energy is the zeroth component.
-      REAL*8 SQRTS,MATELEM(3),BORNELEM,AO2PI           ! sqrt(s)= center of mass energy 
-      REAL*8 PIN(0:3), POUT(0:3)
-      CHARACTER*120 BUFF(NEXTERNAL)
-      REAL*8 PREC_FOUND
-      INTEGER RETURNCODE, HUNDREDS, TENS, UNITS
-!     
-!     EXTERNAL
-!     
-      REAL*8 DOT
-      EXTERNAL DOT
-      
-!-----
-!     BEGIN CODE
-!-----
-!     
-!---  INITIALIZATION CALLS
-!     
-!---  Call to initialize the values of the couplings, masses and widths 
-!     used in the evaluation of the matrix element. The primary parameters of the
-!     models are read from Cards/param_card.dat. The secondary parameters are calculated
-!     in Source/MODEL/couplings.f. The values are stored in common blocks that are listed
-!     in coupl.inc .
-
-      call setpara('param_card.dat')  !first call to setup the paramaters
-      include "pmass.inc"             !set up masses
-
-      AO2PI=G**2/(8.D0*(3.14159265358979323846d0**2))
-
-      write(*,*) 'AO2PI=',AO2PI
-!---  Now use a simple multipurpose PS generator (RAMBO) just to get a 
-!     RANDOM set of four momenta of given masses pmass(i) to be used to evaluate 
-!     the madgraph matrix-element.       
-!     Alternatevely, here the user can call or set the four momenta at his will, see below.
-!           
-      IF(nincoming.EQ.1) THEN
-         SQRTS=PMASS(1)
-      ELSE
-         SQRTS=1000d0              !CMS energy in GEV
-      ENDIF
-
-      call printout()
-
-      do K=1,NPSPOINTS
-
-      if(READPS) then
-        open(967, file="PS.input", err=976, status='OLD', action='READ')
-        do i=1,NExternal
-          read(967,*,end=978) P(0,i),P(1,i),P(2,i),P(3,i)
-        enddo
-       goto 978
-  976  continue
-         stop 'Could not read the PS.input phase-space point.'
-  978  continue
-       close(967)
-      else
-        if ((nincoming.eq.2).and.((nexternal - nincoming .eq.1))) then
-          if (pmass(3).eq.0.0d0) then
-              stop 'Cannot generate 2>1 kin. config. with m3=0.0d0'
-          else
-              ! deal with the case of only one particle in the final
-              ! state
-               p(0,1) = pmass(3)/2d0
-               p(1,1) = 0d0
-               p(2,1) = 0d0
-               p(3,1) = pmass(3)/2d0
-               if (pmass(1).GT.0d0) then
-                 p(3,1) = dsqrt(pmass(3)**2/4d0 - pmass(1)**2)
-               endif
-               p(0,2) = pmass(3)/2d0
-               p(1,2) = 0d0
-               p(2,2) = 0d0
-               p(3,2) = -pmass(3)/2d0
-               if (pmass(2) > 0d0) then
-                 p(3,2) = -dsqrt(pmass(3)**2/4d0 - pmass(1)**2)
-               endif
-               p(0,3) = pmass(3)
-               p(1,3) = 0d0
-               p(2,3) = 0d0
-               p(3,3) = 0d0
-          endif
-        else
-          CALL GET_MOMENTA(SQRTS,PMASS,P)
-        endif
-      endif
-
-      do i=0,3
-        PIN(i)=0.0d0
-        do j=1,nincoming
-          PIN(i)=PIN(i)+p(i,j)
-        enddo
-      enddo
-
-!---  In standalone mode, always use sqrt_s as the renormalization
-!     scale.
-      SQRTS=dsqrt(dabs(DOT(PIN(0),PIN(0))))
-      MU_R=SQRTS
-
-!---  Update the couplings with the new MU_R
-      CALL UPDATE_AS_PARAM()
-
-!-- Possible tweaks to the run
-!     CALL FORCE_STABILITY_CHECK(.TRUE.)
-
-!     
-!     Now we can call the matrix element!
-!
-      CALL SMATRIX(P,BORNELEM)
-      CALL SLOOPMATRIX_THRES(P,MATELEM,
-     &-1.0d0,PREC_FOUND,RETURNCODE)
-!
-!        write the information on the four momenta 
-!
-      if (K.eq.NPSPOINTS) then
-      write (*,*)
-      write (*,*) " Phase space point:"
-      write (*,*) "--------------------------------------------------",
-     &"---------------------------"
-      write (*,*)  "n        E             px             py         ",
-     &"pz               m "
-      do i=1,nexternal
-         write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i), 
-     .dsqrt(dabs(DOT(p(0,i),p(0,i))))
-      enddo
-      write (*,*) "--------------------------------------------------",
-     &"---------------------------"
-      UNITS=MOD(RETURNCODE,10)
-      TENS=(MOD(RETURNCODE,100)-UNITS)/10
-      HUNDREDS=(RETURNCODE-TENS*10-UNITS)/100
-      if (HUNDREDS.eq.1) then
-        if (TENS.eq.3.or.TENS.eq.4) then
-          write(*,*) 'Unknown numerical stability because ',
-     &' MadLoop is in the initialization stage.'
-        else
-          write(*,*) 'Unknown numerical stability, check ',
-     &'CTModeRun value in MadLoopParams.dat.'
-        endif
-      elseif (HUNDREDS.eq.2) then
-        write(*,*) 'Stable kinematic configuration (SPS).'
-      elseif (HUNDREDS.eq.3) then
-        write(*,*) 'Unstable kinematic configuration (UPS).'
-        write(*,*) 'Quadruple precision rescue successful.'
-      elseif (HUNDREDS.eq.4) then
-        write(*,*) 'Exceptional kinematic configuration (EPS).'
-        write(*,*) 'Both double an quadruple precision computations',
-     %' are unstable.'
-      endif
-      if (TENS.eq.2.or.TENS.eq.4) then
-        write(*,*) 'Quadruple precision computation used.'
-      endif
-      if (HUNDREDS.ne.1) then
-        if (PREC_FOUND.gt.0.0d0) then
-          write(*,'(1x,a23,1x,1e10.2)') 'Relative accuracy     =',
-     &PREC_FOUND
-        elseif (PREC_FOUND.eq.0.0d0) then
-          write(*,'(1x,a23,1x,1e10.2,1x,a30)') 'Relative accuracy     ='
-     &,PREC_FOUND,'(i.e. beyond double precision)'
-        else
-          write(*,*) 'Estimated accuracy could not be',
-     &' computed for an unknown reason.'
-        endif
-      endif
-      write (*,*) "---------------------------------------------------",
-     &"--------------------------"
-      write (*,*) "Matrix element born   = ", BORNELEM, " GeV^",
-     &-(2*nexternal-8)
-      write (*,*) "Matrix element finite = ", MATELEM(1), " GeV^",
-     &-(2*nexternal-8)
-      write (*,*) "Matrix element 1eps   = ", MATELEM(2), " GeV^",
-     &-(2*nexternal-8)
-      write (*,*) "Matrix element 2eps   = ", MATELEM(3), " GeV^",
-     &-(2*nexternal-8)
-      write (*,*) "---------------------------------------------------",
-     &"--------------------------"
-      write (*,*) "finite / (born*ao2pi) = ", MATELEM(1)/BORNELEM/AO2PI
-      write (*,*) "1eps   / (born*ao2pi) = ", MATELEM(2)/BORNELEM/AO2PI
-      write (*,*) "2eps   / (born*ao2pi) = ", MATELEM(3)/BORNELEM/AO2PI
-      write (*,*) "---------------------------------------------------",
-     &"--------------------------"
-
-      open(69, file="result.dat", err=976, action='WRITE')
-      do i=1,nexternal      
-          write (69,'(a2,1x,5e25.15)') 'PS',P(0,i),P(1,i),P(2,i),P(3,i)
-      enddo
-      write (69,'(a3,1x,i2)') 'EXP',-(2*nexternal-8)
-      write (69,'(a4,1x,1e25.15)') 'BORN',BORNELEM
-      write (69,'(a3,1x,1e25.15)') 'FIN',MATELEM(1)/BORNELEM/AO2PI
-      write (69,'(a4,1x,1e25.15)') '1EPS',MATELEM(2)/BORNELEM/AO2PI
-      write (69,'(a4,1x,1e25.15)') '2EPS',MATELEM(3)/BORNELEM/AO2PI
-      write (69,'(a3,1x,1e10.2)')  'ACC',PREC_FOUND
-      write (69,'(a7,1x,i3)')     'RETCODE',RETURNCODE
-      write (69,*) 'Export_Format Default'      
-      close(69)
-      else
-          write (*,*) "PS Point #",K," done."
-      endif
-      enddo
-
-!!
-!!      Copy down here (or read in) the four momenta as a string. 
-!!      
-!!
-!      buff(1)=" 1   0.5630480E+04  0.0000000E+00  0.0000000E+00  0.5630480E+04"
-!      buff(2)=" 2   0.5630480E+04  0.0000000E+00  0.0000000E+00 -0.5630480E+04"
-!      buff(3)=" 3   0.5466073E+04  0.4443190E+03  0.2446331E+04 -0.4864732E+04"
-!      buff(4)=" 4   0.8785819E+03 -0.2533886E+03  0.2741971E+03  0.7759741E+03"
-!      buff(5)=" 5   0.4916306E+04 -0.1909305E+03 -0.2720528E+04  0.4088757E+04"
-!!
-!!      Here the k,E,px,py,pz are read from the string into the momenta array.
-!!      k=1,2          : incoming
-!!      k=3,nexternal  : outgoing
-!!
-!      do i=1,nexternal
-!         read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i)
-!      enddo
-!
-!!- print the momenta out
-!
-!      do i=1,nexternal
-!         write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i), 
-!     .dsqrt(dabs(DOT(p(0,i),p(0,i))))
-!      enddo
-!
-!      CALL SLOOPMATRIX(P,MATELEM)
-!
-!      write (*,*) "-------------------------------------------------"
-!      write (*,*) "Matrix element = ", MATELEM(1), " GeV^",
-!      &-(2*nexternal-8)      
-!      write (*,*) "-------------------------------------------------"
-
-      end
-      
-        
-        
-        
-         double precision function dot(p1,p2)
-C****************************************************************************
-C     4-Vector Dot product
-C****************************************************************************
-      implicit none
-      double precision p1(0:3),p2(0:3)
-      dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
-      end
-
-
-        SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P)
-C---- auxiliary function to change convention between madgraph and rambo
-c---- four momenta.         
-        IMPLICIT NONE
-        INCLUDE "nexternal.inc"
-C        ARGUMENTS
-        REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT
-C         LOCAL
-         INTEGER I
-         REAL*8 etot2,mom,m1,m2,e1,e2
-
-         ETOT2=energy**2
-         m1=pmass(1)
-         m2=pmass(2)
-         mom=(Etot2**2 - 2*Etot2*m1**2 + m1**4 -
-     -        2*Etot2*m2**2 - 2*m1**2*m2**2 + m2**4)/(4.*Etot2)
-         mom=dsqrt(mom)
-         e1=DSQRT(mom**2+m1**2)
-         e2=DSQRT(mom**2+m2**2)
-c         write (*,*) e1+e2,mom
-
-         if(nincoming.eq.2) then
-
-            P(0,1)=e1
-            P(1,1)=0d0
-            P(2,1)=0d0
-            P(3,1)=mom
-
-            P(0,2)=e2
-            P(1,2)=0d0
-            P(2,2)=0d0
-            P(3,2)=-mom
-             
-            call rambo(nexternal-2,energy,pmass(3),prambo,WGT)
-            DO I=3, NEXTERNAL
-               P(0,I)=PRAMBO(4,I-2)      
-               P(1,I)=PRAMBO(1,I-2)
-               P(2,I)=PRAMBO(2,I-2)
-               P(3,I)=PRAMBO(3,I-2)      
-            ENDDO
-             
-          elseif(nincoming.eq.1) then 
-             
-             P(0,1)=energy
-             P(1,1)=0d0
-             P(2,1)=0d0
-             P(3,1)=0d0
-             
-             call rambo(nexternal-1,energy,pmass(2),prambo,WGT)
-             DO I=2, NEXTERNAL
-                P(0,I)=PRAMBO(4,I-1)      
-                P(1,I)=PRAMBO(1,I-1)
-                P(2,I)=PRAMBO(2,I-1)
-                P(3,I)=PRAMBO(3,I-1)      
-             ENDDO
-          endif
-          
-        RETURN
-        END
-      
-
-      SUBROUTINE RAMBO(N,ET,XM,P,WT)
-***********************************************************************
-*                       RAMBO                                         *
-*    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)                  *
-*                                                                     *
-*    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR                *
-*    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING                 *
-*    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS                      *
-*    -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90)    *
-*                                                                     *
-*    N  = NUMBER OF PARTICLES                                         *
-*    ET = TOTAL CENTRE-OF-MASS ENERGY                                 *
-*    XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming )                 *
-*    P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) )            *
-*    WT = WEIGHT OF THE EVENT                                         *
-***********************************************************************
-      IMPLICIT REAL*8(A-H,O-Z)
-      INCLUDE "nexternal.inc"
-      DIMENSION XM(NEXTERNAL-NINCOMING),P(4,NEXTERNAL-NINCOMING)
-      DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4),
-     .   B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING),
-     .   E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5)
-      SAVE ACC,ITMAX,IBEGIN,IWARN
-      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
-*
-* INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
-      IF(IBEGIN.NE.0) GOTO 103
-      IBEGIN=1
-      TWOPI=8.*DATAN(1.D0)
-      PO2LOG=LOG(TWOPI/4.)
-      Z(2)=PO2LOG
-      DO 101 K=3,(NEXTERNAL-NINCOMING-1)
-  101 Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
-      DO 102 K=3,(NEXTERNAL-NINCOMING-1)
-  102 Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
-*
-* CHECK ON THE NUMBER OF PARTICLES
-  103 IF(N.GT.1.AND.N.LT.101) GOTO 104
-      PRINT 1001,N
-      STOP
-*
-* CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
-  104 XMT=0.
-      NM=0
-      DO 105 I=1,N
-      IF(XM(I).NE.0.D0) NM=NM+1
-  105 XMT=XMT+ABS(XM(I))
-      IF(XMT.LE.ET) GOTO 201
-      PRINT 1002,XMT,ET
-      STOP
-*
-* THE PARAMETER VALUES ARE NOW ACCEPTED
-*
-* GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
-  201 DO 202 I=1,N
-         r1=rn(1)
-      C=2.*r1-1.
-      S=SQRT(1.-C*C)
-      F=TWOPI*RN(2)
-      r1=rn(3)
-      r2=rn(4)
-      Q(4,I)=-LOG(r1*r2)
-      Q(3,I)=Q(4,I)*C
-      Q(2,I)=Q(4,I)*S*COS(F)
-  202 Q(1,I)=Q(4,I)*S*SIN(F)
-*
-* CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
-      DO 203 I=1,4
-  203 R(I)=0.
-      DO 204 I=1,N
-      DO 204 K=1,4
-  204 R(K)=R(K)+Q(K,I)
-      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
-      DO 205 K=1,3
-  205 B(K)=-R(K)/RMAS
-      G=R(4)/RMAS
-      A=1./(1.+G)
-      X=ET/RMAS
-*
-* TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
-      DO 207 I=1,N
-      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
-      DO 206 K=1,3
-  206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
-  207 P(4,I)=X*(G*Q(4,I)+BQ)
-*
-* CALCULATE WEIGHT AND POSSIBLE WARNINGS
-      WT=PO2LOG
-      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
-      IF(WT.GE.-180.D0) GOTO 208
-      IF(IWARN(1).LE.5) PRINT 1004,WT
-      IWARN(1)=IWARN(1)+1
-  208 IF(WT.LE. 174.D0) GOTO 209
-      IF(IWARN(2).LE.5) PRINT 1005,WT
-      IWARN(2)=IWARN(2)+1
-*
-* RETURN FOR WEIGHTED MASSLESS MOMENTA
-  209 IF(NM.NE.0) GOTO 210
-* RETURN LOG OF WEIGHT
-      WT=WT
-      RETURN
-*
-* MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
-  210 XMAX=SQRT(1.-(XMT/ET)**2)
-      DO 301 I=1,N
-      XM2(I)=XM(I)**2
-  301 P2(I)=P(4,I)**2
-      ITER=0
-      X=XMAX
-      ACCU=ET*ACC
-  302 F0=-ET
-      G0=0.
-      X2=X*X
-      DO 303 I=1,N
-      E(I)=SQRT(XM2(I)+X2*P2(I))
-      F0=F0+E(I)
-  303 G0=G0+P2(I)/E(I)
-      IF(ABS(F0).LE.ACCU) GOTO 305
-      ITER=ITER+1
-      IF(ITER.LE.ITMAX) GOTO 304
-      PRINT 1006,ITMAX
-      GOTO 305
-  304 X=X-F0/(X*G0)
-      GOTO 302
-  305 DO 307 I=1,N
-      V(I)=X*P(4,I)
-      DO 306 K=1,3
-  306 P(K,I)=X*P(K,I)
-  307 P(4,I)=E(I)
-*
-* CALCULATE THE MASS-EFFECT WEIGHT FACTOR
-      WT2=1.
-      WT3=0.
-      DO 308 I=1,N
-      WT2=WT2*V(I)/E(I)
-  308 WT3=WT3+V(I)**2/E(I)
-      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
-*
-* RETURN FOR  WEIGHTED MASSIVE MOMENTA
-      WT=WT+WTM
-      IF(WT.GE.-180.D0) GOTO 309
-      IF(IWARN(3).LE.5) PRINT 1004,WT
-      IWARN(3)=IWARN(3)+1
-  309 IF(WT.LE. 174.D0) GOTO 310
-      IF(IWARN(4).LE.5) PRINT 1005,WT
-      IWARN(4)=IWARN(4)+1
-* RETURN LOG OF WEIGHT
-  310 WT=WT
-      RETURN
-*
- 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
- 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
-     . ' SMALLER THAN TOTAL ENERGY =',D15.6)
- 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
- 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
- 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
-     . ' DESIRED ACCURACY =',D15.6)
-      END
-
-      FUNCTION RN(IDUMMY)
-      REAL*8 RN,RAN
-      SAVE INIT
-      DATA INIT /1/
-      IF (INIT.EQ.1) THEN
-        INIT=0
-        CALL RMARIN(1802,9373)
-      END IF
-*
-  10  CALL RANMAR(RAN)
-      IF (RAN.LT.1D-16) GOTO 10
-      RN=RAN
-*
-      END
-
-
-
-      SUBROUTINE RANMAR(RVEC)
-*     -----------------
-* Universal random number generator proposed by Marsaglia and Zaman
-* in report FSU-SCRI-87-50
-* In this version RVEC is a double precision variable.
-      IMPLICIT REAL*8(A-H,O-Z)
-      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
-      COMMON/ RASET2 / IRANMR,JRANMR
-      SAVE /RASET1/,/RASET2/
-      UNI = RANU(IRANMR) - RANU(JRANMR)
-      IF(UNI .LT. 0D0) UNI = UNI + 1D0
-      RANU(IRANMR) = UNI
-      IRANMR = IRANMR - 1
-      JRANMR = JRANMR - 1
-      IF(IRANMR .EQ. 0) IRANMR = 97
-      IF(JRANMR .EQ. 0) JRANMR = 97
-      RANC = RANC - RANCD
-      IF(RANC .LT. 0D0) RANC = RANC + RANCM
-      UNI = UNI - RANC
-      IF(UNI .LT. 0D0) UNI = UNI + 1D0
-      RVEC = UNI
-      END
- 
-      SUBROUTINE RMARIN(IJ,KL)
-*     -----------------
-* Initializing routine for RANMAR, must be called before generating
-* any pseudorandom numbers with RANMAR. The input values should be in
-* the ranges 0<=ij<=31328 ; 0<=kl<=30081
-      IMPLICIT REAL*8(A-H,O-Z)
-      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
-      COMMON/ RASET2 / IRANMR,JRANMR
-      SAVE /RASET1/,/RASET2/
-* This shows correspondence between the simplified input seeds IJ, KL
-* and the original Marsaglia-Zaman seeds I,J,K,L.
-* To get the standard values in the Marsaglia-Zaman paper (i=12,j=34
-* k=56,l=78) put ij=1802, kl=9373
-      I = MOD( IJ/177 , 177 ) + 2
-      J = MOD( IJ     , 177 ) + 2
-      K = MOD( KL/169 , 178 ) + 1
-      L = MOD( KL     , 169 )
-      DO 300 II = 1 , 97
-        S =  0D0
-        T = .5D0
-        DO 200 JJ = 1 , 24
-          M = MOD( MOD(I*J,179)*K , 179 )
-          I = J
-          J = K
-          K = M
-          L = MOD( 53*L+1 , 169 )
-          IF(MOD(L*M,64) .GE. 32) S = S + T
-          T = .5D0*T
-  200   CONTINUE
-        RANU(II) = S
-  300 CONTINUE
-      RANC  =   362436D0 / 16777216D0
-      RANCD =  7654321D0 / 16777216D0
-      RANCM = 16777213D0 / 16777216D0
-      IRANMR = 97
-      JRANMR = 33
-      END
-
-
-
-
-
-

=== removed file 'Template/loop_material/StandAlone/SubProcesses/check_sa_loop_induced.f'
--- Template/loop_material/StandAlone/SubProcesses/check_sa_loop_induced.f	2013-12-13 07:04:16 +0000
+++ Template/loop_material/StandAlone/SubProcesses/check_sa_loop_induced.f	1970-01-01 00:00:00 +0000
@@ -1,574 +0,0 @@
-      PROGRAM DRIVER
-!**************************************************************************
-!     THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT.
-!     IT USES A SIMPLE PHASE SPACE GENERATOR
-!     Fabio Maltoni - 3rd Febraury 2007
-!**************************************************************************
-      IMPLICIT NONE
-!     
-!     CONSTANTS  
-!     
-      REAL*8 ZERO
-      PARAMETER (ZERO=0D0)
-
-      LOGICAL READPS
-      PARAMETER (READPS = .FALSE.)
-
-      INTEGER NPSPOINTS
-      PARAMETER (NPSPOINTS = 4)
-
-!     
-!     INCLUDE FILES
-!     
-!---  the include file with the values of the parameters and masses      
-      INCLUDE "coupl.inc"
-!---  integer nexternal ! number particles (incoming+outgoing) in the me 
-      INCLUDE "nexternal.inc" 
-!---  particle masses
-      REAL*8 PMASS(NEXTERNAL)      
-!---  integer    n_max_cg
-      INCLUDE "ngraphs.inc"     !how many diagrams (could be useful to know...)
-
-!     
-!     LOCAL
-!     
-      INTEGER I,J,K
-      REAL*8 P(0:3,NEXTERNAL)   ! four momenta. Energy is the zeroth component.
-      REAL*8 SQRTS,MATELEM(3)           ! sqrt(s)= center of mass energy 
-      REAL*8 PIN(0:3), POUT(0:3)
-      CHARACTER*120 BUFF(NEXTERNAL)
-      REAL*8 PREC_FOUND
-      INTEGER RETURNCODE, HUNDREDS, TENS, UNITS
-!     
-!     EXTERNAL
-!     
-      REAL*8 DOT
-      EXTERNAL DOT
-      
-!-----
-!     BEGIN CODE
-!-----
-!     
-!---  INITIALIZATION CALLS
-!     
-!---  Call to initialize the values of the couplings, masses and widths 
-!     used in the evaluation of the matrix element. The primary parameters of the
-!     models are read from Cards/param_card.dat. The secondary parameters are calculated
-!     in Source/MODEL/couplings.f. The values are stored in common blocks that are listed
-!     in coupl.inc .
-
-      call setpara('param_card.dat')  !first call to setup the paramaters
-      include "pmass.inc"             !set up masses
-
-!---  Now use a simple multipurpose PS generator (RAMBO) just to get a 
-!     RANDOM set of four momenta of given masses pmass(i) to be used to evaluate 
-!     the madgraph matrix-element.       
-!     Alternatevely, here the user can call or set the four momenta at his will, see below.
-!           
-      IF(nincoming.EQ.1) THEN
-         SQRTS=PMASS(1)
-      ELSE
-         SQRTS=1000d0              !CMS energy in GEV
-      ENDIF
-
-      call printout()
-
-      do K=1,NPSPOINTS
-
-      if(READPS) then
-        open(967, file="PS.input", err=976, status='OLD', action='READ')
-        do i=1,NExternal
-          read(967,*,end=978) P(0,i),P(1,i),P(2,i),P(3,i)
-        enddo
-       goto 978
-  976  continue
-         stop 'Could not read the PS.input phase-space point.'
-  978  continue
-       close(967)
-      else
-        if ((nincoming.eq.2).and.((nexternal - nincoming .eq.1))) then
-          if (pmass(3).eq.0.0d0) then
-              stop 'Cannot generate 2>1 kin. config. with m3=0.0d0'
-          else
-              ! deal with the case of only one particle in the final
-              ! state
-               p(0,1) = pmass(3)/2d0
-               p(1,1) = 0d0
-               p(2,1) = 0d0
-               p(3,1) = pmass(3)/2d0
-               if (pmass(1).GT.0d0) then
-                 p(3,1) = dsqrt(pmass(3)**2/4d0 - pmass(1)**2)
-               endif
-               p(0,2) = pmass(3)/2d0
-               p(1,2) = 0d0
-               p(2,2) = 0d0
-               p(3,2) = -pmass(3)/2d0
-               if (pmass(2) > 0d0) then
-                 p(3,2) = -dsqrt(pmass(3)**2/4d0 - pmass(1)**2)
-               endif
-               p(0,3) = pmass(3)
-               p(1,3) = 0d0
-               p(2,3) = 0d0
-               p(3,3) = 0d0
-          endif
-        else
-          CALL GET_MOMENTA(SQRTS,PMASS,P)
-        endif
-      endif
-
-      do i=0,3
-        PIN(i)=0.0d0
-        do j=1,nincoming
-          PIN(i)=PIN(i)+p(i,j)
-        enddo
-      enddo
-
-!---  In standalone mode, always use sqrt_s as the renormalization
-!     scale.
-      SQRTS=dsqrt(dabs(DOT(PIN(0),PIN(0))))
-
-      MU_R=SQRTS
-
-!---  Update the couplings with the new MU_R
-      CALL UPDATE_AS_PARAM()
-
-!-- Possible tweaks to the run
-!     CALL FORCE_STABILITY_CHECK(.TRUE.)
-
-!     
-!     Now we can call the matrix element!
-!
-      CALL SLOOPMATRIX_THRES(P,MATELEM,-1.0d0,
-     &PREC_FOUND,RETURNCODE)
-!
-!        write the information on the four momenta 
-!
-      if (K.eq.NPSPOINTS) then
-      write (*,*)
-      write (*,*) " Phase space point:"
-      write (*,*) "--------------------------------------------------",
-     &"---------------------------"
-      write (*,*)  "n        E             px             py         ",
-     &"pz               m "
-      do i=1,nexternal
-         write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i), 
-     .dsqrt(dabs(DOT(p(0,i),p(0,i))))
-      enddo
-      write (*,*) "--------------------------------------------------",
-     &"---------------------------"
-      UNITS=MOD(RETURNCODE,10)
-      TENS=(MOD(RETURNCODE,100)-UNITS)/10
-      HUNDREDS=(RETURNCODE-TENS*10-UNITS)/100
-      if (HUNDREDS.eq.1) then
-        if (TENS.eq.3.or.TENS.eq.4) then
-          write(*,*) 'Unknown numerical stability because ',
-     &' MadLoop is in the initialization stage.'
-        else
-          write(*,*) 'Unknown numerical stability, due to ',
-     &'the value of CTModeRun set in MadLoopParams.dat.'
-        endif
-      elseif (HUNDREDS.eq.2) then
-        write(*,*) 'Stable kinematic configuration (SPS).'
-      elseif (HUNDREDS.eq.3) then
-        write(*,*) 'Unstable kinematic configuration (UPS).'
-        write(*,*) 'Quadruple precision rescue successful.'
-      elseif (HUNDREDS.eq.4) then
-        write(*,*) 'Exceptional kinematic configuration (EPS).'
-        write(*,*) 'Both double an quadruple precision computations',
-     %' are unstable.'
-      endif
-      if (TENS.eq.2.or.TENS.eq.4) then
-        write(*,*) 'Quadruple precision computation used.'
-      endif
-      if (HUNDREDS.ne.1) then
-        if (PREC_FOUND.ge.0.0d0) then
-          write(*,'(1x,a23,1x,1e10.2)') 'Relative accuracy     =',
-     &PREC_FOUND
-        elseif (PREC_FOUND.eq.0.0d0) then
-          write(*,'(1x,a23,1x,1e10.2,1x,a30)') 'Relative accuracy     ='
-     &,PREC_FOUND,'(i.e. beyond double precision)'
-        else    
-          write(*,*) 'Estimated accuracy could not be',
-     &' computed for an unknown reason.'
-        endif
-      endif
-      write (*,*) "This is a loop induced process, so only the ",
-     &"unnormalized finite part is output"
-      write (*,*) " here. Be aware that all loops",
-     &" are expected to be UV-finite as no "
-      write (*,*) " renormalization prescription is considered. "
-      write (*,*) "---------------------------------------------------",
-     &"--------------------------"
-      write (*,*) "Matrix element finite = ", MATELEM(1), " GeV^",
-     &-(2*nexternal-8)
-      write (*,*) "---------------------------------------------------",
-     &"--------------------------"
-
-      open(69, file="result.dat", err=976, action='WRITE')
-      do i=1,nexternal      
-          write (69,'(a2,1x,5e25.15)') 'PS',P(0,i),P(1,i),P(2,i),P(3,i)
-      enddo
-      write (69,'(a3,1x,i2)') 'EXP',-(2*nexternal-8)
-      write (69,'(a4,1x,1e25.15)') 'BORN',0.0d0
-      write (69,'(a3,1x,1e25.15)') 'FIN',MATELEM(1)
-      write (69,'(a4,1x,1e25.15)') '1EPS',0.0d0
-      write (69,'(a4,1x,1e25.15)') '2EPS',0.0d0
-      write (69,'(a3,1x,1e10.2)')  'ACC',PREC_FOUND
-      write (69,'(a7,1x,i3)')     'RETCODE',RETURNCODE
-      write (69,*) 'Export_Format LoopInduced'
-      close(69)
-      else
-          write (*,*) "PS Point #",K," done."
-      endif
-      enddo
-
-!!
-!!      Copy down here (or read in) the four momenta as a string. 
-!!      
-!!
-!      buff(1)=" 1   0.5630480E+04  0.0000000E+00  0.0000000E+00  0.5630480E+04"
-!      buff(2)=" 2   0.5630480E+04  0.0000000E+00  0.0000000E+00 -0.5630480E+04"
-!      buff(3)=" 3   0.5466073E+04  0.4443190E+03  0.2446331E+04 -0.4864732E+04"
-!      buff(4)=" 4   0.8785819E+03 -0.2533886E+03  0.2741971E+03  0.7759741E+03"
-!      buff(5)=" 5   0.4916306E+04 -0.1909305E+03 -0.2720528E+04  0.4088757E+04"
-!!
-!!      Here the k,E,px,py,pz are read from the string into the momenta array.
-!!      k=1,2          : incoming
-!!      k=3,nexternal  : outgoing
-!!
-!      do i=1,nexternal
-!         read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i)
-!      enddo
-!
-!!- print the momenta out
-!
-!      do i=1,nexternal
-!         write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i), 
-!     .dsqrt(dabs(DOT(p(0,i),p(0,i))))
-!      enddo
-!
-!      CALL SLOOPMATRIX(P,MATELEM)
-!
-!      write (*,*) "-------------------------------------------------"
-!      write (*,*) "Matrix element = ", MATELEM(1), " GeV^",
-!      &-(2*nexternal-8)      
-!      write (*,*) "-------------------------------------------------"
-
-      end
-      
-        
-        
-        
-         double precision function dot(p1,p2)
-C****************************************************************************
-C     4-Vector Dot product
-C****************************************************************************
-      implicit none
-      double precision p1(0:3),p2(0:3)
-      dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
-      end
-
-
-        SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P)
-C---- auxiliary function to change convention between madgraph and rambo
-c---- four momenta.         
-        IMPLICIT NONE
-        INCLUDE "nexternal.inc"
-C        ARGUMENTS
-        REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT
-C         LOCAL
-         INTEGER I
-         REAL*8 etot2,mom,m1,m2,e1,e2
-
-         ETOT2=energy**2
-         m1=pmass(1)
-         m2=pmass(2)
-         mom=(Etot2**2 - 2*Etot2*m1**2 + m1**4 -
-     -        2*Etot2*m2**2 - 2*m1**2*m2**2 + m2**4)/(4.*Etot2)
-         mom=dsqrt(mom)
-         e1=DSQRT(mom**2+m1**2)
-         e2=DSQRT(mom**2+m2**2)
-c         write (*,*) e1+e2,mom
-
-         if(nincoming.eq.2) then
-
-            P(0,1)=e1
-            P(1,1)=0d0
-            P(2,1)=0d0
-            P(3,1)=mom
-
-            P(0,2)=e2
-            P(1,2)=0d0
-            P(2,2)=0d0
-            P(3,2)=-mom
-             
-            call rambo(nexternal-2,energy,pmass(3),prambo,WGT)
-            DO I=3, NEXTERNAL
-               P(0,I)=PRAMBO(4,I-2)      
-               P(1,I)=PRAMBO(1,I-2)
-               P(2,I)=PRAMBO(2,I-2)
-               P(3,I)=PRAMBO(3,I-2)      
-            ENDDO
-             
-          elseif(nincoming.eq.1) then 
-             
-             P(0,1)=energy
-             P(1,1)=0d0
-             P(2,1)=0d0
-             P(3,1)=0d0
-             
-             call rambo(nexternal-1,energy,pmass(2),prambo,WGT)
-             DO I=2, NEXTERNAL
-                P(0,I)=PRAMBO(4,I-1)      
-                P(1,I)=PRAMBO(1,I-1)
-                P(2,I)=PRAMBO(2,I-1)
-                P(3,I)=PRAMBO(3,I-1)      
-             ENDDO
-          endif
-          
-        RETURN
-        END
-      
-
-      SUBROUTINE RAMBO(N,ET,XM,P,WT)
-***********************************************************************
-*                       RAMBO                                         *
-*    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)                  *
-*                                                                     *
-*    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR                *
-*    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING                 *
-*    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS                      *
-*    -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90)    *
-*                                                                     *
-*    N  = NUMBER OF PARTICLES                                         *
-*    ET = TOTAL CENTRE-OF-MASS ENERGY                                 *
-*    XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming )                 *
-*    P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) )            *
-*    WT = WEIGHT OF THE EVENT                                         *
-***********************************************************************
-      IMPLICIT REAL*8(A-H,O-Z)
-      INCLUDE "nexternal.inc"
-      DIMENSION XM(NEXTERNAL-NINCOMING),P(4,NEXTERNAL-NINCOMING)
-      DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4),
-     .   B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING),
-     .   E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5)
-      SAVE ACC,ITMAX,IBEGIN,IWARN
-      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
-*
-* INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
-      IF(IBEGIN.NE.0) GOTO 103
-      IBEGIN=1
-      TWOPI=8.*DATAN(1.D0)
-      PO2LOG=LOG(TWOPI/4.)
-      Z(2)=PO2LOG
-      DO 101 K=3,(NEXTERNAL-NINCOMING-1)
-  101 Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
-      DO 102 K=3,(NEXTERNAL-NINCOMING-1)
-  102 Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
-*
-* CHECK ON THE NUMBER OF PARTICLES
-  103 IF(N.GT.1.AND.N.LT.101) GOTO 104
-      PRINT 1001,N
-      STOP
-*
-* CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
-  104 XMT=0.
-      NM=0
-      DO 105 I=1,N
-      IF(XM(I).NE.0.D0) NM=NM+1
-  105 XMT=XMT+ABS(XM(I))
-      IF(XMT.LE.ET) GOTO 201
-      PRINT 1002,XMT,ET
-      STOP
-*
-* THE PARAMETER VALUES ARE NOW ACCEPTED
-*
-* GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
-  201 DO 202 I=1,N
-         r1=rn(1)
-      C=2.*r1-1.
-      S=SQRT(1.-C*C)
-      F=TWOPI*RN(2)
-      r1=rn(3)
-      r2=rn(4)
-      Q(4,I)=-LOG(r1*r2)
-      Q(3,I)=Q(4,I)*C
-      Q(2,I)=Q(4,I)*S*COS(F)
-  202 Q(1,I)=Q(4,I)*S*SIN(F)
-*
-* CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
-      DO 203 I=1,4
-  203 R(I)=0.
-      DO 204 I=1,N
-      DO 204 K=1,4
-  204 R(K)=R(K)+Q(K,I)
-      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
-      DO 205 K=1,3
-  205 B(K)=-R(K)/RMAS
-      G=R(4)/RMAS
-      A=1./(1.+G)
-      X=ET/RMAS
-*
-* TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
-      DO 207 I=1,N
-      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
-      DO 206 K=1,3
-  206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
-  207 P(4,I)=X*(G*Q(4,I)+BQ)
-*
-* CALCULATE WEIGHT AND POSSIBLE WARNINGS
-      WT=PO2LOG
-      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
-      IF(WT.GE.-180.D0) GOTO 208
-      IF(IWARN(1).LE.5) PRINT 1004,WT
-      IWARN(1)=IWARN(1)+1
-  208 IF(WT.LE. 174.D0) GOTO 209
-      IF(IWARN(2).LE.5) PRINT 1005,WT
-      IWARN(2)=IWARN(2)+1
-*
-* RETURN FOR WEIGHTED MASSLESS MOMENTA
-  209 IF(NM.NE.0) GOTO 210
-* RETURN LOG OF WEIGHT
-      WT=WT
-      RETURN
-*
-* MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
-  210 XMAX=SQRT(1.-(XMT/ET)**2)
-      DO 301 I=1,N
-      XM2(I)=XM(I)**2
-  301 P2(I)=P(4,I)**2
-      ITER=0
-      X=XMAX
-      ACCU=ET*ACC
-  302 F0=-ET
-      G0=0.
-      X2=X*X
-      DO 303 I=1,N
-      E(I)=SQRT(XM2(I)+X2*P2(I))
-      F0=F0+E(I)
-  303 G0=G0+P2(I)/E(I)
-      IF(ABS(F0).LE.ACCU) GOTO 305
-      ITER=ITER+1
-      IF(ITER.LE.ITMAX) GOTO 304
-      PRINT 1006,ITMAX
-      GOTO 305
-  304 X=X-F0/(X*G0)
-      GOTO 302
-  305 DO 307 I=1,N
-      V(I)=X*P(4,I)
-      DO 306 K=1,3
-  306 P(K,I)=X*P(K,I)
-  307 P(4,I)=E(I)
-*
-* CALCULATE THE MASS-EFFECT WEIGHT FACTOR
-      WT2=1.
-      WT3=0.
-      DO 308 I=1,N
-      WT2=WT2*V(I)/E(I)
-  308 WT3=WT3+V(I)**2/E(I)
-      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
-*
-* RETURN FOR  WEIGHTED MASSIVE MOMENTA
-      WT=WT+WTM
-      IF(WT.GE.-180.D0) GOTO 309
-      IF(IWARN(3).LE.5) PRINT 1004,WT
-      IWARN(3)=IWARN(3)+1
-  309 IF(WT.LE. 174.D0) GOTO 310
-      IF(IWARN(4).LE.5) PRINT 1005,WT
-      IWARN(4)=IWARN(4)+1
-* RETURN LOG OF WEIGHT
-  310 WT=WT
-      RETURN
-*
- 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
- 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
-     . ' SMALLER THAN TOTAL ENERGY =',D15.6)
- 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
- 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
- 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
-     . ' DESIRED ACCURACY =',D15.6)
-      END
-
-      FUNCTION RN(IDUMMY)
-      REAL*8 RN,RAN
-      SAVE INIT
-      DATA INIT /1/
-      IF (INIT.EQ.1) THEN
-        INIT=0
-        CALL RMARIN(1802,9373)
-      END IF
-*
-  10  CALL RANMAR(RAN)
-      IF (RAN.LT.1D-16) GOTO 10
-      RN=RAN
-*
-      END
-
-
-
-      SUBROUTINE RANMAR(RVEC)
-*     -----------------
-* Universal random number generator proposed by Marsaglia and Zaman
-* in report FSU-SCRI-87-50
-* In this version RVEC is a double precision variable.
-      IMPLICIT REAL*8(A-H,O-Z)
-      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
-      COMMON/ RASET2 / IRANMR,JRANMR
-      SAVE /RASET1/,/RASET2/
-      UNI = RANU(IRANMR) - RANU(JRANMR)
-      IF(UNI .LT. 0D0) UNI = UNI + 1D0
-      RANU(IRANMR) = UNI
-      IRANMR = IRANMR - 1
-      JRANMR = JRANMR - 1
-      IF(IRANMR .EQ. 0) IRANMR = 97
-      IF(JRANMR .EQ. 0) JRANMR = 97
-      RANC = RANC - RANCD
-      IF(RANC .LT. 0D0) RANC = RANC + RANCM
-      UNI = UNI - RANC
-      IF(UNI .LT. 0D0) UNI = UNI + 1D0
-      RVEC = UNI
-      END
- 
-      SUBROUTINE RMARIN(IJ,KL)
-*     -----------------
-* Initializing routine for RANMAR, must be called before generating
-* any pseudorandom numbers with RANMAR. The input values should be in
-* the ranges 0<=ij<=31328 ; 0<=kl<=30081
-      IMPLICIT REAL*8(A-H,O-Z)
-      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
-      COMMON/ RASET2 / IRANMR,JRANMR
-      SAVE /RASET1/,/RASET2/
-* This shows correspondence between the simplified input seeds IJ, KL
-* and the original Marsaglia-Zaman seeds I,J,K,L.
-* To get the standard values in the Marsaglia-Zaman paper (i=12,j=34
-* k=56,l=78) put ij=1802, kl=9373
-      I = MOD( IJ/177 , 177 ) + 2
-      J = MOD( IJ     , 177 ) + 2
-      K = MOD( KL/169 , 178 ) + 1
-      L = MOD( KL     , 169 )
-      DO 300 II = 1 , 97
-        S =  0D0
-        T = .5D0
-        DO 200 JJ = 1 , 24
-          M = MOD( MOD(I*J,179)*K , 179 )
-          I = J
-          J = K
-          K = M
-          L = MOD( 53*L+1 , 169 )
-          IF(MOD(L*M,64) .GE. 32) S = S + T
-          T = .5D0*T
-  200   CONTINUE
-        RANU(II) = S
-  300 CONTINUE
-      RANC  =   362436D0 / 16777216D0
-      RANCD =  7654321D0 / 16777216D0
-      RANCM = 16777213D0 / 16777216D0
-      IRANMR = 97
-      JRANMR = 33
-      END
-
-
-
-
-
-

=== removed file 'Template/loop_material/StandAlone/SubProcesses/makefile'
--- Template/loop_material/StandAlone/SubProcesses/makefile	2012-11-09 17:13:53 +0000
+++ Template/loop_material/StandAlone/SubProcesses/makefile	1970-01-01 00:00:00 +0000
@@ -1,29 +0,0 @@
-include ../../Source/make_opts
-
-LIBDIR = ../../lib/
-PROG   = check
-STABCHECKDRIVER = StabilityCheckDriver
-LINKLIBS =  -L$(LIBDIR) -lcts -ldhelas -lmodel
-LIBS =  $(LIBDIR)libcts.$(libext) $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libmodel.$(libext)
-PROCESS= loop_matrix.o improve_ps.o\
-		 $(patsubst %.f,%.o,$(wildcard born_matrix.f)) \
-		 loop_num.o CT_interface.o \
-		 $(patsubst %.f,%.o,$(wildcard polynomial.f)) \
-		 $(patsubst %.f,%.o,$(wildcard MadLoopParamReader.f)) \
-		 $(patsubst %.f,%.o,$(wildcard helas_calls*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard jamp?_calls_*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard mp_born_amps_and_wfs.f)) \
-		 $(patsubst %.f,%.o,$(wildcard mp_compute_loop_coefs.f)) \
-		 $(patsubst %.f,%.o,$(wildcard mp_helas_calls*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard coef_construction_*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard loop_CT_calls_*.f)) \
-		 $(patsubst %.f,%.o,$(wildcard mp_coef_construction_*.f))
-
-$(PROG): check_sa.o $(PROCESS) makefile $(LIBS)
-	$(FC) $(FFLAGS) -o $(PROG) check_sa.o $(PROCESS) $(LINKLIBS)
-
-$(STABCHECKDRIVER):  StabilityCheckDriver.o $(PROCESS) makefile $(LIBS)
-	$(FC) $(FFLAGS) -o $(STABCHECKDRIVER) StabilityCheckDriver.o $(PROCESS) $(LINKLIBS)
-
-clean:
-		@rm -f *.o

=== added file 'Template/loop_material/StandAlone/SubProcesses/makefile.inc'
--- Template/loop_material/StandAlone/SubProcesses/makefile.inc	1970-01-01 00:00:00 +0000
+++ Template/loop_material/StandAlone/SubProcesses/makefile.inc	2014-08-26 07:56:12 +0000
@@ -0,0 +1,74 @@
+ifeq ($(wildcard ../../Source/make_opts),)
+    ROOT = ..
+else
+    ROOT = ../..
+endif
+
+
+if [ -a ../../Source/make_opts]; then ROOT = ../..; fi;
+
+if [ -a ../Source/make_opts]; then ROOT = ../; fi;
+
+include $(ROOT)/Source/make_opts
+LIBDIR = $(ROOT)/lib/
+
+PROG   = check
+OLP    = OLP
+STABCHECKDRIVER = StabilityCheckDriver
+CHECK_SA_BORN_SPLITORDERS = check_sa_born_splitOrders
+LINKLIBS =  -L$(LIBDIR) -lcts -ldhelas -lmodel %(link_tir_libs)s
+LIBS =  $(LIBDIR)libcts.$(libext) $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libmodel.$(libext) %(tir_libs)s
+PROCESS= MadLoopParamReader.o MadLoopCommons.o \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard loop_matrix.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard improve_ps.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard born_matrix.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard CT_interface.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard loop_num.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard polynomial.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard helas_calls*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard jamp?_calls_*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_born_amps_and_wfs.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_compute_loop_coefs.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_helas_calls*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard coef_construction_*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard loop_CT_calls_*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard mp_coef_construction_*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard TIR_interface.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard GOLEM_interface.f))
+
+OLP_PROCESS= MadLoopParamReader.o MadLoopCommons.o \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/loop_matrix.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/improve_ps.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/born_matrix.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/CT_interface.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/loop_num.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/polynomial.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/helas_calls*.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/jamp?_calls_*.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/mp_born_amps_and_wfs.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/mp_compute_loop_coefs.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/mp_helas_calls*.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/coef_construction_*.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/loop_CT_calls_*.f)) \
+                 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/mp_coef_construction_*.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/TIR_interface.f)) \
+		 $(patsubst %(dotf)s,%(doto)s,$(wildcard P*/GOLEM_interface.f))
+
+$(PROG): check_sa.o $(PROCESS) makefile $(LIBS)
+	$(FC) $(FFLAGS) -o $(PROG) check_sa.o $(PROCESS) $(LINKLIBS)
+
+$(STABCHECKDRIVER):  StabilityCheckDriver.o $(PROCESS) makefile $(LIBS)
+	$(FC) $(FFLAGS) -o $(STABCHECKDRIVER) StabilityCheckDriver.o $(PROCESS) $(LINKLIBS)
+
+# The program below is not essential but just an helpful one to run the born only
+$(CHECK_SA_BORN_SPLITORDERS):  check_sa_born_splitOrders.o $(patsubst %(dotf)s,%(doto)s,$(wildcard *born_matrix.f)) makefile $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libmodel.$(libext)
+	$(FC) $(FFLAGS) -o $(CHECK_SA_BORN_SPLITORDERS) check_sa_born_splitOrders.o $(patsubst %(dotf)s,%(doto)s,$(wildcard *born_matrix.f)) -L$(LIBDIR) -ldhelas -lmodel
+
+%(doto)s : %(dotf)s
+	$(FC) $(FFLAGS) -c $< -o $@ %(tir_include)s
+
+$(OLP): $(OLP_PROCESS) $(LIBS)
+	$(FC) -shared $(OLP_PROCESS) -o libMadLoop.$(dylibext) $(LINKLIBS)
+
+clean:
+		@rm -f *.o

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2014-08-06 16:48:15 +0000
+++ UpdateNotes.txt	2014-08-26 07:56:12 +0000
@@ -1,5 +1,6 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
+<<<<<<< TREE
 2.1.3()XX/XX/XX)OM: Fix various compatibility problem created by the LHEFv3 version (Thanks to S. Brochet)
                 OM: Fix a bug for MadSpin in gridpack mode
                 OM: Add a routine to check the validity of LHE file (check_event command)
@@ -8,6 +9,26 @@
                 OM: Change colorflow information of LHE file in presence of two epsilon_ijk
                     since PY8 was not able to handle such flow in that format.
 
+=======
+2.2.0(21/07/14) VH: General mixed order corrections in MadLoop (only), including QCD/EW/QED and 
+                    the UFO@NLO model 'loop_qcd_qed_sm'.
+			    VH: Re-design of both the tree and MadLoop matrix elements output to compute 
+				    contributions of different coupling order combinations independently from one another.
+			 VH+HS: Tensor integral reduction as implemented in PJFry and IREGI readily available and 
+			        interfaced to MadLoop's output.
+			    VH: Re-structuring of MadLoop's standalone output so as to easily create a single dynamic 
+				    library including many processes at once. Useful for interfacing MadLoop to other MC's 
+					and already working with Sherpa.
+			 VH+HS: This branch contains all the fixes for proper treatment of the latest BSM@NLO models 
+			        produced by FeynRules@NLO. In particular, the fixed related to the presence of majorana 
+					particles in loop ME's.
+
+          PT: Redesigned shower_card.dat and eliminated modbos options for Herwig6          
+          RF: Change the SPINUP information in the NLO LHEF from 0 to 9 (i.e. sum over helicities)
+          RF: Fixed a bug in the check on the determination of the conflicting BWs.
+          MZ: enabled LHAPDF6 interface 
+
+>>>>>>> MERGE-SOURCE
 2.1.2(03/07/14) OM: Fix a bug in ALOHA in presence of customized propagator (Thanks Saurabh)
                 OM: Fixing some compilation issue with MadWeight (Thanks A. Pin)
                 OM: Fixing a bug preventing MadWidth to run due to the model prefixing (depending

=== removed symlink 'VERSION'
=== target was u'madgraph/VERSION'
=== modified file 'aloha/aloha_lib.py'
--- aloha/aloha_lib.py	2014-01-08 20:14:42 +0000
+++ aloha/aloha_lib.py	2014-08-26 07:56:12 +0000
@@ -938,6 +938,7 @@
                     if not veto or not scalar.contains(veto):
                         scalar = scalar.simplify()
                         prefactor = 1
+
                         if hasattr(scalar, 'vartype') and scalar.prefactor not in  [1,-1]:
                             prefactor = scalar.prefactor
                             scalar.prefactor = 1
@@ -1048,7 +1049,16 @@
             else:
                 raise self.LorentzObjectRepresentationError("There is no key of (0,) in representation.")                    
         else:
-            self[(0,)] = representation
+            if isinstance(representation,dict):
+                try:
+                    self[(0,)] = representation[(0,)]
+                except Exception:
+                    if representation:
+		   	            raise LorentzObjectRepresentation.LorentzObjectRepresentationError("There is no key of (0,) in representation.")
+                    else:
+			            self[(0,)] = 0
+            else:
+                self[(0,)] = representation
 
     def __str__(self):
         """ string representation """

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2014-08-04 14:56:54 +0000
+++ aloha/aloha_writers.py	2014-08-26 07:56:12 +0000
@@ -33,7 +33,6 @@
 
             
     def __init__(self, abstract_routine, dirpath):
-        
         if aloha.loop_mode:
             self.momentum_size = 4
         else:
@@ -53,21 +52,19 @@
         self.routine = abstract_routine
         self.tag = self.routine.tag
         self.name = name
+
         self.particles =  [self.type_to_variable[spin] for spin in \
                           abstract_routine.spins]
-        
+
         self.offshell = abstract_routine.outgoing # position of the outgoing in particle list        
         self.outgoing = self.offshell             # expected position for the argument list
         if 'C%s' %((self.outgoing + 1) // 2) in self.tag:
             #flip the outgoing tag if in conjugate
             self.outgoing = self.outgoing + self.outgoing % 2 - (self.outgoing +1) % 2
-        
         self.outname = '%s%s' % (self.particles[self.outgoing -1], \
                                                                self.outgoing)
-        
         #initialize global helper routine
         self.declaration = Declaration_list()
-
                                    
                                        
     def pass_to_HELAS(self, indices, start=0):
@@ -200,7 +197,8 @@
         # couplings
         if  couplings is None:
             detected_couplings = [name for type, name in self.declaration if name.startswith('COUP')]
-            detected_couplings.sort()
+            coup_sort = lambda x,y: int(x[4:])-int(y[4:])
+            detected_couplings.sort(coup_sort)
             if detected_couplings:
                 couplings = detected_couplings
             else:
@@ -449,7 +447,7 @@
                    'im': 'imag(%s)',
                    'cmath.sqrt':'sqrt(dble(%s))', 
                    'sqrt': 'sqrt(dble(%s))',
-                   'complexconjugate': 'conjg(%s)',
+                   'complexconjugate': 'conjg(dcmplx(%s))',
                    '/' : '{0}/(%s)'.format(one),
                    'pow': '(%s)**(%s)',
                    'log': 'log(dble(%s))',
@@ -689,7 +687,7 @@
                     if number.imag == 1:
                         out = 'CI'
                     elif number.imag == -1:
-                        out = '-CI'
+			            out = '-CI'
                     else: 
                         out = '%s * CI' % self.change_number_format(number.imag)
             else:
@@ -728,8 +726,9 @@
         keys.sort(sort_fct)
         for name in keys:
             fct, objs = self.routine.fct[name]
-            format = ' %s = %s\n' % (name, self.get_fct_format(fct))
-            try:
+
+	    format = ' %s = %s\n' % (name, self.get_fct_format(fct))
+	    try:
                 text = format % ','.join([self.write_obj(obj) for obj in objs])
             except TypeError:
                 text = format % tuple([self.write_obj(obj) for obj in objs])
@@ -909,10 +908,10 @@
         # position of the outgoing in particle list        
         self.l_id = [int(c[1:]) for c in abstract_routine.tag if c[0] == 'L'][0] 
         self.l_helas_id = self.l_id   # expected position for the argument list
-        if 'C%s' %((self.outgoing + 1) // 2) in abstract_routine.tag:
+        if 'C%s' %((self.l_id + 1) // 2) in abstract_routine.tag:
             #flip the outgoing tag if in conjugate
-            self.l_helas_id += self.l_id % 2 - (self.l_id +1) % 2 
-        
+            self.l_helas_id += self.l_id % 2 - (self.l_id +1) % 2
+         
 
     def define_expression(self):
         """Define the functions in a 100% way """
@@ -1015,30 +1014,25 @@
         """define a list with the string of object required as incoming argument"""
 
         conjugate = [2*(int(c[1:])-1) for c in self.tag if c[0] == 'C']
-        call_arg = [('list_complex', 'P%s'% self.l_helas_id)] #incoming argument of the routine
+        call_arg = []
+        #incoming argument of the routine
+        call_arg.append( ('list_complex', 'P%s'% self.l_helas_id) )
+        
         self.declaration.add(call_arg[0])
         
         for index,spin in enumerate(self.particles):
-            if self.offshell == index + 1:
+            if self.outgoing == index + 1:
                 continue
             if self.l_helas_id == index + 1:
                 continue
-            
-            if index in conjugate:
-                index2, spin2 = index+1, self.particles[index+1]
-                call_arg.append(('complex','%s%d' % (spin2, index2 +1)))
-                #call_arg.append('%s%d' % (spin, index +1)) 
-            elif index-1 in conjugate:
-                index2, spin2 = index-1, self.particles[index-1]
-                call_arg.append(('complex','%s%d' % (spin2, index2 +1)))
-            else:
-                call_arg.append(('complex','%s%d' % (spin, index +1)))
+            call_arg.append(('complex','%s%d' % (spin, index +1)))
             self.declaration.add(('list_complex', call_arg[-1][-1])) 
         
         # couplings
         if couplings is None:
             detected_couplings = [name for type, name in self.declaration if name.startswith('COUP')]
-            detected_couplings.sort()
+            coup_sort = lambda x,y: int(x[4:])-int(y[4:])
+            detected_couplings.sort(coup_sort)
             if detected_couplings:
                 couplings = detected_couplings
             else:
@@ -1059,7 +1053,7 @@
                 self.declaration.add(('double','W%s' % self.outgoing))
             
         self.call_arg = call_arg
-                
+        
         return call_arg
 
     def get_momenta_txt(self):
@@ -1315,9 +1309,9 @@
                    'im': 'imag(%s)',
                    'cmath.sqrt':'sqrt(%s)', 
                    'sqrt': 'sqrt(%s)',
-                   'complexconjugate': 'conj(%s)',
+                   'complexconjugate': 'conj(dcmplx(%s))',
                    '/' : '{0}/%s'.format(one),
-                   'abs': 'abs(%s)' 
+                   'abs': 'abs(%s)'
                    }
             
         if fct in self.fct_format:
@@ -2098,7 +2092,6 @@
 class WriterFactory(object):
     
     def __new__(cls, data, language, outputdir, tags):
-        
         language = language.lower()
         if isinstance(data.expr, aloha_lib.SplitCoefficient):
             assert language == 'fortran'
@@ -2106,7 +2099,6 @@
                 return ALOHAWriterForFortranLoopQP(data, outputdir)
             else:
                 return ALOHAWriterForFortranLoop(data, outputdir)
-        
         if language == 'fortran':
             if 'MP' in tags:
                 return ALOHAWriterForFortranQP(data, outputdir)

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2014-08-04 14:56:54 +0000
+++ aloha/create_aloha.py	2014-08-26 07:56:12 +0000
@@ -85,7 +85,6 @@
         
     def write(self, output_dir, language='Fortran', mode='self', combine=True,**opt):
         """ write the content of the object """
-
         writer = aloha_writers.WriterFactory(self, language, output_dir, self.tag)
         text = writer.write(mode=mode, **opt)
         if combine:
@@ -94,8 +93,7 @@
                     text = tuple([old.__add__(new)  for old, new in zip(text, 
                              writer.write_combined(grouped, mode=mode+'no_include', **opt))])
                 else:
-                    text += writer.write_combined(grouped, mode=mode+'no_include', **opt)
-                
+                    text += writer.write_combined(grouped, mode=mode+'no_include', **opt)        
         if aloha.mp_precision and 'MP' not in self.tag:
             self.tag.append('MP')
             text += self.write(output_dir, language, mode, **opt)
@@ -288,7 +286,7 @@
             aloha_lib.KERNEL.use_tag = set(self.kernel_tag) 
         else:
             lorentz = copy.copy(self.routine_kernel)
-            aloha_lib.KERNEL.use_tag = set(self.kernel_tag) 
+            aloha_lib.KERNEL.use_tag = set(self.kernel_tag)
         for (i, spin ) in enumerate(self.spins):   
             id = i + 1
             #Check if this is the outgoing particle
@@ -478,8 +476,11 @@
     
     def compute_loop_coefficient(self, lorentz, outgoing):
         
-                    
-        l_in = [int(tag[1:]) for tag in self.tag][0]
+
+        l_in = [int(tag[1:]) for tag in self.tag if tag.startswith('L')][0]
+        if (l_in + 1) // 2 in self.conjg:
+            #flip the outgoing tag if in conjugate                                                                                                                                         
+            l_in = l_in + l_in % 2 - (l_in +1) % 2                    
         assert l_in != outgoing, 'incoming Open Loop can not be the outcoming one'
         
         # modify the expression for the momenta
@@ -820,12 +821,10 @@
                                     self[(realname, outgoing)].add_combine(m)
                                 except Exception,error:
                                     self[(realname, self.symmetries[lorentz.name][outgoing])].add_combine(m)          
-                    
-                    
+                       
         if save:
             self.save()
     
-    
     def add_Lorentz_object(self, lorentzlist):
         """add a series of Lorentz structure created dynamically"""
         
@@ -844,12 +843,15 @@
         self.look_for_symmetries()
         # reorganize the data (in order to use optimization for a given lorentz
         #structure
+        aloha.loop_mode = False
+	    # self.explicit_combine = False
         request = {}
 
         for list_l_name, tag, outgoing in data:
             #allow tag to have integer for retro-compatibility
             all_tag = tag[:]
             conjugate = [i for i in tag if isinstance(i, int)]
+            
             tag =  [i for i in tag if isinstance(i, str) and not i.startswith('P')]
             tag = tag + ['C%s'%i for i in conjugate]             
             tag = tag + [i for i in all_tag if isinstance(i, str) and  i.startswith('P')] 
@@ -859,7 +861,7 @@
             if loop:
                 aloha.loop_mode = True
                 self.explicit_combine = True
-
+                       
             for l_name in list_l_name:
                 try:
                     request[l_name][conjugate].append((outgoing,tag))
@@ -868,7 +870,7 @@
                         request[l_name][conjugate] = [(outgoing,tag)]
                     except Exception:
                         request[l_name] = {conjugate: [(outgoing,tag)]}
-                        
+                           
         # Loop on the structure to build exactly what is request
         for l_name in request:
             lorentz = eval('self.model.lorentz.%s' % l_name)
@@ -898,7 +900,8 @@
                     conjg_builder = builder.define_conjugate_builder(conjg)
                     # Compute routines
                     self.compute_aloha(conjg_builder, symmetry=lorentz.name,
-                                        routines=routines)
+                                         routines=routines)
+            
         
         # Build mutiple lorentz call
         for list_l_name, tag, outgoing in data:
@@ -940,8 +943,10 @@
                         conjg_builder = builder.define_conjugate_builder(conjg)
                         # Compute routines
                         self.compute_aloha(conjg_builder, symmetry=lorentz.name,
-                                        routines=routines)             
-                        
+                                        routines=routines)
+                      
+  
+                            
     def compute_aloha(self, builder, symmetry=None, routines=None, tag=[]):
         """ define all the AbstractRoutine linked to a given lorentz structure
         symmetry authorizes to use the symmetry of anoter lorentz structure.
@@ -971,7 +976,7 @@
                 wavefunction = builder.compute_routine(outgoing, tag)
                 #Store the information
                 self.set(realname, outgoing, wavefunction)
-            
+          
 
     def compute_aloha_without_kernel(self, builder, symmetry=None, routines=None):
         """define all the AbstractRoutine linked to a given lorentz structure
@@ -992,10 +997,9 @@
 
     def write(self, output_dir, language):
         """ write the full set of Helicity Routine in output_dir"""
-
         for abstract_routine in self.values():
             abstract_routine.write(output_dir, language)
-        
+
         for routine in self.external_routines:
             self.locate_external(routine, language, output_dir)
 

=== modified file 'input/.mg5_configuration_default.txt'
--- input/.mg5_configuration_default.txt	2014-06-30 11:59:30 +0000
+++ input/.mg5_configuration_default.txt	2014-08-26 07:56:12 +0000
@@ -90,7 +90,7 @@
 # run_mode = 2
 
 #! Cluster Type [pbs|sge|condor|lsf|ge|slurm] Use for cluster run only
-#!  And cluster queue
+#!  And cluster queue (or partition for slurm)
 # cluster_type = condor
 # cluster_queue = madgraph
 
@@ -151,6 +151,20 @@
 #! Set what OLP to use for the loop ME generation
 # OLP = MadLoop
 
+#! Set the PJFry++ lib libpjfry.a directory
+#! if auto: try to find it automatically on the system (default)
+#! if '' or None: not use PJFry++
+#! if pjfry=/PATH/TO/pjfry/lib: use that specific installation path for PJFry++
+# pjfry = auto
+
+#! Set the Golem95 lib libgolem.a directory
+#! It only supports version higher than 1.3.0
+#! if auto: try to find it automatically on the system (default)
+#! if '' or None: not use Golem95
+#! if golem=/PATH/TO/golem/lib: use that speficif installation path for Golem95
+# golem = auto
+
+
 #! Set how MadLoop dependencies (such as CutTools) should be handled
 #!  > external : ML5 places a link to the MG5_aMC-wide libraries
 #!  > internal : ML5 copies all dependencies in the output so that it is independent
@@ -163,3 +177,4 @@
 # syscalc_path = ./SysCalc
 
 
+

=== removed file 'input/proc_card_default.dat'
--- input/proc_card_default.dat	2013-12-09 01:31:26 +0000
+++ input/proc_card_default.dat	1970-01-01 00:00:00 +0000
@@ -1,40 +0,0 @@
-#************************************************************
-#*                        MadGraph 5                        *
-#*                                                          *
-#*                *                       *                 *
-#*                  *        * *        *                   *
-#*                    * * * * 5 * * * *                     *
-#*                  *        * *        *                   *
-#*                *                       *                 *
-#*                                                          *
-#*                                                          *
-#*    The MadGraph Development Team - Please visit us at    *
-#*    https://server06.fynu.ucl.ac.be/projects/madgraph     *
-#*                                                          *
-#************************************************************
-#*                                                          *
-#*               Command File for MadGraph 5                *
-#*                                                          *
-#*     run as ./bin/mg5  filename                           *
-#*                                                          *
-#************************************************************
-import model sm
-# Define multiparticle labels
-define p = g u c d s u~ c~ d~ s~
-define j = g u c d s u~ c~ d~ s~
-define l+ = e+ mu+
-define l- = e- mu-
-define vl = ve vm vt
-define vl~ = ve~ vm~ vt~
-# Specify process(es) to run
-generate p p > e- ve~ @1
-add process p p > e- ve~ j @2
-add process p p > t t~ @3
-# Output processes to MadEvent directory
-output 
-# This will create a directory PROC_$MODELNAME_$X
-# If you want to specify the path/name of the directory use
-# output PATH
-
-# To generate events, you can go to the created directory and 
-# run ./bin/generate_events
\ No newline at end of file

=== modified file 'madgraph/VERSION'
--- madgraph/VERSION	2014-07-03 13:37:06 +0000
+++ madgraph/VERSION	2014-08-26 07:56:12 +0000
@@ -1,5 +1,5 @@
-version = 2.1.2
-date = 2014-07-03
+version = 2.2.0
+date = 2014-07-21
 
 
 

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2014-06-29 15:15:00 +0000
+++ madgraph/core/base_objects.py	2014-08-26 07:56:12 +0000
@@ -94,7 +94,6 @@
         """Set the value of the property name. First check if value
         is a valid value for the considered property. Return True if the
         value has been correctly set, False otherwise."""
-
         if not __debug__ or force:
             self[name] = value
             return True
@@ -345,6 +344,16 @@
            in argument given the model specified. It is very fast for usual models"""
            
         for int in model['interactions'].get_type('base'):
+            # We discard the interactions with more than one type of orders
+            # contributing because it then doesn't necessarly mean that this
+            # particle (self) is charged under the group corresponding to the
+            # coupling order 'order'. The typical example is in SUSY which 
+            # features a ' photon-gluon-squark-antisquark ' interaction which
+            # has coupling orders QED=1, QCD=1 and would induce the photon
+            # to be considered as a valid particle to circulate in a loop of
+            # type "QCD".
+            if len(int.get('orders'))>1:
+                continue
             if order in int.get('orders').keys() and self.get('pdg_code') in \
               [part.get('pdg_code') for part in int.get('particles')]:
                 return True
@@ -385,7 +394,24 @@
             return - self['color']
         else:
             return self['color']
-
+        
+    def get_charge(self):
+        """Return the charge code with a correct minus sign"""
+
+        if not self['is_part']:
+            return - self['charge']
+        else:
+            return self['charge']
+
+    def get_anti_charge(self):
+        """Return the charge code of the antiparticle with a correct minus sign
+        """
+
+        if self['is_part']:
+            return - self['charge']
+        else:
+            return self['charge']
+        
     def get_name(self):
         """Return the name if particle, antiname if antiparticle"""
 
@@ -453,10 +479,28 @@
         corresponding particle (first one in the list), with the 
         is_part flag set accordingly. None otherwise."""
         
+        assert isinstance(name, str)
+        
         part = self.find_name(name)
         if not part:
+            # Then try to look for a particle with that PDG
+            try:
+                pdg = int(name)
+            except ValueError:
+                return None
+
+            for p in self:
+                if p.get_pdg_code()==pdg:
+                    part = copy.copy(p)
+                    part.set('is_part', True)
+                    return part
+                elif p.get_anti_pdg_code()==pdg:
+                    part = copy.copy(p)
+                    part.set('is_part', False)
+                    return part
+
             return None
-        part = copy.copy(part)     
+        part = copy.copy(part)
           
         if part.get('name') == name:
             part.set('is_part', True)
@@ -1930,6 +1974,13 @@
                 res.append(leg.get('id'))
 
         return res
+    
+    def sort(self,*args, **opts):
+        """Match with FKSLegList"""
+        Opts=copy.copy(opts)
+        if 'pert' in Opts.keys():
+            del Opts['pert']
+        return super(LegList,self).sort(*args, **Opts)
 
 
 #===============================================================================
@@ -2159,6 +2210,25 @@
         coupling_orders['WEIGHTED'] = weight
         self.set('orders', coupling_orders)
 
+    def pass_squared_order_constraints(self, diag_multiplier, squared_orders,
+                                                               sq_orders_types):
+        """ Returns wether the contributiong consisting in the current diagram 
+        multiplied by diag_multiplier passes the *positive* squared_orders 
+        specified ( a dictionary ) of types sq_order_types (a dictionary whose 
+        values are the relational operator used to define the constraint of the
+        order in key)."""
+        
+        for order, value in squared_orders.items():
+            if value<0:
+                continue
+            combined_order = self.get_order(order) + \
+                                                diag_multiplier.get_order(order)
+            if ( sq_orders_types[order]=='==' and combined_order != value ) or \
+               ( sq_orders_types[order] in ['=', '<='] and combined_order > value) or \
+               ( sq_orders_types[order]=='>' and combined_order <= value) :
+                return False
+        return True
+
     def get_order(self, order):
         """Return the order of this diagram. It returns 0 if it is not present."""
 
@@ -2299,6 +2369,40 @@
 
         return max_order
 
+    def apply_negative_sq_order(self, ref_diag_list, order, value, order_type):
+        """ This function returns a fitlered version of the diagram list self
+        which satisfy the negative squared_order constraint 'order' with negative
+        value 'value' and of type 'order_type', assuming that the diagram_list
+        it must be squared against is 'reg_diag_list'. It also returns the
+        new postive target squared order which correspond to this negative order
+        constraint. Example: u u~ > d d~ QED^2<=-2 means that one wants to
+        pick terms only up to the the next-to-leading order contributiong in QED,
+        which is QED=2 in this case, so that target_order=4 is returned."""
+        
+        # First we must compute all contributions to that order
+        target_order = min(ref_diag_list.get_order_values(order))+\
+                                  min(self.get_order_values(order))+2*(-value-1)
+        
+        new_list = self.apply_positive_sq_orders(ref_diag_list, 
+                                       {order:target_order}, {order:order_type})
+        
+        return new_list, target_order
+        
+    def apply_positive_sq_orders(self, ref_diag_list, sq_orders, sq_order_types):
+        """ This function returns a filtered version of self which contain
+        only the diagram which satisfy the positive squared order constraints
+        sq_orders of type sq_order_types and assuming that the diagrams are
+        multiplied with those of the reference diagram list ref_diag_list."""
+                
+        new_diag_list = DiagramList()
+        for tested_diag in self:
+            for ref_diag in ref_diag_list:
+                if tested_diag.pass_squared_order_constraints(ref_diag,
+                                                      sq_orders,sq_order_types):
+                    new_diag_list.append(tested_diag)
+                    break
+        return new_diag_list
+
     def get_min_order(self,order):
         """ Return the order of the diagram in the list with the mimimum coupling
         order for the coupling specified """
@@ -2312,6 +2416,19 @@
 
         return min_order
 
+    def get_order_values(self, order):
+        """ Return the list of possible values appearing in the diagrams of this
+        list for the order given in argument """
+
+        values=set([])
+        for diag in self:
+            if order in diag['orders'].keys():
+                values.add(diag['orders'][order])
+            else:
+                values.add(0)  
+
+        return list(values)
+
 #===============================================================================
 # Process
 #===============================================================================
@@ -2353,10 +2470,25 @@
         # gives the upper bound for the total weighted order of the
         # squared amplitude.
         self['squared_orders'] = {}
+        # The squared order (sqorders) constraints above can either be upper 
+        # bound (<=) or exact match (==) depending on how they were specified
+        # in the user input. This choice is stored in the dictionary below.
+        # Notice that the upper bound is the default
+        self['sqorders_types'] = {}
         self['has_born'] = True
         # The NLO_mode is always None for a tree-level process and can be
         # 'all', 'real', 'virt' for a loop process.
         self['NLO_mode'] = 'tree'
+        # The user might want to have the individual matrix element evaluations
+        # for specific values of the coupling orders. The list below specifies
+        # what are the coupling names which need be individually treated.
+        # For example, for the process p p > j j [] QED=2 (QED=2 is 
+        # then a squared order constraint), then QED will appear in the 
+        # 'split_orders' list so that the subroutine in matrix.f return the
+        # evaluation of the matrix element individually for the pure QCD 
+        # contribution 'QCD=4 QED=0', the pure interference 'QCD=2 QED=2' and
+        # the pure QED contribution of order 'QCD=0 QED=4'.
+        self['split_orders'] = []
 
     def filter(self, name, value):
         """Filter for valid process property values."""
@@ -2369,6 +2501,24 @@
         if name in ['orders', 'overall_orders','squared_orders']:
             Interaction.filter(Interaction(), 'orders', value)
 
+        if name == 'sqorders_types':
+            if not isinstance(value, dict):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid dictionary" % str(value)
+            for order in value.keys()+value.values():
+                if not isinstance(order, str):
+                    raise self.PhysicsObjectError, \
+                          "%s is not a valid string" % str(value)
+
+        if name == 'split_orders':
+            if not isinstance(value, list):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid list" % str(value)
+            for order in value:
+                if not isinstance(order, str):
+                    raise self.PhysicsObjectError, \
+                          "%s is not a valid string" % str(value)
+
         if name == 'model':
             if not isinstance(value, Model):
                 raise self.PhysicsObjectError, \
@@ -2471,12 +2621,29 @@
 
         return super(Process, self).set(name, value) # call the mother routine
 
+    def get_squared_order_type(self, order):
+        """ Return what kind of squared order constraint was specified for the
+        order 'order'."""
+
+        if order in self['sqorders_types'].keys():
+            return self['sqorders_types'][order]
+        else:
+            # Default behavior '=' is interpreted as upper bound '<='
+            return '='
+
     def get(self, name):
         """Special get for legs_with_decays"""
         
         if name == 'legs_with_decays':
             self.get_legs_with_decays()
 
+        if name == 'sqorders_types':
+            # We must make sure that there is a type for each sqorder defined
+            for order in self['squared_orders'].keys():
+                if order not in self['sqorders_types']:
+                    # Then assign its type to the default '='
+                    self['sqorders_types'][order]='='
+
         return super(Process, self).get(name) # call the mother routine
 
     def get_sorted_keys(self):
@@ -2486,8 +2653,8 @@
                 'model', 'id', 'required_s_channels', 
                 'forbidden_onsh_s_channels', 'forbidden_s_channels',
                 'forbidden_particles', 'is_decay_chain', 'decay_chains',
-                'legs_with_decays',
-                'perturbation_couplings', 'has_born', 'NLO_mode']
+                'legs_with_decays', 'perturbation_couplings', 'has_born', 
+                'NLO_mode','split_orders']
 
     def nice_string(self, indent=0, print_weighted = True):
         """Returns a nicely formated string about current process
@@ -2519,7 +2686,8 @@
         # Add orders
         if self['orders']:
             mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
-              for key in self['orders'] if print_weighted or key!='WEIGHTED']) + ' '
+              for key in self['orders'] if (print_weighted or key!='WEIGHTED') \
+              and not key in self['squared_orders'].keys()]) + ' '
 
         # Add perturbation_couplings
         if self['perturbation_couplings']:
@@ -2531,9 +2699,11 @@
             mystr = mystr + '] '
 
         # Add squared orders
-        if self['perturbation_couplings'] and self['squared_orders']:
-            mystr = mystr + " ".join([key + '=' + repr(self['squared_orders'][key]) \
-              for key in self['squared_orders'] if print_weighted or key!='WEIGHTED']) + ' ' 
+        if self['squared_orders']:
+            mystr = mystr + " ".join([key + '^2%s%d'%\
+                (self.get_squared_order_type(key),self['squared_orders'][key]) \
+              for key in self['squared_orders'].keys() \
+                                    if print_weighted or key!='WEIGHTED']) + ' '
 
         # Add forbidden s-channels
         if self['forbidden_onsh_s_channels']:
@@ -2786,6 +2956,20 @@
 
     # Helper functions
 
+    def are_negative_orders_present(self):
+        """ Check iteratively that no coupling order constraint include negative
+        values."""
+
+        if any(val<0 for val in self.get('orders').values()+\
+                                           self.get('squared_orders').values()):
+            return True
+        
+        for procdef in self['decay_chains']:
+            if procdef.are_negative_orders_present():
+                return True
+
+        return False
+
     def are_decays_perturbed(self):
         """ Check iteratively that the decayed processes are not perturbed """
         
@@ -2794,6 +2978,14 @@
                 return True
         return False
     
+    def decays_have_squared_orders(self):
+        """ Check iteratively that the decayed processes are not perturbed """
+        
+        for procdef in self['decay_chains']:
+            if procdef['squared_orders']!={} or procdef.decays_have_squared_orders():
+                return True
+        return False
+    
     def get_ninitial(self):
         """Gives number of initial state particles"""
 
@@ -2945,12 +3137,23 @@
         # Ensure that expansion orders are taken into account
         expansion_orders = self.get('model').get('expansion_order')
         orders = self.get('orders')
+        sq_orders = self.get('squared_orders')
         
         tmp = [(k,v) for (k,v) in expansion_orders.items() if 0 < v < 99]
         for (k,v) in tmp:  
             if k in orders:
                 if v < orders[k]:
-                    logger.warning('''The coupling order (%s=%s) specified is larger than the one allowed 
+                    if k in sq_orders.keys() and \
+                                             (sq_orders[k]>v or sq_orders[k]<0):
+                        logger.warning(
+'''The process with the squared coupling order (%s^2%s%s) specified can potentially 
+recieve contributions with powers of the coupling %s larger than the maximal 
+value allowed by the model builder (%s). Hence, MG5_aMC sets the amplitude order
+for that coupling to be this maximal one. '''%(k,self.get('sqorders_types')[k],
+                                             self.get('squared_orders')[k],k,v))
+                    else:
+                        logger.warning(
+'''The coupling order (%s=%s) specified is larger than the one allowed 
              by the model builder. The maximal value allowed is %s. 
              We set the %s order to this value''' % (k,orders[k],v,k))
                     orders[k] = v
@@ -3113,7 +3316,7 @@
 
         return max_order_now, particles, hierarchy
 
-    def nice_string(self, indent=0):
+    def nice_string(self, indent=0, print_weighted=False):
         """Returns a nicely formated string about current process
         content"""
 
@@ -3174,6 +3377,12 @@
                 mystr = mystr + order + ' '
             mystr = mystr + '] '
 
+        if self['squared_orders']:
+            mystr = mystr + " ".join([key + '^2%s%d'%\
+                (self.get_squared_order_type(key),self['squared_orders'][key]) \
+              for key in self['squared_orders'].keys() \
+                                    if print_weighted or key!='WEIGHTED']) + ' '
+
         # Remove last space
         mystr = mystr[:-1]
 
@@ -3192,6 +3401,30 @@
 
         return mystr
 
+    def get_process_with_legs(self, LegList):
+        """ Return a Process object which has the same properties of this 
+            ProcessDefinition but with the specified LegList as legs attribute. 
+            """
+            
+        return Process({\
+            'legs': LegList,
+            'model':self.get('model'),
+            'id': self.get('id'),
+            'orders': self.get('orders'),
+            'sqorders_types': self.get('sqorders_types'),
+            'squared_orders': self.get('squared_orders'),
+            'has_born': self.get('has_born'),
+            'required_s_channels': self.get('required_s_channels'),
+            'forbidden_onsh_s_channels': self.get('forbidden_onsh_s_channels'),            
+            'forbidden_s_channels': self.get('forbidden_s_channels'),
+            'forbidden_particles': self.get('forbidden_particles'),
+            'perturbation_couplings': self.get('perturbation_couplings'),
+            'is_decay_chain': self.get('is_decay_chain'),
+            'overall_orders': self.get('overall_orders'),
+            'split_orders': self.get('split_orders'),
+            'NLO_mode': self.get('NLO_mode')
+            })
+            
     def get_process(self, initial_state_ids, final_state_ids):
         """ Return a Process object which has the same properties of this 
             ProcessDefinition but with the specified given leg ids. """
@@ -3207,19 +3440,9 @@
         for i, fs_id in enumerate(final_state_ids):
             assert fs_id in my_fsids[i]
         
-        return Process({\
-            'legs': LegList(\
+        return self.get_process_with_legs(LegList(\
                [Leg({'id': id, 'state':False}) for id in initial_state_ids] + \
-               [Leg({'id': id, 'state':True}) for id in final_state_ids]),
-            'model':self.get('model'),
-            'id': self.get('id'),
-            'orders': self.get('orders'),
-            'required_s_channels': self.get('required_s_channels'),
-            'forbidden_s_channels': self.get('forbidden_s_channels'),
-            'forbidden_particles': self.get('forbidden_particles'),
-            'perturbation_couplings': self.get('perturbation_couplings'),
-            'is_decay_chain': self.get('is_decay_chain'),
-            'overall_orders': self.get('overall_orders')})
+               [Leg({'id': id, 'state':True}) for id in final_state_ids]))
 
     def __eq__(self, other):
         """Overloading the equality operator, so that only comparison

=== modified file 'madgraph/core/color_algebra.py'
--- madgraph/core/color_algebra.py	2013-11-12 23:52:05 +0000
+++ madgraph/core/color_algebra.py	2014-08-26 07:56:12 +0000
@@ -165,6 +165,36 @@
         return None
 
 #===============================================================================
+# ColorOne
+#===============================================================================
+class ColorOne(ColorObject):
+    """The one of the color object"""
+
+    def __init__(self, *args):
+        """Check for no index"""
+        
+        assert len(args) == 0 , "ColorOne objects must have no index!"
+        
+        super(ColorOne, self).__init__()
+
+    def simplify(self):
+        """"""
+        assert len(self)==0, "There is argument(s) in color object ColorOne."
+        col_str = ColorString()
+        col_str.coeff = fractions.Fraction(1, 1)
+        return ColorFactor([col_str])
+
+
+    def pair_simplify(self, col_obj):
+        """Implement ColorOne product simplification"""
+
+        if any(isinstance(col_obj, c_type) for c_type in [Tr,T,f,d,ColorOne]):
+            col_str = ColorString([col_obj])
+            return ColorFactor([col_str])
+        return None
+
+
+#===============================================================================
 # T
 #===============================================================================
 class T(ColorObject):
@@ -705,6 +735,10 @@
 
         ret_list = [(col_obj.__class__.__name__, tuple(col_obj)) \
                         for col_obj in self]
+
+        if not ret_list and self.coeff:
+            ret_list=[("ColorOne",tuple([]))]
+
         ret_list.sort()
         self.immutable = tuple(ret_list)
 

=== modified file 'madgraph/core/color_amp.py'
--- madgraph/core/color_amp.py	2013-11-12 23:52:05 +0000
+++ madgraph/core/color_amp.py	2014-08-26 07:56:12 +0000
@@ -67,10 +67,12 @@
             min_index, res_dict = self.add_vertex(vertex, diagram, model,
                             repl_dict, res_dict, min_index)
 
-        # Return an empty list if all entries are empty
-        if all([cs == color_algebra.ColorString() \
-                        for cs in res_dict.values()]):
-            res_dict = {}
+        # if the process has no QCD particles
+        # Return a list filled with ColorOne if all entries are empty ColorString()
+        empty_colorstring = color_algebra.ColorString()
+        if all(cs == empty_colorstring for cs in res_dict.values()):
+            res_dict = dict((key, color_algebra.ColorString(
+                               [color_algebra.ColorOne()])) for key in res_dict)
                     
         return res_dict
 
@@ -369,7 +371,6 @@
         indices) associated to my_color_string. Take a list of the external leg
         color octet state indices as an input. Returns only the leading N 
         contribution!"""
-
         # Create a new color factor to allow for simplification
         my_cf = color_algebra.ColorFactor([my_color_string])
 
@@ -453,7 +454,6 @@
             # Rebuild a color string from a CB entry
             col_str = color_algebra.ColorString()
             col_str.from_immutable(col_basis_entry)
-
             for (leg_num, leg_repr) in repr_dict.items():
                 # By default, assign a (0,0) color flow
                 res_dict[leg_num] = [0, 0]
@@ -548,7 +548,7 @@
 
         self.col_matrix_fixed_Nc = {}
         self.inverted_col_matrix = {}
-
+        
         self._col_basis1 = col_basis
         if col_basis2:
             self._col_basis2 = col_basis2
@@ -574,7 +574,6 @@
                     enumerate(sorted(self._col_basis1.keys())):
             for i2, struct2 in \
                     enumerate(sorted(self._col_basis2.keys())):
-
                 # Only scan upper right triangle if symmetric
                 if is_symmetric and i2 < i1:
                     continue
@@ -699,13 +698,18 @@
         appearing in struct1. Assumes internal summed indices are negative."""
 
         # First, determines what is the smallest index appearing in struct1
-        min_index = min(reduce(operator.add,
-                                [list(elem[1]) for elem in struct1])) - 1
+        #list2 = reduce(operator.add,[list(elem[1]) for elem in struct1])
+        list2 = sum((list(elem[1]) for elem in struct1),[])
+        if not list2: 
+            min_index = -1
+        else:
+           min_index = min(list2) - 1
+
         # Second, determines the summed indices in struct2 and create a 
         # replacement dictionary
         repl_dict = {}
-        list2 = reduce(operator.add,
-                       [list(elem[1]) for elem in struct1])
+        #list2 = reduce(operator.add,
+        #               [list(elem[1]) for elem in struct1])
         for summed_index in list(set([i for i in list2 \
                                       if list2.count(i) == 2])):
             repl_dict[summed_index] = min_index
@@ -732,5 +736,8 @@
     @staticmethod
     def lcmm(*args):
         """Return lcm of args."""
-        return reduce(ColorMatrix.lcm, args)
+        if args:
+            return reduce(ColorMatrix.lcm, args)
+        else:
+            return 1
 

=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2014-06-29 15:15:00 +0000
+++ madgraph/core/diagram_generation.py	2014-08-26 07:56:12 +0000
@@ -387,7 +387,6 @@
 
     def __init__(self, argument=None):
         """Allow initialization with Process"""
-
         if isinstance(argument, base_objects.Process):
             super(Amplitude, self).__init__()
             self.set('process', argument)
@@ -507,7 +506,6 @@
         DiagramList by the function. This is controlled by the argument
         returndiag.
         """
-
         process = self.get('process')
         model = process.get('model')
         legs = process.get('legs')
@@ -744,13 +742,20 @@
                 newleg = copy.copy(vert.get('legs').pop(-1))
                 newleg.set('onshell', False)
                 vert.get('legs').append(newleg)
-                
-        # Set diagrams to res
-        self['diagrams'] = res
 
         # Set actual coupling orders for each diagram
         for diagram in res:
             diagram.calculate_orders(model)
+            
+        # Filter the diagrams according to the squared coupling order
+        # constraints and possible the negative one. Remember that OrderName=-n
+        # means that the user wants to include everything up to the N^(n+1)LO
+        # contribution in that order and at most one order can be restricted
+        # in this way. We shall do this only if the diagrams are not asked to
+        # be returned, as it is the case for NLO because it this case the
+        # interference are not necessarily among the diagrams generated here only.
+        if not returndiag and len(res)>0:
+            res = self.apply_squared_order_constraints(res)
 
         # Replace final id=0 vertex if necessary
         if not process.get('is_decay_chain'):
@@ -784,7 +789,10 @@
         self.trim_diagrams(diaglist=res)
 
         # Sort process legs according to leg number
-        self.get('process').get('legs').sort()
+        pertur = 'QCD'
+        if self.get('process')['perturbation_couplings']:
+            pertur = sorted(self.get('process')['perturbation_couplings'])[0]
+        self.get('process').get('legs').sort(pert=pertur)
 
         # Set diagrams to res if not asked to be returned
         if not returndiag:
@@ -793,6 +801,47 @@
         else:
            return not failed_crossing, res
 
+    def apply_squared_order_constraints(self, diag_list):
+        """Applies the user specified squared order constraints on the diagram
+        list in argument."""
+
+        res = copy.copy(diag_list)                  
+
+        # Iterate the filtering since the applying the constraint on one
+        # type of coupling order can impact what the filtering on a previous
+        # one (relevant for the '==' type of constraint).
+        while True:   
+            new_res = res.apply_positive_sq_orders(res, 
+                                          self['process'].get('squared_orders'), 
+                                              self['process']['sqorders_types'])
+            # Exit condition
+            if len(res)==len(new_res):
+                break
+            elif (len(new_res)>len(res)):
+                raise MadGraph5Error(
+                 'Inconsistency in function apply_squared_order_constraints().')
+            # Actualizing the list of diagram for the next iteration
+            res = new_res
+
+        # Now treat the negative squared order constraint (at most one)
+        neg_orders = [(order, value) for order, value in \
+                       self['process'].get('squared_orders').items() if value<0]
+        if len(neg_orders)==1:
+            neg_order, neg_value = neg_orders[0]
+            # Now check any negative order constraint
+            res, target_order = res.apply_negative_sq_order(res, neg_order,\
+                  neg_value, self['process']['sqorders_types'][neg_order])
+            # Substitute the negative value to this positive one so that
+            # the resulting computed constraints appears in the print out
+            # and at the output stage we no longer have to deal with 
+            # negative valued target orders
+            self['process']['squared_orders'][neg_order]=target_order
+        elif len(neg_orders)>1:
+            raise InvalidCmd('At most one negative squared order constraint'+\
+                                   ' can be specified, not %s.'%str(neg_orders))
+
+        return res
+
     def create_diagram(self, vertexlist):
         """ Return a Diagram created from the vertex list. This function can be
             overloaded by daughter classes."""
@@ -908,7 +957,9 @@
         orders subtracted. If coupling_orders is not given, return
         None (which counts as success).
         WEIGHTED is a special order, which corresponds to the sum of
-        order hierarchys for the couplings."""
+        order hierarchies for the couplings.
+        We ignore negative constraints as these cannot be taken into
+        account on the fly but only after generation."""
 
         if not coupling_orders:
             return None
@@ -922,7 +973,8 @@
             for coupling in inter.get('orders').keys():
                 # Note that we don't consider a missing coupling as a
                 # constraint
-                if coupling in present_couplings:
+                if coupling in present_couplings and \
+                                                 present_couplings[coupling]>=0:
                     # Reduce the number of couplings that are left
                     present_couplings[coupling] -= \
                              inter.get('orders')[coupling]
@@ -930,7 +982,8 @@
                         # We have too many couplings of this type
                         return False
             # Now check for WEIGHTED, i.e. the sum of coupling hierarchy values
-            if 'WEIGHTED' in present_couplings:
+            if 'WEIGHTED' in present_couplings and \
+                                               present_couplings['WEIGHTED']>=0:
                 weight = sum([model.get('order_hierarchy')[c]*n for \
                               (c,n) in inter.get('orders').items()])
                 present_couplings['WEIGHTED'] -= weight
@@ -1124,7 +1177,7 @@
     def trim_diagrams(self, decay_ids=[], diaglist=None):
         """Reduce the number of legs and vertices used in memory.
         When called by a diagram generation initiated by LoopAmplitude, 
-        this function should no trim the diagrams in the attribute 'diagrams'
+        this function should not trim the diagrams in the attribute 'diagrams'
         but rather a given list in the 'diaglist' argument."""
 
         legs = []
@@ -1164,7 +1217,6 @@
                     diagram.get('vertices')[ivx] = vertices[index]
                 except ValueError:
                     vertices.append(vertex)
-        
 
 #===============================================================================
 # AmplitudeList
@@ -1627,10 +1679,14 @@
                                  process_definition.get('perturbation_couplings'),
                               'squared_orders': \
                                  process_definition.get('squared_orders'),
+                              'sqorders_types': \
+                                 process_definition.get('sqorders_types'),                              
                               'overall_orders': \
                                  process_definition.get('overall_orders'),
                               'has_born': \
-                                 process_definition.get('has_born')
+                                 process_definition.get('has_born'),
+                              'split_orders': \
+                                 process_definition.get('split_orders')                                
                                  })
                 fast_proc = \
                           array.array('i',[leg.get('id') for leg in legs])
@@ -1854,6 +1910,12 @@
                                  process_definition.get('required_s_channels'),
                               'forbidden_onsh_s_channels': \
                                  process_definition.get('forbidden_onsh_s_channels'),
+                              'sqorders_types': \
+                                 process_definition.get('sqorders_types'),
+                              'squared_orders': \
+                                 process_definition.get('squared_orders'),
+                              'split_orders': \
+                                 process_definition.get('split_orders'), 
                               'forbidden_s_channels': \
                                  process_definition.get('forbidden_s_channels'),
                               'forbidden_particles': \
@@ -1861,7 +1923,9 @@
                               'is_decay_chain': \
                                  process_definition.get('is_decay_chain'),
                               'overall_orders': \
-                                 process_definition.get('overall_orders')})
+                                 process_definition.get('overall_orders'),
+                              'split_orders': \
+                                 process_definition.get('split_orders')})
 
                     # Check for couplings with given expansion orders
                     process.check_expansion_orders()

=== modified file 'madgraph/core/drawing.py'
--- madgraph/core/drawing.py	2013-11-13 02:21:10 +0000
+++ madgraph/core/drawing.py	2014-08-26 07:56:12 +0000
@@ -92,8 +92,6 @@
         assert isinstance(vertex, VertexPoint), 'The begin point should be a ' + \
                  'Vertex_Point object'
 
-        assert vertex is not self.end
-
         self.begin = vertex
         vertex.add_line(self)
         return
@@ -104,8 +102,6 @@
         assert isinstance(vertex, VertexPoint), 'The end point should be a ' + \
                  'Vertex_Point object'
 
-        assert vertex is not self.begin
-                 
         self.end = vertex
         vertex.add_line(self)
         return
@@ -1797,6 +1793,8 @@
 
         # check if we need curved loop particles
         curved_for_loop = False
+        circled_for_loop = False
+        
         if isinstance(diagram, LoopFeynmanDiagram):
             # If only 2 particle in the loop require that those lines are
             # curved
@@ -1804,10 +1802,17 @@
                 curved_for_loop = True
                 self.curved_part_start = (0, 0)
         
+            # for tadpole DOES NOT CRASH BUT STILL NEED FIXING
+            elif len([l for l in diagram.lineList if l.loop_line]) == 1:
+                circled_for_loop = True
+                self.curved_part_start = (0, 0)
+
         # drawing the particles
-        for line in diagram.lineList:    
-            if not curved_for_loop or not line.loop_line:
+        for line in diagram.lineList:
+            if (not curved_for_loop and not circled_for_loop) or not line.loop_line:
                 self.draw_line(line)
+            elif circled_for_loop:
+                self.draw_circled_line(line)
             else:
                 self.draw_curved_line(line)
                 
@@ -1856,6 +1861,34 @@
             number = line.number
             self.associate_number(line, number)
 
+    # To draw tadpole
+    def draw_circled_line(self, line):
+        """Draw the line information.
+        First, call the method associate the line type [draw_circled_XXXXXX]
+        Then finalize line representation by adding his name."""
+
+        # if 4 point (or more interaction at the beginning/end of the loop
+        # need to reduce curvature
+        if len(line.begin.lines) > 3 or len(line.end.lines) > 3 :
+            cercle = False
+        else:
+            cercle = True
+
+        # Find the type line of the particle [straight, wavy, ...]
+        line_type = line.get_info('line')
+        # Call the routine associate to this type [self.draw_straight, ...]
+        if hasattr(self, 'draw_circled_' + line_type):
+            getattr(self, 'draw_circled_' + line_type)(line, cercle)
+        else:
+            self.draw_circled_straight(line, reduce)
+            
+        # Finalize the line representation with adding the name of the particle
+        name = line.get_name()
+        self.associate_name(line, name)
+        
+        #store begin for helping future curving
+        self.curved_part_start = (line.begin.pos_x, line.begin.pos_y*1.2)
+
     def draw_curved_line(self, line):
         """Draw the line information.
         First, call the method associate the line type [draw_curved_XXXXXX]

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2013-11-29 07:28:53 +0000
+++ madgraph/core/helas_objects.py	2014-08-26 07:56:12 +0000
@@ -89,6 +89,7 @@
         else:
             sorted_tags = sorted([cls(d, model, ninitial) for d in \
                                       amplitude.get('diagrams')])
+
         # Do not use this for loop diagrams as for now the HelasMultiProcess
         # always contain only exactly one loop amplitude.
         if sorted_tags and not isinstance(amplitude, \
@@ -898,12 +899,11 @@
         part = self.get('particle')
         self.set('particle', self.get('antiparticle'))
         self.set('antiparticle', part)
-
+    
     def is_anticommutating_ghost(self):
-        """ Return True if the particle of this wavefunction is a ghost charged
-        under SU(3)"""
-        return self.get('particle').get('ghost') and self.get('color')!=1
-
+        """ Return True if the particle of this wavefunction is a ghost"""
+        return self.get('particle').get('ghost')
+    
     def is_fermion(self):
         return self.get('spin') % 2 == 0
 
@@ -1006,9 +1006,11 @@
         array_rep.append(self['color_key'])
         # Also need to specify if it is a loop wf
         array_rep.append(int(self['is_loop']))
+        
         # Finally, the mother numbers
         array_rep.extend([mother['number'] for \
                           mother in self['mothers']])
+        
         return array_rep
 
     def get_pdg_code(self):
@@ -1039,9 +1041,9 @@
             if particles[1].get_pdg_code() != particles[2].get_pdg_code() \
                    and self.get('pdg_code') == \
                    particles[1].get_anti_pdg_code()\
-                   and self.get('coupling')[0] != '-':
+                   and not self.get('coupling')[0].startswith('-'):
                 # We need a minus sign in front of the coupling
-                self.set('coupling', ['-' + c for c in self.get('coupling')])
+                self.set('coupling', ['-%s'%c for c in self.get('coupling')])
 
     def set_octet_majorana_coupling_sign(self):
         """For octet Majorana fermions, need an extra minus sign in
@@ -1052,12 +1054,9 @@
         if self.get('color') == 8 and \
                self.get_spin_state_number() == -2 and \
                self.get('self_antipart') and \
-               [m.get('color') for m in self.get('mothers')] == [8, 8]:
-            if hasattr(self, 'octet_majorana_flip'):
-                return
-            self.octet_majorana_flip = True
-            self.set('coupling', ['-%s' % c 
-                                              for c in self.get('coupling')])
+               [m.get('color') for m in self.get('mothers')] == [8, 8] and \
+               not self.get('coupling')[0].startswith('-'):
+            self.set('coupling', ['-%s' % c for c in self.get('coupling')])
         
     def set_state_and_particle(self, model):
         """Set incoming/outgoing state according to mother states and
@@ -1260,11 +1259,36 @@
                                 diagram_wavefunctions.insert(i, new_wf)
                                 break
                         else:
-                            diagram_wavefunctions.append(new_wf)
+                            # For loop processes, care is needed since
+                            # some loop wavefunctions in the diag_wfs might have
+                            # the new_wf in their mother, so we want to place 
+                            # new_wf as early as possible in the list.
+                            # We first look if any mother of the wavefunction
+                            # we want to add appears in the diagram_wavefunctions
+                            # list. If it doesn't, max_mother_index is -1.
+                            # If it does, then max_mother_index is the maximum
+                            # index in diagram_wavefunctions of those of the 
+                            # mothers present in this list.
+                            max_mother_index = max([-1]+
+                                [diagram_wavefunctions.index(wf) for wf in 
+                                        mothers if wf in diagram_wavefunctions])
+                            
+                            # We want to insert this new_wf as early as 
+                            # possible in the diagram_wavefunctions list so that
+                            # we are guaranteed that it will be placed *before*
+                            # wavefunctions that have new_wf as a mother.
+                            # We therefore place it at max_mother_index+1.
+                            if max_mother_index<len(diagram_wavefunctions)-1:
+                                new_wf.set('number',diagram_wavefunctions[
+                                              max_mother_index+1].get('number'))
+                            for wf in diagram_wavefunctions[max_mother_index+1:]:
+                                wf.set('number',wf.get('number')+1)
+                            diagram_wavefunctions.insert(max_mother_index+1,
+                                                                         new_wf)
 
             # Set new mothers
             new_wf.set('mothers', mothers)
-
+                    
             # Now flip flow or sign
             if flip_flow:
                 # Flip fermion flow
@@ -1283,7 +1307,7 @@
                 new_wf_number = new_wf.get('number')
                 new_wf = wavefunctions[wavefunctions.index(new_wf)]
                 diagram_wf_numbers = [w.get('number') for w in \
-                                      diagram_wavefunctions]
+                                                          diagram_wavefunctions]
                 index = diagram_wf_numbers.index(new_wf_number)
                 diagram_wavefunctions.pop(index)
                 # We need to decrease the wf number for later
@@ -1293,13 +1317,27 @@
                         wf.set('number', wf.get('number') - 1)
                 # Since we reuse the old wavefunction, reset wf_number
                 wf_number = wf_number - 1
+
                 # Need to replace wavefunction in number_to_wavefunctions
-                # (in case this wavefunction is in another of the dicts)
+                # (in case this wavefunction is in another of the dicts) 
                 for n_to_wf_dict in number_to_wavefunctions:
                     if new_wf in n_to_wf_dict.values():
                         for key in n_to_wf_dict.keys():
                             if n_to_wf_dict[key] == new_wf:
-                                n_to_wf_dict[key] = new_wf
+                                n_to_wf_dict[key] = new_wf          
+ 
+                if self.get('is_loop'):
+                    # fix a bug for the g g > go go g [virt=QCD]
+                    # when there is a wf which is replaced, we need to propagate
+                    # the change in all wavefunction of that diagrams which could
+                    # have this replaced wavefunction in their mothers. This
+                    # plays the role of the 'number_to_wavefunction' dictionary
+                    # used for tree level.
+                    for wf in diagram_wavefunctions:
+                        for i,mother_wf in enumerate(wf.get('mothers')):
+                            if mother_wf.get('number')==new_wf_number:
+                                wf.get('mothers')[i]=new_wf
+
             except ValueError:
                 pass
 
@@ -1464,7 +1502,15 @@
                     
         #fixed argument
         for i, coup in enumerate(self.get_with_flow('coupling')):
-            output['coup%d'%i] = coup
+            # We do not include the - sign in front of the coupling of loop
+            # wavefunctions (only the loop ones, the tree ones are treated normally)
+            # in the non optimized output because this sign was already applied to
+            # the coupling passed in argument when calling the loop amplitude.
+            if not OptimizedOutput and self.get('is_loop'):
+                output['coup%d'%i] = coup[1:] if coup.startswith('-') else coup  
+            else:
+                output['coup%d'%i] = coup
+              
         output['out'] = self.get('me_id') - flip
         output['M'] = self.get('mass')
         output['W'] = self.get('width')
@@ -1543,10 +1589,9 @@
 
         # Sort according to spin and flow direction
         res.sort()
-#        if not self['is_loop']:
         res.append(self.get_spin_state_number())
         res.append(self.find_outgoing_number())
-#        else:
+
         if self['is_loop']:
             res.append(self.get_loop_index())
             if not self.get('mothers'):
@@ -2075,6 +2120,7 @@
                                    wf_number,
                                    force_flip_flow=False,
                                    number_to_wavefunctions=[]):
+
         """Check for clashing fermion flow (N(incoming) !=
         N(outgoing)). If found, we need to trace back through the
         mother structure (only looking at fermions), until we find a
@@ -2267,6 +2313,76 @@
 
         return conjugates
 
+    
+    def check_wavefunction_numbers_order(self, applyChanges=False, raiseError=True):
+        """ This function only serves as an internal consistency check to 
+        make sure that when setting the 'wavefunctions' attribute of the
+        diagram, their order is consistent, in the sense that all mothers 
+        of any given wavefunction appear before that wavefunction.
+        This function returns True if there was no change and the original 
+        wavefunction list was consistent and False otherwise.
+        The option 'applyChanges' controls whether the function should substitute
+        the original list (self) with the new corrected one. For now, this function
+        is only used for self-consistency checks and the changes are not applied."""
+    
+        if len(self)<2:
+            return True