Merge lp:~maddevelopers/mg5amcnlo/nlo_reweighting into lp:~maddevelopers/mg5amcnlo/2.3.4

[Olivier Mattelaer]

Hi,

This branch

1) includes change from Rikkert corresponding to the new format of the reweight information
2) include a new flag in the NLO run_card allowing to keep/discard the reweight information
3) include the possibility to perform a proper NLO reweighting

This implementation follows closely the one done in the unleashed_reweighting (accepted) branch for the strategy/interface. The details on how to use the branch is on this page:
https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Reweight

Basic validations are presented in this talk:
https://cp3.irmp.ucl.ac.be/projects/madgraph/raw-attachment/wiki/Reweight/15_12_11_ReWeighting_EFT.pdf

I would propose that Rikkert focus on the change in
fks_singular/ handling_lhe_events as well as lhe_parser.py
While Valentin can focus on the computation of the loop with f2py.

Cheers and Thanks,

Olivier

PS: A couple of unittest/acceptance tests are failing but they are all already failing in 2.3.4.
I'm going to track the responsible(s) of those failing tests in 2.3.4 and will re-merge 2.3.4 when those will be fixed in that branch.

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

347. By Olivier Mattelaer on 2016-01-27

adding missing file

[Rikkert Frederix (frederix)]

Hi Olivier,

It doesn't work. f2py tries to link LHAPDF, even though I didn't specify the code to use LHAPDF and I get a crash that it doesn't fine the -lLHAPDF library.

cheers,
Rik

review: Needs Fixing

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

Yes this mode need lhapdf.
I do not want to play with internal pdf.

I will add this in the documentation and add the code to check if lhapdf is present or not before proposing to the user
to use the module.

Cheers and thanks,

Olivier

> On Feb 3, 2016, at 13:49, Rikkert Frederix <email address hidden> wrote:
>
> Review: Needs Fixing
>
> Hi Olivier,
>
> It doesn't work. f2py tries to link LHAPDF, even though I didn't specify the code to use LHAPDF and I get a crash that it doesn't fine the -lLHAPDF library.
>
> cheers,
> Rik
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/nlo_reweighting/+merge/283769
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/2.3.4.

[Rikkert Frederix (frederix)]

Hi Olivier,

I double checked. I don't know if your solution helps. In fact, I do have LHAPDF installed, it's in the PATH (as well as that I specified the path to lhapdf-config in the mg5_configuration.txt file), and I can even generate the events with LHAPDF. Still, it doesn't work for the reweight module.

Cheers,
Rik

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

348. By Olivier Mattelaer on 2016-02-04

add details for NLO reweighting in the reweight_card

349. By Olivier Mattelaer on 2016-02-06

merge with 2.3.4 (code not working anymore)

350. By Olivier Mattelaer on 2016-02-07

fixing linking problem with lhapdf6

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

Could you retry now?
I have merged with the latest 2.3.4 and fix the way to link to lhapdf.
I'm still trying to understand better how to correctly link to lhapdf but I hope that you should have something working now.

Thanks,

Olivier

[Rikkert Frederix (frederix)]

Hi Olivier,

It still isn't working. Although it doesn't seem to be related to LHAPDF. I get the following compilation error:

Command "launch auto " interrupted with error:
MadGraph5Error : A compilation Error occurs when trying to compile /tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/rw_mevirt/SubProcesses/P1_csx_epve.
 The compilation fails with the following output message:
     cd ..; make OLP_static; cd -
     make[1]: Entering directory `/tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/rw_mevirt/SubProcesses'
     ar rcs libMadLoop.a MadLoopParamReader.o MadLoopCommons.o P0_udx_epve/loop_matrix.o P1_csx_epve/loop_matrix.o P0_udx_epve/improve_ps.o P1_csx_epve/improve_ps.o P0_udx_epve/born_matrix.o P1_csx_epve/born_matrix.o P0_udx_epve/CT_interface.o P1_csx_epve/CT_interface.o P0_udx_epve/loop_num.o P1_csx_epve/loop_num.o P0_udx_epve/polynomial.o P1_csx_epve/polynomial.o P0_udx_epve/mp_compute_loop_coefs.o P1_csx_epve/mp_compute_loop_coefs.o P0_udx_epve/TIR_interface.o P1_csx_epve/TIR_interface.o
     mv libMadLoop.a ../lib/libMadLoop.a
     make[1]: Leaving directory `/tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/rw_mevirt/SubProcesses'
     /tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/rw_mevirt/SubProcesses/P1_csx_epve
     ln -s ../../lib/libMadLoop.a
     ln: failed to create symbolic link ‘./libMadLoop.a’: File exists
     make: *** [../OLP_static] Error 1

 Please try to fix this compilations issue and retry.
 Help might be found at https://answers.launchpad.net/madgraph5.
 If you think that this is a bug, you can report this at https://bugs.launchpad.net/madgraph5
Please report this bug on https://bugs.launchpad.net/madgraph5
More information is found in '/tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/run_01_tag_1_debug.log'.
Please attach this file to your report.

Cheers,
Rik

review: Needs Fixing

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

351. By Olivier Mattelaer on 2016-02-09

make possible mulitple linking to the library even when ln -s crashes

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

Thanks a lot for this beta testing, this is really useful to have.
I guess that this should be fixed now.

Cheers,

Olivier

> On Feb 9, 2016, at 09:56, Rikkert Frederix <email address hidden> wrote:
>
> Review: Needs Fixing
>
> Hi Olivier,
>
> It still isn't working. Although it doesn't seem to be related to LHAPDF. I get the following compilation error:
>
> Command "launch auto " interrupted with error:
> MadGraph5Error : A compilation Error occurs when trying to compile /tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/rw_mevirt/SubProcesses/P1_csx_epve.
> The compilation fails with the following output message:
> cd ..; make OLP_static; cd -
> make[1]: Entering directory `/tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/rw_mevirt/SubProcesses'
> ar rcs libMadLoop.a MadLoopParamReader.o MadLoopCommons.o P0_udx_epve/loop_matrix.o P1_csx_epve/loop_matrix.o P0_udx_epve/improve_ps.o P1_csx_epve/improve_ps.o P0_udx_epve/born_matrix.o P1_csx_epve/born_matrix.o P0_udx_epve/CT_interface.o P1_csx_epve/CT_interface.o P0_udx_epve/loop_num.o P1_csx_epve/loop_num.o P0_udx_epve/polynomial.o P1_csx_epve/polynomial.o P0_udx_epve/mp_compute_loop_coefs.o P1_csx_epve/mp_compute_loop_coefs.o P0_udx_epve/TIR_interface.o P1_csx_epve/TIR_interface.o
> mv libMadLoop.a ../lib/libMadLoop.a
> make[1]: Leaving directory `/tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/rw_mevirt/SubProcesses'
> /tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/rw_mevirt/SubProcesses/P1_csx_epve
> ln -s ../../lib/libMadLoop.a
> ln: failed to create symbolic link ‘./libMadLoop.a’: File exists
> make: *** [../OLP_static] Error 1
>
> Please try to fix this compilations issue and retry.
> Help might be found at https://answers.launchpad.net/madgraph5.
> If you think that this is a bug, you can report this at https://bugs.launchpad.net/madgraph5
> Please report this bug on https://bugs.launchpad.net/madgraph5
> More information is found in '/tuph/t30/space/frederix/MG5_aMC/nlo_reweighting/PROCNLO_loop_sm_0/run_01_tag_1_debug.log'.
> Please attach this file to your report.
>
>
> Cheers,
> Rik
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/nlo_reweighting/+merge/283769
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/2.3.4.

[Rikkert Frederix (frederix)]

Hi Olivier,

I really think you should do some more testing before asking for a review. Now the code runs, but the reweighted cross sections is orders of magnitude away from the correct value.

generate p p > e+ ve [QCD] ; output
launch --> "reweight ON' and 'keep_rwgt_info --> True'

(It gives the warning that the cards are identical, but it runs anyway due to DEBUG mode).

gives that the reweighted cross section is ~4.3 pb, while the original is ~1.0*10^4 pb.

Cheers,
Rik

review: Needs Fixing

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

Actually, such tests was working when I proposed the merging.
This bug was introduced with the latest merge of 2.3.4.
Funny that it did not show up for the EFT test that I did after this merging.
Since the only big change in 2.3.4 is linked to the lhapdf, I guess that's the problem are related to that.

I will fix that, will add a couple of acceptance tests on this functionality to prevent similar future problem and will rerun a battery of tests.

I will let you know when this is ready.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

352. By Olivier Mattelaer on 2016-02-14

reintroducing PDF factor and improve the combination of the weights

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,

So we already did some testing together on skype about this branch and identified an issue arising when using PDF reweighting together with your reweighting. I guess you fixed it by now.

So here are my comments on the functionalities:

a) If I understood correctly, the run_card option

False = keep_rwgt_info ! keep the information for reweighting (huge impact on size file)

does more than just changing what is stored in the intermediate lhe, but also dictates what type of reweighting will be performed. I think that this should be emphasized more, maybe with a comment above this parameter in the card.
Also, my understanding is that keeping these additional weights in the .lhe will increase its filesize by more or less a constant factor independent of the process (actually this factor would become smaller for more complicated processes since one must store more momenta).
I tried the comparison on 10k events for p p > e+ ve and I got a factor a bit less than 3 (both zipped and unzipped). So I think the parenthesis after this factor should read

(leads to an .lhe event file about ~3 times larger)

This would scare people much less than 'a huge factor'.

b) When using GoSam as an OLP, the true reweighting crashes (expectedly so) [but now for the kamikaze reweighting, which is good]. It would be nice if you could handle this earlier an in a nicer way, because it made python crash (!) not even an exception thrown! I guess it is ok an nice to keep the kamikaze reweighting working in this case.

c) With lhapdf I still have the issue of having '-lc++' in make_opts. Did you push your fix? It works fine in 2.3.4

d) generating the square matrix element for reweighting
   ->
   generating the squared matrix element for reweighting

e) It would be nice that when it says 'Running Reweight' it reads instead

   'INFO: Running reweighting (exactly NLO accurate)' for the true reweighting

   and for the kamikaze one, it reminds (once again):

   "INFO: Running reweighting (approximate)"
   "WARNING: To perform exact NLO accurate reweighting, please turn on the run_card parameter 'keep_rwgt_info'"

   Also I would rephrase your existing warning :

Warning: To have NLO accurate reweighting, you need to specify in the
+# run_card (at the time of the NLO sample generation)
+# keep_rwgt_info on "True". If it is on False a LO type of
+# reweighting will be performed (not assuring NLO accuracy)

->

Warning: To insure NLO accuracy in the reweighting procedure, you must set the
+# run_card parameter 'keep_rwgt_info' to '.true.'(when doing event generation).
+# If it is set to '.false.', then and approximate reweighting based on
+# LO matrix element only will be performed (hence losing formal NLO accuracy)

f) I tried to reweight p p > e+ ve [QCD] by setting the WWidth to 10.0 (kamikaze reweighting) I got:

INFO: new cross-section is : 2180.81 pb (indicative error: 14.8007 pb)
INFO: Original cross-section: 10825.316 +- 40.826559 pb
INFO: Computed cross-section:
INFO: 1 : 2180.80629797 +- 14.8006705613 pb

And when doing it the other way around (generating with a width of 10 Gev and reweighting to...

review: Needs Fixing

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

353. By Olivier Mattelaer on 2016-02-15

merge with latest 2.3.4

354. By Olivier Mattelaer on 2016-02-19

- allow the possibility to reweight wihout the virt-trick handling
- allow to run more than one reweighting (VT/standard/LO) at the same time
- adding better support for LHAPDF (automatic downloading of the sample)
- fixing a couple of glitch in the reading of the weight in the shower analysis

355. By Olivier Mattelaer on 2016-02-19

merge with 2.3.4 (336)

[Olivier Mattelaer (olivier-mattelaer)]

@Rik,

So the identity tests was failing due to the lhapdf problem that I face with the merge with Marco branch. For the test, I remove some call to lhapdf and forget to put those back. I have add an acceptance test to check that part of the code.

So you can go back to testing.

@Valentin,

What an extensive review. Thanks so much for this.

> a) (leads to an .lhe event file about ~3 times larger)

Done ( I have put 5)

> I think that this should be emphasized more, maybe with a comment above this parameter in the card.

Don't think this is necessary in the card, if reweighting is on, a warning is set automatically anyway.
Actually, I have modified the behavior such that the default when setting reweighitng on ON is to set that parameter on True automatically.

> b) When using GoSam as an OLP.

Ok Done.

> c) With lhapdf I still have the issue of having '-lc++' in make_opts. Did you
> push your fix? It works fine in 2.3.4

Should be fixed now actually. I did not made the push so far indeed

> d) generating the square matrix element for reweighting
> ->
> generating the squared matrix element for reweighting
>
done

> e) It would be nice that when it says 'Running Reweight' it reads

Well this is a bit tricky. I never know for sure if an exact NLO reweighitng is going to be used since the user has always the possibility to put it on LO mode.
So I have put "INFO: Running reweighting (approximate)" when I'm sure that the kamikaze method is used.

> Also I would rephrase your existing warning :

DONE

> f) Anyway, I think that it might be nice to have a warning message if your
> reweighted relative error becomes say 100 times larger than the one of the
> original cross-section, do you have one?

Well the reweighting can be 100% correct and the error be 100 times larger than the original one.
(if the cross-section is also multiplied by exactly 100) so such simple warning does not work in general.

> g) ==== POWHEG (OR POSSIBLY HERWIG) WILL USE LHAPDFv6 ====

This one is funny, this is probably a misspell from Rik or Marco.
This come from one fortran file.

> terminate called after throwing an instance of 'LHAPDF::ReadError'
> what(): Empty PDF file name given to Info::load
> Abort trap: 6
>
> What is this about :) ??

Crash in the f2py are like that. they make python to crash in a quite bad way...
If we decide to push this f2py idea, we will need a way to recover from those kind of crash but this is ok for the moment.

>If lhapdf is indeed needed for true reweighting, then it must be handled upstream, in the same way
>as when the user tries to turn on PDF reweighting and built-in PDFs are selected.

I have put a switch to automatically do LO reweighitng in that case, with the appropriate warning.

> h)
> Should it already follow the strategy we discussed so as to easily derive the
> HEPMC weight name afterwards? Anyway, this branch also update the shower
> driver so as to keep these weights in the

Ok I have change type= to name= like in the susmentioned banch.
I will not put the diff in the attribute.

> (Detail: also maybe pick a shorter name than 'mg_reweight_1' given that it is
> bound to appear fo each event.)

Done...

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

356. By Olivier Mattelaer on 2016-02-19

apply Valentin modification

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,

On Fri, Feb 19, 2016 at 9:12 AM, Olivier Mattelaer <
<email address hidden>> wrote:

> @Rik,
>
> So the identity tests was failing due to the lhapdf problem that I face
> with the merge with Marco branch. For the test, I remove some call to
> lhapdf and forget to put those back. I have add an acceptance test to check
> that part of the code.
>
> So you can go back to testing.
>
>
> @Valentin,
>
> What an extensive review. Thanks so much for this.
>
> > a) (leads to an .lhe event file about ~3 times larger)
>
> Done ( I have put 5)
>

Well ok, but 3 was already rounding up aove the factor I had found :).

>
> > I think that this should be emphasized more, maybe with a comment above
> this parameter in the card.
>
> Don't think this is necessary in the card, if reweighting is on, a warning
> is set automatically anyway.
> Actually, I have modified the behavior such that the default when setting
> reweighitng on ON is to set that parameter on True automatically.
>

Ok then. Also because when doing kamikaze the subsequent warning makes it
clear how to do the exact reweighting.

But the problem is that if I set it to false then it will still be
automatically be reassigned to True if I try to launch again. The user will
certainly don't expect his run_card.dat to be changed on his behalf when it
does a new run.
So I'm not sure I like this behavior....

> > b) When using GoSam as an OLP.
>
> Ok Done.
>
>
> > c) With lhapdf I still have the issue of having '-lc++' in make_opts.
> Did you
> > push your fix? It works fine in 2.3.4
>
> Should be fixed now actually. I did not made the push so far indeed

Ok it works now for the kamikaze reweighting with lhapdf reweighting.
However, I tried a kamikaze reweighting *with the default reweighting card*
(so launch only) on p p > e+ ve with lhapdf and pdf reweighting (with lhaid
= 11000) and I got

INFO: new cross-section is : 10827.4 pb (indicative error: 42.5828 pb)
INFO: Original cross-section: 11125.181 +- 42.582827 pb
INFO: Computed cross-section:
INFO: 1 : 10827.4488 +- 42.582827 pb

The cross-sections aren't identical, so this is purely an effect of the
kamikaze reweighting right? I guess I'm missing something because I'd
expect that even in this case you divide and multiply by the same thing, so
you should leave the original weight untouched. Since kamikaze reweighting
is so fast, I would consider to *always* use this "central reweighting" so
as to be able to automatically provide the user with an assessment of how
reliable he can expect the kamikaze reweighting to be for this case.

Unfortunately, it still crashes at compilation time when doing the exact
reweighting with lhapdf reweighting. The error being:

[...]
ld: library not found for -lLHAPDF
collect2: error: ld returned 1 exit status
error: Command "/usr/local/bin/master/gfortran -Wall -g -Wall -g -undefined
dynamic_lookup -bundle
/var/folders/4c/cbnnj8zs0l728myp8kt89j500000gn/T/tmpW9KIcx/var/folders/4c/cbnnj8zs0l728myp8kt89j500000gn/T/tmpW9KIcx/src.macosx-10.8-x86_64-2.7/rwgt2pymodule.o
/var/folders/4c/cbnnj8zs0l728myp8kt89j500000gn/T/tmpW9KIcx/var/folders/4c/cbnnj8zs0l728myp8kt89j500000gn/T/tmpW9KIcx/src.macosx-10.8-x86_64...

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

357. By Olivier Mattelaer on 2016-02-20

try to fix Valentin problem

358. By Olivier Mattelaer on 2016-02-20

remove golem + other fix related to Valentin review

359. By Olivier Mattelaer on 2016-02-21

remove pjfry as well

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

Thanks a lot for those very useful comments.
And for the time on Skype. I'm so happy that it is working on your computer now.

 > But the problem is that if I set it to false then it will still be
> automatically be reassigned to True if I try to launch again. The user will
> certainly don't expect his run_card.dat to be changed on his behalf when it
> does a new run.
> So I'm not sure I like this behavior....

Ok you have a point, I will not do that then.

I have then apply the following.
When the user quit the edition of the card interface.
I check the consistency between the run_card and the reweight_card. and if the user request explicitly the NLO reweighting type, then I force the keep_rwgt_info to True.

I have modified the default reweight_card to always request explicitly the NLO accurate reweighting.
Such that by default, I have the same behavior.

> The cross-sections aren't identical, so this is purely an effect of the
> kamikaze reweighting right?

No you also have an statistical effect. Even if you set all the reweighting factor to be 1
(which is the case of your kamikaze in this case) you will face that the sum of the weights is not equal to the cross-section. This is an effect partly due to the of the negative weights.

Since kamikaze reweighting
> is so fast, I would consider to *always* use this "central reweighting" so
> as to be able to automatically provide the user with an assessment of how
> reliable he can expect the kamikaze reweighting to be for this case.

This is equivalent to simply compute the sum of the weights (up to normalisation) which is ... so much faster than running such kamikaze reweighting. Printing the sum of the weights for the sample is actually an interesting idea.

> Unfortunately, it still crashes at compilation time when doing the exact
> reweighting with lhapdf reweighting. The error being:

Ok we will need to settle that via Skype.

> Ok, that would do I guess. So how would a user do a reweight in LO mode
> with an event that *includes the extra weight*?

You have to put in your reweight_card the line
change mode LO
The fact that the line
"change mode NLO" is now present by default in the reweight_card should clear your confusion on this.

> Also, speaking of this command, its help is a bit minimal, i.e.:
>
> PROCNLO_loop_sm_3>help reweight
> Allow to reweight the events generated with a new choices of model
> parameter.

OK DONE

> And its auto-completion is broken (with a debug statement left in by the
> way).

OK DONE

> Sure this is why I stressed *relative* error. All my comments above
> referred to a *relative* error, so please reconsider my comments with that
> in mind.
> I think this goes along the way of giving more info to the user on the
> reliability of the reweighting

I think that this is too much. They have the information printed it should be the user task to check if it looks correct or not. Actually I do not trust that much the error quoted but this is the best estimator that I have found.

> Did you check this out?
> http://stackoverflow.com/questions/10929162/how-can-i-catch-a-seg-fault-while-
> importing-an-f2py-module

Just did and th...

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

360. By Olivier Mattelaer on 2016-02-21

reduce the logger to have nicer screen info

361. By Olivier Mattelaer on 2016-02-21

use the mur value associate to the couter-events for the loop computation

362. By Olivier Mattelaer on 2016-02-21

merge with latest 2.3.4

363. By Olivier Mattelaer on 2016-02-22

small fix

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier

> Hi Valentin,
>
> Thanks a lot for those very useful comments.
> And for the time on Skype. I'm so happy that it is working on your computer
> now.
>
>
> > But the problem is that if I set it to false then it will still be
> > automatically be reassigned to True if I try to launch again. The user will
> > certainly don't expect his run_card.dat to be changed on his behalf when it
> > does a new run.
> > So I'm not sure I like this behavior....
>
> Ok you have a point, I will not do that then.
>
> I have then apply the following.
> When the user quit the edition of the card interface.
> I check the consistency between the run_card and the reweight_card. and if the
> user request explicitly the NLO reweighting type, then I force the
> keep_rwgt_info to True.
>
> I have modified the default reweight_card to always request explicitly the NLO
> accurate reweighting.
> Such that by default, I have the same behavior.

I like the default reweight card much better this way.
However, I would change this:

        - For LO sample only LO mode is allowed.
        - NLO runs both NLO_basic and NLO_VT reweighting methods.
        - Read the online manual to check the difference between those modes.

to

        - When reweighting an .lhe event file with LO accuracy, only the LO mode is allowed.
        - When the .lhe file reweighted is at NLO accuracy, then all other modes are allowed;
          read about their definition @ https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Reweight

and also:

change mode NLO # Do nothing on LO sample.

->

change mode NLO # No reweighting will be performed if the reweighted .lhe file is at LO accuracy.

Anyway, now I correctly have this happening:

INFO: NLO reweighting is on ON. Automatically set keep_rwgt_info to True
INFO: modify parameter keep_rwgt_info of the run_card.dat to .true.

which is nice but there is still a problem with your dependence on lhapdf with the exact NLO reweighting. You must really check for this at the beginning and take action correspondingly.

I don't know what strategy you want to adopt here, but I would suggest to automatically switch to kamikaze reweighting if lhapdf is not used, but with a big warning. Or 'crash' the run and tell the user to either turn on the kamikaze reweighting or to turn on lhapdf.

>
> > The cross-sections aren't identical, so this is purely an effect of the
> > kamikaze reweighting right?
>
> No you also have an statistical effect. Even if you set all the reweighting
> factor to be 1
> (which is the case of your kamikaze in this case) you will face that the sum
> of the weights is not equal to the cross-section. This is an effect partly due
> to the of the negative weights.

Ok, I understand now.

> Since kamikaze reweighting
> > is so fast, I would consider to *always* use this "central reweighting" so
> > as to be able to automatically provide the user with an assessment of how
> > reliable he can expect the kamikaze reweighting to be for this case.
>
> This is equivalent to simply compute the sum of the weights (up to
> normalisation) which is ... so much faster than running such kamikaze
> reweighting. Printing the sum of t...

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

> I like the default reweight card much better this way.
> However, I would change this:

Not really happy for those proposal.
So I have made another version which should be clearer (and stress that both NLO_VT and NLO_basic are both NLO accurate).

> which is nice but there is still a problem with your dependence on lhapdf with
> the exact NLO reweighting. You must really check for this at the beginning and
> take action correspondingly.

Ok this is handle in a nicer way and the pdf set are automatically downloaded now.

> What about doing a trial import with this:
>
> subprocess.call("python -c 'import module_or_modules_in_question'")
>
> allowing you to nicely handle a crash there? I know this seems clunky but I
> think it is ok and really much better than a python crash. I mean MG5aMC
> crashing with no message whatsoever is really scary for the user and looks
> bad; we should really avoid it at all cost if possible!

I can do that, but this is obviously very limited since all C++/fortran error at run time will not be catched (like Lhapdf do not find the pdfset/... even if python finds it/...)
So I think that the nicer method would be to run the full reweighting in a separate subprocess.
I will do that.

> You have to call this function:
>
> MadLoopInitialize.init_MadLoop in madevent_interface.py
>
> yourself. But it really is a one-liner, so it should really be easy.
> It is best to initialize it first so that you have a trustworthy accuracy
> estimate for the first couple of events you reweight.
>

I tried but it seems ... not working:
root: Error while executing make in /Users/omatt/Documents/eclipse/nlo_reweighting/PROCNLO_loop_sm_85/rw_mevirt/SubProcesses/P1_csx_epve
root: Failed at running the process in /Users/omatt/Documents/eclipse/nlo_reweighting/PROCNLO_loop_sm_85/rw_mevirt/SubProcesses/P1_csx_epve.

> I don't think so (or I don't understand what you meant). The ratio of loop ME
> for the original and new hypothesis can for sure depend on the renormalization
> scale, so it is not enough that the numerator and denominator have the same
> scale, it must *also* be the same scale as used in the original event I think

DONE

Thanks a lot,

Olivier

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

> I like the default reweight card much better this way.
> However, I would change this:

Not really happy for those proposal.
So I have made another version which should be clearer (and stress that both NLO_VT and NLO_basic are both NLO accurate).

> which is nice but there is still a problem with your dependence on lhapdf with
> the exact NLO reweighting. You must really check for this at the beginning and
> take action correspondingly.

Ok this is handle in a nicer way and the pdf set are automatically downloaded now.

> What about doing a trial import with this:
>
> subprocess.call("python -c 'import module_or_modules_in_question'")
>
> allowing you to nicely handle a crash there? I know this seems clunky but I
> think it is ok and really much better than a python crash. I mean MG5aMC
> crashing with no message whatsoever is really scary for the user and looks
> bad; we should really avoid it at all cost if possible!

I can do that, but this is obviously very limited since all C++/fortran error at run time will not be catched (like Lhapdf do not find the pdfset/... even if python finds it/...)
So I think that the nicer method would be to run the full reweighting in a separate subprocess.
I will do that.

> You have to call this function:
>
> MadLoopInitialize.init_MadLoop in madevent_interface.py
>
> yourself. But it really is a one-liner, so it should really be easy.
> It is best to initialize it first so that you have a trustworthy accuracy
> estimate for the first couple of events you reweight.
>

I tried but it seems ... not working:
root: Error while executing make in /Users/omatt/Documents/eclipse/nlo_reweighting/PROCNLO_loop_sm_85/rw_mevirt/SubProcesses/P1_csx_epve
root: Failed at running the process in /Users/omatt/Documents/eclipse/nlo_reweighting/PROCNLO_loop_sm_85/rw_mevirt/SubProcesses/P1_csx_epve.

> I don't think so (or I don't understand what you meant). The ratio of loop ME
> for the original and new hypothesis can for sure depend on the renormalization
> scale, so it is not enough that the numerator and denominator have the same
> scale, it must *also* be the same scale as used in the original event I think

DONE

Thanks a lot,

Olivier

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

364. By Olivier Mattelaer on 2016-02-23

fixing automatic download of pdf sets + update reweight card

365. By Olivier Mattelaer on 2016-02-23

allow to run MadEvent/aMC@NLO in bypassing RunWeb if run from the interface

366. By Olivier Mattelaer on 2016-02-23

ensure that scan can be use with './bin/madevent launch -f'

367. By Olivier Mattelaer on 2016-02-23

run reweighting in a separate process + print the sum of event weight for NLO run

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

Just to say that
1) I run now the reweight in a separate process. Such that the main code do not crash.
2) I print the sum of the weights in the final report.
3) The initialisation is still failing (but the code still works nicely) I left the associated error message in the log such that you can look at it if you want.

Cheers,

Olivier

[Rikkert Frederix (frederix)]

Hi Olivier,

I wanted to go back to testing this branch, but still I can't get it working at NLO. Even with the default setup (only turning the REWEIGHT=ON), the code crashes when starting the reweighting with a IndexError: list index out of range (line 1120 in common_run_interface.py).

B.t.w., what's the difference between NLO_basic and NLO_VT ? The wiki-page is not very enlightening on this subject.

Cheers,
Rik

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

> IndexError: list index out of range (line 1120 in common_run_interface.py).

Did you have the latest version of the code?
In my computer this lines are:

                if any(t in line for t in ['INFO:', 'WARNING:', 'CRITICAL:', 'ERROR:', 'root:']) and \
                   not ‘***********’ in line: #<<<<<<<<<<<<<line 1120
                        print line[:-1].replace(‘INFO’, 'REWEIGTH')
                elif __debug__ and line:
                    logger.debug(line[:-1])

This is just impossible that line 1120 (or 1119) raise an IndexError.
Line 1121: can not raise an index error (potentially a TypeError but not an IndexError.

So my best guess is that we do not have the same version. I’m at version 367

> B.t.w., what's the difference between NLO_basic and NLO_VT ? The wiki-page is not very enlightening on this subject.

They differ on how they reweight some of the contribution
- the born
- the approximate virtual
- the virtual
They have the same treatment for the various counter-term.

NLO_basic reweighs
- the born by the ratio of the borns
- the approximate virtual by the ratio of the borns
- the virtual by the ratio of the virtual

NLO_VT reweighs those three contributions by the ratio of the sum of the Born and the virtual (B+V).

The interest of those two is that this can enlighten the sensitivity to the virt-trick if they differ from each other.

Cheers,

Olivier

> On Feb 25, 2016, at 08:56, Rikkert Frederix <email address hidden> wrote:
>
> Hi Olivier,
>
> I wanted to go back to testing this branch, but still I can't get it working at NLO. Even with the default setup (only turning the REWEIGHT=ON), the code crashes when starting the reweighting with a IndexError: list index out of range (line 1120 in common_run_interface.py).
>
> B.t.w., what's the difference between NLO_basic and NLO_VT ? The wiki-page is not very enlightening on this subject.
>
> Cheers,
> Rik
>
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/nlo_reweighting/+merge/283769
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/2.3.4.

[Rikkert Frederix (frederix)]

Hi Olivier,

Sorry. I had a typo. It's line 1130. I'm using revision 367.

I don't think there is any freedom in doing the reweighting. There can only be one correct way, I think. I believe it is:

1. use the ratio of the Borns for the Born and the counter terms.
2. use the ratio of the virtuals for the virtual and the approximated virtual.

Cheers,
Rik

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

> Sorry. I had a typo. It's line 1130. I’m using revision 367.

Ok, thanks, I will secure that part (This is just to have a nice html output and removing some fake entry).

> I don't think there is any freedom in doing the reweighting. There can only be one correct way, I think.

That’s just not True. Let’s take a simple example (no physics in mind here)

\int(b+v) and you want to reweight to get \int(B+V)
you can easily check that (at least) two reweighting procedure are valid
1)
b term reweighted by B/b
v term reweighted by V/v

2)
b term reweighted by (B+V)/(b+v)
v term reweighted by (B+V)/(b+v)

Both are theoritically correct (in infinite statistic limit).
Only a detail study on how fast they converge (which is likely to be function dependent) can determine which of those two method is the most appropriate.
My feeling is that the second method is more efficient but I actually do not know for sure.

> 1. use the ratio of the Borns for the Born and the counter terms.
> 2. use the ratio of the virtual for the virtual and the approximated virtual.

To take the kind of example above you are trying to reweight:
\int{b+a}+ \int{v-a} = \int(b+v-a+a) and we want to get \int(B+V)
so you can reweight “a” by any factor that you like since the integral does not depend on it at the end.
(So you actually have an infinite number of possible reweighting method)

if we take “a” proportional to “b” then for me the two most natural method are the one I have implemented.
b term reweighted by B/b
v term reweighted by V/v
a term reweighted by A/a=B/b
and
b term reweighted by (B+V)/(b+v)
v term reweighted by (B+V)/(b+v)
a term reweighted by (B+V)/(b+v)

The one that you propose is
b term reweighted by B/b
v term reweighted by V/v
a term reweighted by V/v

The advantage of the first method (which is actually NLO_basic) is that you have \int{b+a} which is re-weighted to \int{B+A} and \int{v+a} which is re-weighted to \int{V+A}.
The second advantage of that method is that this is the most natural one (to my point of view).

The advantage of the two other methods (which can be called NLO_VT and NLO_Rik) is to have “v” and “a” are reweighted by the same factor. Which is arguably an interesting property to avoid the contribution of \int{v-a} to blow up after the reweighting procedure.
I do not think this is that interesting to have both of them implemented in the code (so three NLO accurate one). I would argue that having NLO_VT is slightly more interesting since it applies a global factor to all of them.
Which correspond to the second reweighting method of my first example. Since I have a feeling that it can have some nice mathematical property, I think it is legitimate to prefer that one over the other one.

Cheers,

Olivier

> On Feb 25, 2016, at 12:18, Rikkert Frederix <email address hidden> wrote:
>
> Hi Olivier,
>
> Sorry. I had a typo. It's line 1130. I'm using revision 367.
>
> I don't think there is any freedom in doing the reweighting. There can only be one correct way, I think. I believe it is:
>
> 1. use the ratio of the Borns for the Born and the counter terms....

[Rikkert Frederix (frederix)]

Hi Olivier,

Ok. Now I understand. The problem is that this is not what's in the code... In the code the weight that correspond to the virtual is actually 'v-a'. In other words, 'v' and 'a' can currently not be reweighted independently. Which let me (incorrectly) to think that there is only a single method.

This means that NLO_VT and NLO_Rik are correct, but NLO_Basic does currently not give the correct cross section.

I'm pretty sure that the statistical difference between NLO_VT and NLO_Rik would be pretty small, so I suggest to simply use NLO_VT.

Cheers,
Rik

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

368. By Olivier Mattelaer on 2016-02-25

secure the html treatment

369. By Olivier Mattelaer on 2016-02-26

fixing a problem forbiding to run with scan parameter

370. By Olivier Mattelaer on 2016-02-27

merge with 2.3.4 (including spin2/ninja branch)

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,

So I tried again the very default run, that is

generate p p > e+ ve [NLO]
output
launch

Then I turned off the shower and turned on the reweighting (without turning on LHAPDF).
As expected the code turned keep_rwgt_info to True automatically, which is nice, i.e.

NLO reweighting is on ON. Automatically set keep_rwgt_info to True

but then it crashed, probably because it needs lhadpf?
This time around and thanks to the sandboxing I suggested and that you worked hard to implement, the crash is nice:

[...]
DEBUG: Traceback (most recent call last):
DEBUG: File "/Users/valentin/Documents/Work/MG5/nlo_reweighting/PROCNLO_loop_sm_18/bin/internal/amcatnlo_run_interface.py", line 4816, in <module>
DEBUG: except MadEventAlreadyRunning, error:
DEBUG: NameError: name 'MadEventAlreadyRunning' is not defined
ERROR: Reweighting failed
decay_events -from_cards

Then I went on and turned LHAPDF but got the same error as above. So it looks like you have an issue with the fact that you now run the reweight separately but there is of course still the 'RunWeb' tagger file present, so it crashes. How come you didn't run into the same issue?

Even with the kamikaze reweighting (mode LO), the issue above occurs.
So I need it to be fixed in order to continue the review.

Cheers

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

371. By Olivier Mattelaer on 2016-02-28

remove the wrong raise of Error, forbid the orig weight to be written in the file

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

This has nothing to do with RunWeb. I implemented a hidden method to run one single command bypassing all interactions with RunWeb.

Looks like the error catching was using a undefined error in amcatnlo_runinterface (while it was define in the madevent one).
and that make the code to fail to crash correctly.

I fix it but you will still face the real crash… (that I have no idea what it is in your case)
I doubt that this is linked to lhapdf since you also have it for kamikaze…

Cheers,

Olivier

> On Feb 28, 2016, at 03:29, Valentin Hirschi <email address hidden> wrote:
>
> Hi Olivier,
>
> So I tried again the very default run, that is
>
> generate p p > e+ ve [NLO]
> output
> launch
>
> Then I turned off the shower and turned on the reweighting (without turning on LHAPDF).
> As expected the code turned keep_rwgt_info to True automatically, which is nice, i.e.
>
> NLO reweighting is on ON. Automatically set keep_rwgt_info to True
>
> but then it crashed, probably because it needs lhadpf?
> This time around and thanks to the sandboxing I suggested and that you worked hard to implement, the crash is nice:
>
> [...]
> DEBUG: Traceback (most recent call last):
> DEBUG: File "/Users/valentin/Documents/Work/MG5/nlo_reweighting/PROCNLO_loop_sm_18/bin/internal/amcatnlo_run_interface.py", line 4816, in <module>
> DEBUG: except MadEventAlreadyRunning, error:
> DEBUG: NameError: name 'MadEventAlreadyRunning' is not defined
> ERROR: Reweighting failed
> decay_events -from_cards
>
> Then I went on and turned LHAPDF but got the same error as above. So it looks like you have an issue with the fact that you now run the reweight separately but there is of course still the 'RunWeb' tagger file present, so it crashes. How come you didn't run into the same issue?
>
> Even with the kamikaze reweighting (mode LO), the issue above occurs.
> So I need it to be fixed in order to continue the review.
>
> Cheers
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/nlo_reweighting/+merge/283769
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/2.3.4.

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

372. By Olivier Mattelaer on 2016-02-28

ensure that the order of the weights is always the same

373. By Olivier Mattelaer on 2016-02-28

try to fix Valentin problem

374. By Olivier Mattelaer on 2016-02-29

remove a print statement

375. By Olivier Mattelaer on 2016-03-01

refactorise directory creation (allow for 2 model @NLO) + remove NLO_Basic + add NLO_tree (for loop-improved mode)

376. By Olivier Mattelaer on 2016-03-01

fix a parsing error in presence of second process definition

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

>I'm pretty sure that the statistical difference between NLO_VT and NLO_Rik would be pretty small, > so I suggest to simply use NLO_VT.

I tried to fix the code to have the correct NLO_basic but I failed.
Therefore I implemented your suggestion (and renamed it to NLO).

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

377. By Olivier Mattelaer on 2016-03-01

fix a problem for Loop Induced Reweighting

378. By Olivier Mattelaer on 2016-03-01

secure last fix

379. By Rikkert Frederix on 2016-03-02

fixed the problem with ngluon

380. By Olivier Mattelaer on 2016-03-02

Automatic download of lhapdf set for SysCalc

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,

Let me continue this review

a) My OCD at work again, sorry:
    >> NLO reweighting is on ON. Automatically set keep_rwgt_info to True --> detail: you could remove the 'on' ON.
    >> remove leading space in ' When doing NLO reweighting, MG5aMC canno'
    >> generating the square matrix element for reweighting -> generating squared matrix element for reweighting
    >> REWEIGHT: starts to compute weight for events with the following modification to the param_card:
REWEIGHT: change model loop_qcd_qed_sm
DEBUG: change process p p > e+ ve
DEBUG: change process p p > e+ ve j

why only showing the first line in 'INFO' and the rest in 'DEBUG' ? imho, the whole set of lines should be in INFO.

b) Debug left at:
CRITICAL: w_orig, w_new, w_new/w_orig = 0.393281380656 0.393281380656 1.0 [reweight_interface.py at line 813]

c) I tried to run the default PY8 analysis for p p > e+ ve [QCD] with your kamikaze reweighting 'ON', i.e.
by setting 'py8an_HwU_pp_lvl.o HwU.o' in the analysis parameter of the shower card, and it fails with:

terminate called after throwing an instance of 'std::runtime_error'
  what(): Tried to read a file which does not start with the LesHouchesEvents tag.
Abort trap: 6

This seems to come from the fact that the events are gzipped again before the shower step.
If I unzip them by hand and try again, it runs further but the shower+analysis stops again with the error message:

 WW(1) = 0. Stopping

So maybe the shower analysis must be updated so as to work with your reweighting (check with the improved_scale_pdf_variation branch, maybe this is already done there).
In any case, please make sure that you can run a default shower analysis on top of your reweighting.

d) In the card, I think it would be better to use another example than:

# SPECIFY A PATH OR USE THE SET COMMAND LIKE THIS:
 25 # set sminputs 1 130 # modify 1/alpha_EW

because with the above, you will change the W mass with the default EW input scheme and this is typically a dangerous thing to do.

e) I tried to reweight p p > e+ ve [QCD] in less conventional ways, with:

change mode NLO # Define type of Reweighting. For LO sample this command
                     # has no effect since only "LO" mode is allowed.
change model loop_qcd_qed_sm
change process u d~ > e+ ve QCD=99 QED=99 [virt=QCD QED]
change process c s~ > e+ ve QCD=99 QED=99 [virt=QCD QED] --add
launch

But it crashed with:

DEBUG: Command "reweight run_03 -from_cards" interrupted with error:
DEBUG: InvalidCmd : No particle qcd=99 in model

I think it is important that any valid process syntax gets correctly parsed by the reweighting module.
Then I tried with

change model loop_qcd_qed_sm
change process u d~ > e+ ve
change process c s~ > e+ ve --add
launch

DEBUG: Command "reweight run_03 -from_cards" interrupted with error:
DEBUG: KeyError : ((2, 21), (-11, 1, 12))

I guess I am missing some processes, but it would be nice if the error was more explicit!
So I tried again with something like :

change model loop_qcd_qed_sm
change process p p > e+ ve
change process p p > e+ ve j --add
launch

But this time crashed with:

CRITICAL: w_orig, w_new, w_new/w_orig = 4.3346035354e-05 4.3346...

review: Needs Fixing

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

a) why only showing the first line in 'INFO' and the rest in 'DEBUG' ? imho, the whole set of lines should be in INFO.
The sandboxing, allow very low control on what is printed, so if a line do not start by INFO, it pass in DEBUG.

> b) Debug left at:
> CRITICAL: w_orig, w_new, w_new/w_orig = 0.393281380656 0.393281380656 1.0 [reweight_interface.py at line 813]

Sorry, I was studying a weird behavior today and forget to remove that print statement.

> c) I tried to run the default PY8 analysis for p p > e+ ve [QCD] with your kamikaze reweighting 'ON', i.e.
> by setting 'py8an_HwU_pp_lvl.o HwU.o' in the analysis parameter of the shower card, and it fails with:

It runs fine with HW6. Never succeed to run PY8 analysis… so I will keep it like that.
Paolo told me that the parsing should have be fine….

> d) In the card, I think it would be better to use another example than:

This is fine like that

> e) I think it is important that any valid process syntax gets correctly parsed by the reweighting module.

The problem is that a lower case is applied… The parser is the one of madgraph.
Will check why.
But the first fix is to allow such syntax. (why the coupling order should be case sensitive)

> I guess I am missing some processes, but it would be nice if the error was more explicit!

That is fine like that

> f) The following syntax instruction in the reweight card is wrong:

Thanks. fixed

> g) Can I ask you a favor? Could you add your fortran wrapper that includes all relevant MadLoop functions as well as the f2py meta-codes in the standalone MadLoop output as well? (and make sure there is a target in the makefile to f2py it). I find it very cool and it would be nice to have it available in the standalone output (doesn’t cost anything).

I guess you know the answer.

> h) Related to the point above, the reason why the intialization failed is because you have substituted check_sa.f !!!! So of course now it has not ‘main' program statement and when you try to compile it, you get :

Do not see what you are talking about.

Cheers,

Olivier

> On Mar 3, 2016, at 00:15, Valentin Hirschi <email address hidden> wrote:
>
> Review: Needs Fixing
>
> Hi Olivier,
>
> Let me continue this review
>
> a) My OCD at work again, sorry:
>>> NLO reweighting is on ON. Automatically set keep_rwgt_info to True --> detail: you could remove the 'on' ON.
>>> remove leading space in ' When doing NLO reweighting, MG5aMC canno'
>>> generating the square matrix element for reweighting -> generating squared matrix element for reweighting
>>> REWEIGHT: starts to compute weight for events with the following modification to the param_card:
> REWEIGHT: change model loop_qcd_qed_sm
> DEBUG: change process p p > e+ ve
> DEBUG: change process p p > e+ ve j
>
> why only showing the first line in 'INFO' and the rest in 'DEBUG' ? imho, the whole set of lines should be in INFO.
>
> b) Debug left at:
> CRITICAL: w_orig, w_new, w_new/w_orig = 0.393281380656 0.393281380656 1.0 [reweight_interface.py at line 813]
>
> c) I tried to run the default PY8 analysis for p p > e+ ve [QCD] with you...

[Valentin Hirschi (valentin-hirschi)]

Him
> Hi,
>
> a) why only showing the first line in 'INFO' and the rest in 'DEBUG' ? imho,
> the whole set of lines should be in INFO.
> The sandboxing, allow very low control on what is printed, so if a line do not
> start by INFO, it pass in DEBUG.

Ok, then maybe just print these lines all one by one so that they are prefixed by INFO and therefore correctly handled by the sandboxing? Anyway, up to you.

>
> > b) Debug left at:
> > CRITICAL: w_orig, w_new, w_new/w_orig = 0.393281380656 0.393281380656 1.0
> [reweight_interface.py at line 813]
>
> Sorry, I was studying a weird behavior today and forget to remove that print
> statement.
>
> > c) I tried to run the default PY8 analysis for p p > e+ ve [QCD] with your
> kamikaze reweighting 'ON', i.e.
> > by setting 'py8an_HwU_pp_lvl.o HwU.o' in the analysis parameter of the
> shower card, and it fails with:
>
>
> It runs fine with HW6. Never succeed to run PY8 analysis… so I will keep it
> like that.
> Paolo told me that the parsing should have be fine….

I confirm that an HW6 run or PYTHIA6PT run goes through fine, so this issue is really a problem of the parsing in the py8 driver/analysis which is not really part of this branch, so it will have to be fixed somewhere else I agree. But I think that it is important that it will; hopefully this will come for free when merging with the improved_scale/pdf branch.

>
> > d) In the card, I think it would be better to use another example than:
>
> This is fine like that
>
> > e) I think it is important that any valid process syntax gets correctly
> parsed by the reweighting module.
>
> The problem is that a lower case is applied… The parser is the one of
> madgraph.
> Will check why.
> But the first fix is to allow such syntax. (why the coupling order should be
> case sensitive)

I don't see why they shouldn't be case sensitive. Anyway, there should be no need for a lower_case() here no?

>
> > I guess I am missing some processes, but it would be nice if the error was
> more explicit!
>
> That is fine like that

Ok.
You didn't reply to the other issues I found when changing the model, processes and directory to reuse. So when doing exact NLO reweighting how are these process definition used? How do you reweight differently the Born and virtual in this case?

>
> > f) The following syntax instruction in the reweight card is wrong:
>
> Thanks. fixed
>
>
> > g) Can I ask you a favor? Could you add your fortran wrapper that includes
> all relevant MadLoop functions as well as the f2py meta-codes in the
> standalone MadLoop output as well? (and make sure there is a target in the
> makefile to f2py it). I find it very cool and it would be nice to have it
> available in the standalone output (doesn’t cost anything).
>
> I guess you know the answer.

The answer is yes because it's just a simple 'cp command' to add to the standalone output? Great, thanks.

>
> > h) Related to the point above, the reason why the intialization failed is
> because you have substituted check_sa.f !!!! So of course now it has not
> ‘main' program statement and when you try to compile it, you get :
>
> Do not see what you are talking about.

We sorted that out o...

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

381. By Olivier Mattelaer on 2016-03-03

implement some of the comment of Valentin

[Valentin Hirschi (valentin-hirschi)]

Here is just a summary of the status of my review since it is likely to go idle for a couple of weeks.

a) As far as the use of 'change process' for NLO reweighting is concerned, it seems both powerful but also quite tricky and not fully explained on the wiki when used in the NLO context.
If this is intended to remain a 'hidden' functionality for now, I suggest that whenever a 'change process' is issued for NLO reweighting, a message warns the user that this is an alpha functionality that requires expertise.
Anyway, many of the issues reported in my previous mail have been fixed over our skype session, so as far as 'change process' at NLO goes, things are fine like this, with maybe the message above if you like it too.

b) The 'change_rwgt_dir' option didn't seem to work.

c) The loop initialization is maybe not ideal so if you skip it is fine, but then make sure to call the subroutine
 CALL ML5_0_FORCE_STABILITY_CHECK(.True.)
to make sure that you get a stability estimate even during the first calls which will be during initialization for which there is no stability test by default. The above would enforce it.

d) It is important not to use check_sa.f as being both the fortran driver *and* the f2py wrapper. This should be a separate file, say 'f2py_wrapper.f' that you use to make the matrix<n>.so.

e) Also, it would be nice to add the file 'f2py_wrapper.f' mentioned above to the default MadLoop standalone output. If you don't want to do it, then let me know when you'll have this file and I'll do it myself, no problem.

f) In the MadLoop makefile there is a new rule which seems to be designed to make a static library for just one Pxxx process, i.e.

$(MADLOOP_LIB): $(OLP_PROCESS)
ifeq ($(ROOT), ../..)
 cd ..;make ../lib/libMadLoop.$(libext);cd -
else
 ar rcs libMadLoop.$(libext) $(OLP_PROCESS)
 mv libMadLoop.$(libext) $(MADLOOP_LIB)
endif

I don't think that it makes much sense, also because it conflicts with the rule of the same name in MadLoop_makefile_definitions. There is no real reason not to generate the library for all SubProcesses at once by doing 'make OLP_static' in the Subprocesses directory.
So I think that f2py should be linked directly to the static MadLoop library that includes all the processes and can be compiled in the usual way with the existing rule in MadLoop_makefile_definitions or at least this rule for a MadLoop library of only one subprocess must be renamed so as not to create a conflict with the existing one.

That's all, sorry for being such a picky reviewer; I promise to accept the branch right away once the above points are discussed/addressed ;) !

Cheers

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

382. By Olivier Mattelaer on 2016-03-03

Fix the rwgt_dir option for NLO reweighting

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

> a) As far as the use of 'change process' for NLO reweighting is concerned, it
> seems both powerful but also quite tricky and not fully explained on the wiki
> when used in the NLO context.
> If this is intended to remain a 'hidden' functionality for now, I suggest that
> whenever a 'change process' is issued for NLO reweighting, a message warns the
> user that this is an alpha functionality that requires expertise.
> Anyway, many of the issues reported in my previous mail have been fixed over
> our skype session, so as far as 'change process' at NLO goes, things are fine
> like this, with maybe the message above if you like it too.

Will keep it like that, simple redefinition works. And this options is fully fonctional for LO type of reweighting.

> b) The 'change_rwgt_dir' option didn't seem to work.

Actually, I left that unchecked on purpose, In the idea to do it if someone (=CMS) want to have such mode. But since you point it here. I have taken the time to make it work.
Pretty sure no one will use it but ok.

> c) The loop initialization is maybe not ideal so if you skip it is fine, but
> then make sure to call the subroutine
> CALL ML5_0_FORCE_STABILITY_CHECK(.True.)
> to make sure that you get a stability estimate even during the first calls
> which will be during initialization for which there is no stability test by
> default. The above would enforce it.

Done.
Not fully sure where I have to put such line.
So please check (and if this is not correct please indicate me where you want such line)

> d) It is important not to use check_sa.f as being both the fortran driver
> *and* the f2py wrapper. This should be a separate file, say 'f2py_wrapper.f'
> that you use to make the matrix<n>.so.

Well I do not see the point here. check_sa is suppose to be a working template to present one way to call the loop computation nothing more right? So what is the point.

> e) Also, it would be nice to add the file 'f2py_wrapper.f' mentioned above to
> the default MadLoop standalone output. If you don't want to do it, then let me
> know when you'll have this file and I'll do it myself, no problem.

ok I put the file, but then you will need to modify the makefile since now it is pointless to have a file that you can not compile...

> f) In the MadLoop makefile there is a new rule which seems to be designed to
> make a static library for just one Pxxx process, i.e.

Here the best is that you perform the modification that you want directly.
It will be more efficient and less error prone. Especially since you have to change it anyway because of the above file.

> That's all, sorry for being such a picky reviewer; I promise to accept the
> branch right away once the above points are discussed/addressed ;) !

Thanks to you.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

383. By Olivier Mattelaer on 2016-03-03

last modif request by Valentin

384. By Valentin Hirschi on 2016-03-03

1. Fixed the makefile targets for generating the MadLoop library. It is still not very elegant, but at least now there
   is some logic and no conflict of rules with the one $(MADLOOPLIB) used in aMC@NLO runs.
2. No longer modifying check_sa.f in the reweighting but doing all with f2py_wrapper.py
3. Added a f2py driver for MadLoop standalone output which compiles itself when simply doing '.\check_sa.py', so neat!!

385. By Valentin Hirschi on 2016-03-03

1. Fixed a type in 'check_sa.py.inc' template file

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,
> Hi,
>
>
> > a) As far as the use of 'change process' for NLO reweighting is concerned,
> it
> > seems both powerful but also quite tricky and not fully explained on the
> wiki
> > when used in the NLO context.
> > If this is intended to remain a 'hidden' functionality for now, I suggest
> that
> > whenever a 'change process' is issued for NLO reweighting, a message warns
> the
> > user that this is an alpha functionality that requires expertise.
> > Anyway, many of the issues reported in my previous mail have been fixed over
> > our skype session, so as far as 'change process' at NLO goes, things are
> fine
> > like this, with maybe the message above if you like it too.
>
> Will keep it like that, simple redefinition works. And this options is fully
> fonctional for LO type of reweighting.

I think that if 'change process' is only partially working at NLO, there must be warning of some sort, or ideally MG5aMC should detect and crash on the cases which are not supported (see my comments for point b. below).
But ok, I agree that this is quite of a marginal use, so it is up to you.

>
> > b) The 'change_rwgt_dir' option didn't seem to work.
>
> Actually, I left that unchecked on purpose, In the idea to do it if someone
> (=CMS) want to have such mode. But since you point it here. I have taken the
> time to make it work.
> Pretty sure no one will use it but ok.

I am strongly against having (advertised, i.e. non-hidden) features of the code which are partially functional. This is dangerous on a longer term, and even for options rarely used like 'change_rwgt_dir'.
My philosophy is that MG5aMC should always either support an option/functionality correctly or stop immediately saying that this option is not supported.
Any in-between behavior is dangerous and should be prohibited. So I am happy that you settled for fixing it (rather than forbidding it).
And I don't agree that nobody will use it, one can save a lot of time this way and I am really happy that there is a way to control what, when and where stuff should be reused. Thanks.

> > c) The loop initialization is maybe not ideal so if you skip it is fine, but
> > then make sure to call the subroutine
> > CALL ML5_0_FORCE_STABILITY_CHECK(.True.)
> > to make sure that you get a stability estimate even during the first calls
> > which will be during initialization for which there is no stability test by
> > default. The above would enforce it.
>
> Done.
> Not fully sure where I have to put such line.
> So please check (and if this is not correct please indicate me where you want
> such line)

With my fix of the points below, the initialization does go through now. But I understand why you might still want not to use it.
Anyway, I checked that where you placed this call is indeed correct, so it is ok.

>
> > d) It is important not to use check_sa.f as being both the fortran driver
> > *and* the f2py wrapper. This should be a separate file, say 'f2py_wrapper.f'
> > that you use to make the matrix<n>.so.
>
> Well I do not see the point here. check_sa is suppose to be a working template
> to present one way to call the loop computation nothing more right? So what is
> the point...

review: Approve

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

386. By Valentin Hirschi on 2016-03-04

1. Fixed estheatics in check_sa.py

387. By Olivier Mattelaer on 2016-03-05

better comparison of energy-momentum conservation in ML. Fix NLO_tree (loop_improved)

388. By Olivier Mattelaer on 2016-03-05

merge with latest 2.3.4

389. By Olivier Mattelaer on 2016-03-08

fix a problem with alpha_s and with the merging of two fks particle+ add security if particle are too soft/colinear (then use born ratio)

390. By Olivier Mattelaer on 2016-03-09

allow to specify the tag of the reweighting label

391. By Olivier Mattelaer on 2016-03-10

improve check to 0 momenta in MadLoop (re-introduce some absolute lower band)

392. By Olivier Mattelaer on 2016-03-11

1) Change the handling of make_opts to fix the -lstdc++ problem related to f2py
2) fix various problem with the possibility to choose the name of the weight_tag in the reweithing
3) Change some lhapdf interface to be compatible with lhapdf5

[Rikkert Frederix (frederix)]

Hi Olivier,

I approve this branch. Just, quickly fix the order for the linking in the makefile as we found out by e-mail.

Cheers,
Rik

review: Approve

[lp:~maddevelopers/mg5amcnlo/nlo_reweighting]

393. By Valentin Hirschi on 2016-03-16

1. Fixed the issue in testIO_ProcOutputIOTests. Didn't update the IOTest because you have to review your changes in the ref files.

=== modified file 'Template/Common/Cards/reweight_card_default.dat'
--- Template/Common/Cards/reweight_card_default.dat	2015-08-11 12:18:25 +0000
+++ Template/Common/Cards/reweight_card_default.dat	2016-03-16 17:59:03 +0000
@@ -1,13 +1,9 @@
 #*************************************************************************
-#                       Reweight Module                                  *
+#                          Reweight Module                               *
+#                Matrix-Element reweighting at LO/NLO                    *
+#  Mattelaer Olivier                                    arxiv:xxxx.xxxx  *  
 #*************************************************************************
-#
-# Reweighting by Leading Order (tree-level / loop-induced) matrix
-# elements
-#
-# Caution: this is not (in general) NLO accurate, even when starting
-#          with NLO accurate events.
-#
+#          
 # Note: 
 #   1) the value of alpha_s will be used from the event so the value in
 #      the param_card is not taken into account.
@@ -17,15 +13,32 @@
 #      width leads to very inaccurate result). In such a case, separate 
 #      event generation runs are needed.
 #
+#************************************************************************
+# ENTER YOUR COMMANDS BELOW. 
+#************************************************************************
+
+change mode NLO    # Define type of Reweighting. For LO sample this command
+                   # has no effect since only "LO" mode is allowed.
+                   
+launch
+#  SPECIFY A PATH OR USE THE SET COMMAND LIKE THIS:
+# set sminputs 1 130 # modify 1/alpha_EW
+
+
+
+#************************************************************************
+# Manual:      https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Reweight
+#************************************************************************
 #
 # Example of (standard) code for the computation of two weights:
 #
-# launch                    ! tag to start the computat'n of the first weight
+# launch                    ! tag to start the computation of the first weight
 #  set BLOCKNAME ID VALUE   ! rule to modify the current param_card
 #  set BLOCKNAME ID VALUE   ! rule to modify the current param_card
 # launch                    ! start to compute a second weight
 #  /home/Cards/param_card_2.dat ! you can also enter a path to a valid card
 #
+# Note: The command to specify the parameter are AFTER the associated "launch" 
 #
 # Possible options: 
 #   You can enter one of the following lines to customize the reweighting 
@@ -35,17 +48,22 @@
 #        with. In this case you need to provide the path to a correct 
 #        param_card for the new model; you cannot modify the original one 
 #        with the 'set' command.
-# change process DEF [-add]: change the process by which you reweight. 
+# change process DEF [--add]: change the process by which you reweight. 
 #        The initial and final state particles of the new process should 
 #        be exactly identical to the ones in the original process.
 # change helicity False: perform the reweighting by helicity summed 
 #        matrix-elements even if the events have been written with a 
-#        single helicity state.
-#        
-#*************************************************************************
-# ENTER YOUR COMMANDS BELOW. 
-#*************************************************************************
+#        single helicity state.       
+# change mode XXX: change the type of reweighting performed.
+#        allowed values: LO, NLO, LO+NLO
+#        - This command has no effect for reweighting an .lhe event file with LO accuracy. 
+#          In that case LO mode is always used (whatever entry is set).
+#        - When the .lhe file reweighted is at NLO accuracy, then all modes are allowed.
+#          * "LO" is an approximate leading order method
+#          * "NLO" is the NLO accurate method
+#          * "LO+NLO" runs both
+#       - "NLO" and "LO+NLO" modes requires 'keep_rwgt_info' equals True (run_card.dat)
+#          If the reweighting is done at generation level this parameter will
+#          automatically be set on True.
+#************************************************************************
 
-launch
-#  SPECIFY A PATH OR USE THE SET COMMAND LIKE THIS:
-# set sminputs 1 130 # modify 1/alpha_EW

=== modified file 'Template/LO/Source/make_opts'
--- Template/LO/Source/make_opts	2016-03-04 16:16:23 +0000
+++ Template/LO/Source/make_opts	2016-03-16 17:59:03 +0000
@@ -1,18 +1,35 @@
-# Environment variables
+DEFAULT_F2PY_COMPILER=f2py
+DEFAULT_F_COMPILER=gfortran
+MACFLAG=-mmacosx-version-min=10.7
+DEFAULT_CPP_COMPILER=clang
+STDLIB=-lc++
+STDLIB_FLAG=-stdlib=libc++
+#end_of_make_opts_variables
+# Rest of the makefile
 
 ifeq ($(origin FFLAGS),undefined)
-  FFLAGS= -O -w -fbounds-check -fPIC
-  #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
-endif
+FFLAGS= -O -w -fbounds-check -fPIC
+#FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
+endif
+
+
+# REMOVE MACFLAG IF NOT ON MAC OR FOR F2PY
+UNAME := $(shell uname -s)
+ifdef f2pymode
+MACFLAG=
+else
+ifneq ($(UNAME), Darwin)
+MACFLAG=
+endif
+endif
+
 
 ifeq ($(origin CXXFLAGS),undefined)
-  #CXXFLAGS= -g 
-  CXXFLAGS= -O -stdlib=libc++ -mmacosx-version-min=10.7
+CXXFLAGS= -O -stdlib=libc++ $(MACFLAG)
 endif
 
 ifeq ($(origin CFLAGS),undefined)
-  #CFLAGS= -g 
-  CFLAGS= -O -stdlib=libc++ -mmacosx-version-min=10.7
+CFLAGS= -O -stdlib=libc++ $(MACFLAG)
 endif
 
 # Increase the number of allowed charcters in a Fortran line
@@ -20,57 +37,57 @@
 
 # Set FC unless it's defined by an environment variable
 ifeq ($(origin FC),default)
-   FC=gfortran
+FC=gfortran
 endif
 ifeq ($(origin F2PY), undefined)
-   F2PY=f2py
+F2PY=f2py
 endif
 
 UNAME := $(shell uname -s)
 ifeq ($(origin LDFLAGS), undefined)
-LDFLAGS=-lc++  -mmacosx-version-min=10.7
+LDFLAGS=$(STDLIB) $(MACFLAG)
 endif
 
 # Options: dynamic, lhapdf
 # Option dynamic
 
 ifeq ($(UNAME), Darwin)
-  dylibext=dylib
+dylibext=dylib
 else
-  dylibext=so  
+dylibext=so
 endif
 
 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
+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
-  CXXFLAGS += $(shell $(lhapdf) --cppflags)
-  alfas_functions=alfas_functions_lhapdf
-  llhapdf=-lLHAPDF
+CXXFLAGS += $(shell $(lhapdf) --cppflags)
+alfas_functions=alfas_functions_lhapdf
+llhapdf+= -lLHAPDF
 else
-  alfas_functions=alfas_functions
-  llhapdf=
-endif
+alfas_functions=alfas_functions
+llhapdf=
+endif
\ No newline at end of file

=== modified file 'Template/LO/bin/madevent'
--- Template/LO/bin/madevent	2015-03-14 02:21:38 +0000
+++ Template/LO/bin/madevent	2016-03-16 17:59:03 +0000
@@ -164,7 +164,10 @@
             cmd_line.run_cmd('quit')
     elif args:
         # a single command is provided
-        cmd_line = cmd_interface.MadEventCmd()
+        if options.web:
+            cmd_line = cmd_interface.MadEventCmd()
+        else:
+            cmd_line = cmd_interface.MadEventCmdShell()
         if not hasattr(cmd_line, 'do_%s' % args[0]):
             if parser_error:
                 print parser_error

=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2015-11-02 16:43:45 +0000
+++ Template/NLO/Cards/run_card.dat	2016-03-16 17:59:03 +0000
@@ -81,6 +81,7 @@
  %(muf1_ref_fixed)s   = muF1_ref_fixed   ! fixed fact reference scale for pdf1
  %(muf2_ref_fixed)s   = muF2_ref_fixed   ! fixed fact reference scale for pdf2
  %(dynamical_scale_choice)s = dynamical_scale_choice ! Choose one of the preselected dynamical choices
+ %(keep_rwgt_info)s = keep_rwgt_info ! keep the information for reweighting (leads to an .lhe event file about ~5 times larger)
 #***********************************************************************
 # Renormalization and factorization scales (advanced and NLO options)  *
 #***********************************************************************

=== modified file 'Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_V.f'
--- Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_V.f	2015-12-04 15:47:32 +0000
+++ Template/NLO/MCatNLO/HWAnalyzer/hw6an_HwU_pp_V.f	2016-03-16 17:59:03 +0000
@@ -129,7 +129,7 @@
       ENDIF
 c
 c CHOOSE IDENT=24 FOR W+, IDENT=-24 FOR W-, IDENT=23 FOR Z0
-      IDENT=24
+      IDENT=23
 C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT''S A POWER-SUPPRESSED
 C EFFECT, SO THROW THE EVENT AWAY
       IF(SIGN(1.D0,PHEP(3,4)).EQ.SIGN(1.D0,PHEP(3,5)))THEN

=== modified file 'Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f'
--- Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f	2015-06-12 07:34:40 +0000
+++ Template/NLO/MCatNLO/srcHerwig/madfks_hwlhin.f	2016-03-16 17:59:03 +0000
@@ -127,13 +127,16 @@
             read(iunit,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
           enddo
           read(iunit,'(a)')string
-        elseif(jwgtinfo.eq.9)then
+        elseif(abs(jwgtinfo).eq.9)then
           if (numscales.eq.0 .and. numPDFpairs.eq.0) then
              write (*,*) 'event file not correct format'
              write(*,*)'FATAL ERROR #3 IN UPEVNT'
              stop
           endif
-          read(iunit,'(a)')string
+          string = ''
+          do while(index(string,'</event>').eq.0.and.index(string,'<rwgt>').eq.0)
+            read(iunit,'(a)')string
+          enddo
           wgtref=XWGTUP/MQQ
           do i=1,numscales
              do j=1,numscales
@@ -293,7 +296,7 @@
             ENDDO
             nwgt=nwgt+numPDFpairs
             numPDFpairs=numPDFpairs/2
-         ELSEIF( INDEX(STRING,"<weightgroup type='mg_reweighting'").ne.0
+         ELSEIF( INDEX(STRING,"<weightgroup name='mg_reweighting'").ne.0
      $           .and.STRING(1:1).ne.'#') then
             DO WHILE (.TRUE.)
                READ(61,'(a)') STRING

=== modified file 'Template/NLO/MCatNLO/srcPythia/madfks_pylhin.f'
--- Template/NLO/MCatNLO/srcPythia/madfks_pylhin.f	2015-06-12 07:34:40 +0000
+++ Template/NLO/MCatNLO/srcPythia/madfks_pylhin.f	2016-03-16 17:59:03 +0000
@@ -130,13 +130,16 @@
             read(iunit,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
           enddo
           read(iunit,'(a)')string
-        elseif(jwgtinfo.eq.9)then
+        elseif(abs(jwgtinfo).eq.9)then
           if (numscales.eq.0 .and. numPDFpairs.eq.0) then
              write (*,*) 'event file not correct format'
              write(*,*)'FATAL ERROR #3 IN UPEVNT'
              stop
           endif
-          read(iunit,'(a)')string
+          string = ''
+          do while(index(string,'</event>').eq.0.and.index(string,'<rwgt>').eq.0)
+            read(iunit,'(a)')string
+          enddo
           wgtref=XWGTUP/MQQ
           do i=1,numscales
              do j=1,numscales
@@ -295,7 +298,7 @@
             ENDDO
             nwgt=nwgt+numPDFpairs
             numPDFpairs=numPDFpairs/2
-         ELSEIF( INDEX(STRING,"<weightgroup type='mg_reweighting'").ne.0
+         ELSEIF( INDEX(STRING,"<weightgroup name='mg_reweighting'").ne.0
      $           .and.STRING(1:1).ne.'#') then
             DO WHILE (.TRUE.)
                READ(61,'(a)') STRING

=== modified file 'Template/NLO/Source/make_opts.inc'
--- Template/NLO/Source/make_opts.inc	2016-02-12 00:05:24 +0000
+++ Template/NLO/Source/make_opts.inc	2016-03-16 17:59:03 +0000
@@ -1,19 +1,38 @@
-# Environment variables
+STDLIB = -lc++
+STDLIB_FLAG = -stdlib=libc++
+MACFLAG = -mmacosx-version-min=10.7
+DEFAULT_F_COMPILER = gfortran
+DEFAULT_F2PY_COMPILER = f2py
+DEFAULT_CPP_COMPILER = g++
+#end_of_make_opts_variables
 
+# Rest of the makefile
 ifeq ($(origin FFLAGS),undefined)
   #FFLAGS= -O -w
   FFLAGS = -O -fno-automatic
   #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
 endif
 
+# REMOVE MACFLAG IF NOT ON MAC 
+UNAME := $(shell uname -s)
+ifneq ($(UNAME), Darwin)
+    MACFLAG=  
+endif
+
+
 ifeq ($(origin CXXFLAGS),undefined)
   #CXXFLAGS= -g 
-  CXXFLAGS = -O 
+  CXXFLAGS = -O $(STDLIB_FLAG) $(MACFLAG)
 endif
 
 ifeq ($(origin CFLAGS),undefined)
   #CFLAGS= -g 
-  CFLAGS = -O 
+  CFLAGS = -O  $(STDLIB_FLAG) $(MACFLAG)
+endif
+
+# REMOVE MACFLAG IF IN F2PY mode
+ifdef f2pymode
+    MACFLAG= 
 endif
 
 # Increase the number of allowed charcters in a Fortran line
@@ -22,21 +41,21 @@
 # Set FC unless it's defined by an environment variable
 ifeq ($(origin FC),default)
 #  FC=g77
-  FC=gfortran
+  FC=$(DEFAULT_F_COMPILER)
 endif
 # Set F2PY unless it's defined by an environment variable
 ifeq ($(origin F2PY),undefined)
-  F2PY=f2py
+  F2PY=$(DEFAULT_F2PY_COMPILER)
 endif
 
 # Set CXX unless it's defined by an environment variable
 ifeq ($(origin CXX),default)
-  CXX=g++
+  CXX=$(DEFAULT_CPP_COMPILER)
 endif
 
 UNAME := $(shell uname -s)
 ifeq ($(origin LDFLAGS), undefined)
-     LDFLAGS=-lc++
+     LDFLAGS=$(STDLIB)
 endif
 # Options: dynamic, lhapdf
 
@@ -73,7 +92,7 @@
 ifdef lhapdf
   CXXFLAGS += $(shell $(lhapdf) --cppflags)
   alfas_functions=alfas_functions_lhapdf
-  llhapdf=-lLHAPDF 
+  llhapdf+=-lLHAPDF 
   reweight_xsec_events_pdf_dummy=
 else
   alfas_functions=alfas_functions

=== modified file 'Template/NLO/Source/run.inc'
--- Template/NLO/Source/run.inc	2015-09-14 18:56:28 +0000
+++ Template/NLO/Source/run.inc	2016-03-16 17:59:03 +0000
@@ -73,9 +73,10 @@
 c
 c     Reweight info
 c
-      logical do_rwgt_scale, do_rwgt_pdf
+      logical do_rwgt_scale, do_rwgt_pdf, keep_rwgt_info
       integer pdf_set_min, pdf_set_max
       double precision rw_Fscale_down, rw_Fscale_up, rw_Rscale_down, rw_Rscale_up
       common/to_rwgt/ do_rwgt_scale, rw_Fscale_down, rw_Fscale_up, rw_Rscale_down, rw_Rscale_up,
-     # do_rwgt_pdf, pdf_set_min, pdf_set_max 
+     # do_rwgt_pdf, pdf_set_min, pdf_set_max,
+     # keep_rwgt_info 
 

=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
--- Template/NLO/SubProcesses/BinothLHA.f	2015-08-18 08:10:27 +0000
+++ Template/NLO/SubProcesses/BinothLHA.f	2016-03-16 17:59:03 +0000
@@ -3,20 +3,7 @@
 c Given the Born momenta, this is the Binoth-Les Houches interface file
 c that calls the OLP and returns the virtual weights. For convenience
 c also the born_wgt is passed to this subroutine.
-c
-C************************************************************************
-c WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING
-C************************************************************************
-c The Born in MadFKS -- and therefore also the virtual!-- should have a
-c slightly adapted identical particle symmetry factor. The normal
-c virtual weight as coming from the OLP should be divided by the number
-c of gluons in the corresponding real-emission process (i.e.  the number
-c of gluons in the Born plus one). This factor is passed to this
-c subroutine in /numberofparticles/ common block, as "ngluons". So,
-c divided virt_wgt by dble(ngluons) to get the correct virtual to be
-c used in MadFKS. The born_wgt that is passed to this subroutine has
-c already been divided by this factor.
-C************************************************************************
+cC************************************************************************
 c
       implicit none
       include "nexternal.inc"
@@ -34,11 +21,6 @@
       double precision, allocatable :: virt_wgts_hel(:,:)
       double precision mu,ao2pi,conversion,alpha_S
       save conversion
-      double precision fkssymmetryfactor,fkssymmetryfactorBorn,
-     &     fkssymmetryfactorDeg
-      integer ngluons,nquarks(-6:6)
-      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
-     &                         fkssymmetryfactorDeg,ngluons,nquarks
       logical firsttime,firsttime_conversion
       data firsttime,firsttime_conversion /.true.,.true./
       logical firsttime_run
@@ -112,9 +94,9 @@
          call sloopmatrix_thres(p, virt_wgts, tolerance, accuracies,
      $        ret_code)
          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)
+         virt_wgt= virt_wgts(1,0)
+         single  = virt_wgts(2,0)
+         double  = virt_wgts(3,0)
       else
          tolerance=PrecisionVirtualAtRunTime
 c Just set the accuracy found to a positive value as it is not specified
@@ -123,9 +105,9 @@
             call sloopmatrix_thres(p,virt_wgts,tolerance,accuracies
      $           ,ret_code)
             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)
+            virt_wgt= virt_wgts(1,0)
+            single  = virt_wgts(2,0)
+            double  = virt_wgts(3,0)
          elseif (mc_hel.eq.1) then
 c Use the Born helicity amplitudes to sample the helicities of the
 c virtual as flat as possible
@@ -148,14 +130,13 @@
             call sloopmatrixhel_thres(p,hel(ihel),virt_wgts_hel
      $           ,tolerance,accuracies,ret_code)
             prec_found = accuracies(0)
+c Average over initial state helicities
             virt_wgt = virt_wgt + virt_wgts_hel(1,0)*dble(goodhel(ihel))
-     $           /volh/4d0/dble(ngluons)
+     $           /volh/4d0
             single   = single   + virt_wgts_hel(2,0)*dble(goodhel(ihel))
-     $           /volh/4d0/dble(ngluons)
+     $           /volh/4d0
             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)
+     $           /volh/4d0
             if (nincoming.ne.2) then
                write (*,*)
      &              'Cannot do MC over helicities for 1->N processes'

=== modified file 'Template/NLO/SubProcesses/BinothLHA_OLP.f'
--- Template/NLO/SubProcesses/BinothLHA_OLP.f	2015-02-19 07:39:52 +0000
+++ Template/NLO/SubProcesses/BinothLHA_OLP.f	2016-03-16 17:59:03 +0000
@@ -5,18 +5,6 @@
 c also the born_wgt is passed to this subroutine.
 c
 C************************************************************************
-c WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING
-C************************************************************************
-c The Born in MadFKS -- and therefore also the virtual!-- should have a
-c slightly adapted identical particle symmetry factor. The normal
-c virtual weight as coming from the OLP should be divided by the number
-c of gluons in the corresponding real-emission process (i.e.  the number
-c of gluons in the Born plus one). This factor is passed to this
-c subroutine in /numberofparticles/ common block, as "ngluons". So,
-c divided virt_wgt by dble(ngluons) to get the correct virtual to be
-c used in MadFKS. The born_wgt that is passed to this subroutine has
-c already been divided by this factor.
-C************************************************************************
 c
       implicit none
       include "nexternal.inc"
@@ -29,11 +17,6 @@
       double precision virt_wgt,born_wgt,double,single,born,virt_wgts(4)
       double precision mu,ao2pi,conversion,alpha_S
       save conversion
-      double precision fkssymmetryfactor,fkssymmetryfactorBorn,
-     &     fkssymmetryfactorDeg
-      integer ngluons,nquarks(-6:6)
-      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
-     &                         fkssymmetryfactorDeg,ngluons,nquarks
       logical firsttime_pole,firsttime_conversion,firsttime_init
       data firsttime_pole,firsttime_conversion,firsttime_init
      &     /.true.,.true.,.true./
@@ -75,18 +58,17 @@
 c======================================================================
 c Replace the following line with the call to the one-loop code you wish
 c to use. virt_wgts contains finite part, single and double pole and the
-c Born. To understand why we need to divide them by ngluons, read the
-c warning message at the beginning of this subroutine
+c Born. 
 c
       if (firsttime_init) then
          call BinothLHAInit()
          firsttime_init=.false.
       endif
       call OLP_EvalSubProcess(proc_label,p,mu_r,alpha_S,virt_wgts)
-      double  = virt_wgts(1)/dble(ngluons)
-      single  = virt_wgts(2)/dble(ngluons)
-      virt_wgt= virt_wgts(3)/dble(ngluons)
-      born    = virt_wgts(4)/dble(ngluons)
+      double  = virt_wgts(1)
+      single  = virt_wgts(2)
+      virt_wgt= virt_wgts(3)
+      born    = virt_wgts(4)
 c======================================================================
 
 c======================================================================

=== modified file 'Template/NLO/SubProcesses/BinothLHA_user.f'
--- Template/NLO/SubProcesses/BinothLHA_user.f	2015-02-19 07:39:52 +0000
+++ Template/NLO/SubProcesses/BinothLHA_user.f	2016-03-16 17:59:03 +0000
@@ -5,18 +5,6 @@
 c also the born_wgt is passed to this subroutine.
 c
 C************************************************************************
-c WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING
-C************************************************************************
-c The Born in MadFKS -- and therefore also the virtual!-- should have a
-c slightly adapted identical particle symmetry factor. The normal
-c virtual weight as coming from the OLP should be divided by the number
-c of gluons in the corresponding real-emission process (i.e.  the number
-c of gluons in the Born plus one). This factor is passed to this
-c subroutine in /numberofparticles/ common block, as "ngluons". So,
-c divided virt_wgt by dble(ngluons) to get the correct virtual to be
-c used in MadFKS. The born_wgt that is passed to this subroutine has
-c already been divided by this factor.
-C************************************************************************
 c
       implicit none
       include "nexternal.inc"
@@ -27,11 +15,6 @@
       double precision virt_wgt,born_wgt,double,single,virt_wgts(3)
       double precision mu,ao2pi,conversion,alpha_S
       save conversion
-      double precision fkssymmetryfactor,fkssymmetryfactorBorn,
-     &     fkssymmetryfactorDeg
-      integer ngluons,nquarks(-6:6)
-      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
-     &                         fkssymmetryfactorDeg,ngluons,nquarks
       logical firsttime,firsttime_conversion
       data firsttime,firsttime_conversion /.true.,.true./
       integer           isum_hel
@@ -62,13 +45,11 @@
 c======================================================================
 c Replace the following line with the call to the one-loop code you 
 c wish to use. virt_wgts contains finite part, single and double pole
-c To understand why we need to divide them by ngluons, read the warning
-c message at the beginning of this subroutine
 c      
 c      call sloopmatrix(p, virt_wgts)
-c      virt_wgt= virt_wgts(1)/dble(ngluons)
-c      single  = virt_wgts(2)/dble(ngluons)
-c      double  = virt_wgts(3)/dble(ngluons)
+c      virt_wgt= virt_wgts(1)
+c      single  = virt_wgts(2)
+c      double  = virt_wgts(3)
 c
 c======================================================================
 

=== modified file 'Template/NLO/SubProcesses/c_weight.inc'
--- Template/NLO/SubProcesses/c_weight.inc	2015-02-13 20:03:34 +0000
+++ Template/NLO/SubProcesses/c_weight.inc	2016-03-16 17:59:03 +0000
@@ -4,20 +4,21 @@
       parameter (max_contr=128,max_wgt=256,max_iproc=100)
       logical H_event(max_contr)
       integer nFKS(max_contr),itype(max_contr),QCDpower(max_contr)
-     $     ,pdg(nexternal,max_contr),pdg_uborn(nexternal,max_contr)
+     $     ,pdg(nexternal,0:max_contr),pdg_uborn(nexternal,0:max_contr)
      $     ,parton_pdg_uborn(nexternal,max_iproc,max_contr)
      $     ,parton_pdg(nexternal,max_iproc,max_contr),icontr,iwgt
      $     ,plot_id(max_contr),niproc(max_contr),parton_pdf(nexternal
      $     ,max_iproc,max_contr),icontr_sum(0:max_iproc,max_contr)
      $     ,icontr_picked,iproc_picked
       double precision momenta(0:3,nexternal,max_contr),momenta_m(0:3
-     &     ,nexternal,max_contr),wgt(3,max_contr),bjx(2,max_contr)
-     &     ,scales2(3,max_contr),g_strong(max_contr),wgts(max_wgt
-     &     ,max_contr),parton_iproc(max_iproc,max_contr)
-     &     ,y_bst(max_contr),plot_wgts(max_wgt,max_contr)
-     &     ,shower_scale(max_contr),unwgt(max_iproc,max_contr)
-      common /c_weights/momenta,momenta_m,wgt,wgts,plot_wgts
-     &     ,parton_iproc,bjx,scales2,g_strong,y_bst,shower_scale,unwgt
-     &     ,parton_pdg,nFKS,itype,QCDpower,pdg,pdg_uborn
-     &     ,parton_pdg_uborn,plot_id,niproc,icontr,iwgt,icontr_sum
-     &     ,icontr_picked,iproc_picked,H_event
+     $     ,nexternal,2,max_contr),wgt(3,max_contr),wgt_ME_tree(2
+     $     ,max_contr),bjx(2,max_contr),scales2(3,max_contr)
+     $     ,g_strong(max_contr),wgts(max_wgt,max_contr)
+     $     ,parton_iproc(max_iproc,max_contr),y_bst(max_contr)
+     $     ,plot_wgts(max_wgt,max_contr),shower_scale(max_contr)
+     $     ,unwgt(max_iproc,max_contr)
+      common /c_weights/momenta,momenta_m,wgt,wgt_ME_tree,wgts,plot_wgts
+     $     ,parton_iproc,bjx,scales2,g_strong,y_bst,shower_scale,unwgt
+     $     ,parton_pdg,nFKS,itype,QCDpower,pdg,pdg_uborn
+     $     ,parton_pdg_uborn,plot_id,niproc,icontr,iwgt,icontr_sum
+     $     ,icontr_picked,iproc_picked,H_event

=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2015-10-20 13:26:25 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2016-03-16 17:59:03 +0000
@@ -182,9 +182,9 @@
      1 accuracies,return_code)
           accuracy=accuracies(0)
 
-          finite = virt_wgts(1,0)/dble(ngluons)
-          single = virt_wgts(2,0)/dble(ngluons)
-          double = virt_wgts(3,0)/dble(ngluons)
+          finite = virt_wgts(1,0)
+          single = virt_wgts(2,0)
+          double = virt_wgts(3,0)
 
 C         If MadLoop was still in initialization mode, then skip this
 C         point for the checks

=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2015-10-20 13:26:25 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2016-03-16 17:59:03 +0000
@@ -395,6 +395,8 @@
       common /to_abrv/ abrv
       integer iappl
       common /for_applgrid/ iappl
+      double precision       wgt_ME_born,wgt_ME_real
+      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
       if (ifl.ne.0) then
          write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
          stop 1
@@ -410,6 +412,8 @@
       virt_wgt_mint=0d0
       born_wgt_mint=0d0
       virtual_over_born=0d0
+      wgt_me_born=0d0
+      wgt_me_real=0d0
       if (ickkw.eq.-1) H1_factor_virt=0d0
       if (ickkw.eq.3) call set_FxFx_scale(0,p)
       call update_vegas_x(xx,x)
@@ -452,6 +456,8 @@
          MC_int_wgt=1d0/vol
       endif
       do iFKS=nFKS_min,nFKS_max
+         wgt_me_born=0d0
+         wgt_me_real=0d0
          jac=MC_int_wgt
          call update_fks_dir(iFKS,iconfig)
          call generate_momenta(ndim,iconfig,jac,x,p)

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2016-03-16 17:59:03 +0000
@@ -792,6 +792,8 @@
       logical               only_virt
       integer         imode
       common /c_imode/imode,only_virt
+      double precision       wgt_ME_born,wgt_ME_real
+      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
       sigintF=0d0
 c Find the nFKSprocess for which we compute the Born-like contributions
       if (firsttime) then
@@ -818,6 +820,8 @@
          born_wgt_mint=0d0
          virtual_over_born=0d0
          MCcntcalled=.false.
+         wgt_me_real=0d0
+         wgt_me_born=0d0
          if (ickkw.eq.3) call set_FxFx_scale(0,p)
          call update_vegas_x(xx,x)
          call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
@@ -866,6 +870,8 @@
          if(sum.eq.0) calculatedBorn=.false.
          nbody=.false.
          do i=1,proc_map(proc_map(0,1),0)
+            wgt_me_real=0d0
+            wgt_me_born=0d0
             iFKS=proc_map(proc_map(0,1),i)
             call update_fks_dir(iFKS,iconfig)
             jac=1d0/vol1

=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2016-03-16 17:59:03 +0000
@@ -38,7 +38,7 @@
       include 'run.inc'
       include 'timing_variables.inc'
       double precision wgt1,wgt2,wgt3,bsv_wgt,virt_wgt,born_wgt,pi,g2
-     &     ,g22
+     &     ,g22,wgt4
       parameter (pi=3.1415926535897932385d0)
       double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
      $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
@@ -61,6 +61,7 @@
       g2=g**(nint(2*wgtbpower))
       g22=g**(nint(2*wgtbpower+2))
       wgt1=wgtnstmp*f_nb/g22
+      wgt4=wgtnstmp_avgvirt*f_nb/g22
       if (ickkw.eq.3 .and. fxfx_exp_rewgt.ne.0d0) then
          wgt1=wgt1 - fxfx_exp_rewgt*born_wgt*f_nb/g2/(4d0*pi)
       elseif (ickkw.eq.-1) then
@@ -78,6 +79,7 @@
       wgt2=wgtwnstmpmur*f_nb/g22
       wgt3=wgtwnstmpmuf*f_nb/g22
       call add_wgt(3,wgt1,wgt2,wgt3)
+      call add_wgt(15,wgt4,0d0,0d0)
 c Special for the soft-virtual needed for the virt-tricks. The
 c *_wgt_mint variable should be directly passed to the mint-integrator
 c and not be part of the plots nor computation of the cross section.
@@ -955,6 +957,7 @@
 c     type=12: MC subtraction with n-body kin.
 c     type=13: MC subtraction with n+1-body kin.
 c     type=14: virtual corrections
+c     type=15: virt-trick: average born contribution
 c     wgt1 : weight of the contribution not multiplying a scale log
 c     wgt2 : coefficient of the weight multiplying the log[mu_R^2/Q^2]
 c     wgt3 : coefficient of the weight multiplying the log[mu_F^2/Q^2]
@@ -990,6 +993,9 @@
 c        PDG=21 (gluon) code.
 c     If the contribution belongs to an H-event or S-event:
 c        H_event(icontr)
+c     The weight of the born or real-emission matrix element
+c        corresponding to this contribution: wgt_ME_tree. This weight does
+c        include the 'ngluon' correction factor for the Born.
 c
 c Not set in this subroutine, but included in the c_weights common block
 c are the
@@ -1033,6 +1039,15 @@
       double precision    p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
      $                    ,pswgt_cnt(-2:2),jac_cnt(-2:2)
       common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
+      double precision         fkssymmetryfactor,fkssymmetryfactorBorn,
+     &                         fkssymmetryfactorDeg
+      integer                                      ngluons,nquarks(-6:6)
+      common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
+     &                         fkssymmetryfactorDeg,ngluons,nquarks
+      double precision       wgt_ME_born,wgt_ME_real
+      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
+      double precision     iden_comp
+      common /c_iden_comp/ iden_comp
       if (wgt1.eq.0d0 .and. wgt2.eq.0d0 .and. wgt3.eq.0d0) return
 c Check for NaN's and INF's. Simply skip the contribution
       if (wgt1.ne.wgt1) return
@@ -1064,57 +1079,49 @@
 c     Anything else
          QCDpower(icontr)=nint(2*wgtbpower+2)
       endif
+c Compensate for the fact that in the Born matrix elements, we use the
+c identical particle symmetry factor of the corresponding real emission
+c matrix elements
+c IDEN_COMP STUFF NEEDS TO BE UPDATED WHEN MERGING WITH 'FKS_EW' STUFF
+      wgt_ME_tree(1,icontr)=wgt_me_born
+      wgt_ME_tree(2,icontr)=wgt_me_real
+      do i=1,nexternal
+         do j=0,3
+            if (p1_cnt(0,1,0).gt.0d0) then
+               momenta_m(j,i,1,icontr)=p1_cnt(j,i,0)
+            elseif (p1_cnt(0,1,1).gt.0d0) then
+               momenta_m(j,i,1,icontr)=p1_cnt(j,i,1)
+            elseif (p1_cnt(0,1,2).gt.0d0) then
+               momenta_m(j,i,1,icontr)=p1_cnt(j,i,2)
+            else
+               if (i.lt.fks_i_d(nFKSprocess)) then
+                  momenta_m(j,i,1,icontr)=p_born(j,i)
+               elseif(i.eq.fks_i_d(nFKSprocess)) then
+                  momenta_m(j,i,1,icontr)=0d0
+               else
+                  momenta_m(j,i,1,icontr)=p_born(j,i-1)
+               endif
+            endif
+            momenta_m(j,i,2,icontr)=p_ev(j,i)
+         enddo
+      enddo
       if(type.eq.1 .or. type.eq. 8 .or. type.eq.9 .or. type.eq.10 .or.
      &     type.eq.13) then
 c real emission and n+1-body kin. contributions to counter terms and MC
 c subtr term
          do i=1,nexternal
             do j=0,3
-c     'momenta' is the momentum configuration for this contribution;
-c     'momenta_m' is the momentum configuration that was used in the
-c     matrix elements of this contribution.
-               momenta(j,i,icontr)=p_ev(j,i)
-               if (type.eq.1) then
-                  momenta_m(j,i,icontr)=momenta(j,i,icontr)
-               else
-                  if (p1_cnt(0,1,0).gt.0d0) then
-                     momenta_m(j,i,icontr)=p1_cnt(j,i,0)
-                  elseif (p1_cnt(0,1,1).gt.0d0) then
-                     momenta_m(j,i,icontr)=p1_cnt(j,i,1)
-                  elseif (p1_cnt(0,1,2).gt.0d0) then
-                     momenta_m(j,i,icontr)=p1_cnt(j,i,2)
-                  else
-                     write (*,*) 'ERROR in add_wgt: no valid momenta'
-                     stop 1
-                  endif
-               endif
+               momenta(j,i,icontr)=momenta_m(j,i,2,icontr)
             enddo
          enddo
          H_event(icontr)=.true.
       elseif(type.ge.2 .and. type.le.7 .or. type.eq.11 .or. type.eq.12
-     $        .or. type.eq.14)then
+     $        .or. type.eq.14 .or. type.eq.15)then
 c Born, counter term, soft-virtual, or n-body kin. contributions to real
 c and MC subtraction terms.
          do i=1,nexternal
             do j=0,3
-c     'momenta' is the momentum configuration for this contribution;
-c     'momenta_m' is the momentum configuration that was used in the
-c     matrix elements of this contribution.
-               if (p1_cnt(0,1,0).gt.0d0) then
-                  momenta(j,i,icontr)=p1_cnt(j,i,0)
-               elseif (p1_cnt(0,1,1).gt.0d0) then
-                  momenta(j,i,icontr)=p1_cnt(j,i,1)
-               elseif (p1_cnt(0,1,2).gt.0d0) then
-                  momenta(j,i,icontr)=p1_cnt(j,i,2)
-               else
-                  write (*,*) 'ERROR in add_wgt: no valid momenta'
-                  stop 1
-               endif
-               if (type.eq.11) then
-                  momenta_m(j,i,icontr)=p_ev(j,i)
-               else
-                  momenta_m(j,i,icontr)=momenta(j,i,icontr)
-               endif
+               momenta(j,i,icontr)=momenta_m(j,i,1,icontr)
             enddo
          enddo
          H_event(icontr)=.false.
@@ -1494,7 +1501,7 @@
       integer iproc_save(fks_configs),eto(maxproc,fks_configs),
      &     etoi(maxproc,fks_configs),maxproc_found
       common/cproc_combination/iproc_save,eto,etoi,maxproc_found
-      if (icontr.gt.7) then
+      if (icontr.gt.8) then
          write (*,*) 'ERROR: too many applgrid weights. '/
      &        /'Should have at most one of each itype.',icontr
          stop 1
@@ -1536,8 +1543,8 @@
             appl_QES2(2)=scales2(1,i)
             appl_muR2(2)=scales2(2,i)
             appl_muF2(2)=scales2(3,i)
-         elseif (itype(i).eq.3 .or. itype(i).eq.4 .or. itype(i).eq.14)
-     $           then
+         elseif (itype(i).eq.3 .or. itype(i).eq.4 .or. itype(i).eq.14
+     &           .or. itype(i).eq.15)then
 c     virtual, soft-virtual or soft-counter
             appl_w0(2)=appl_w0(2)+wgt(1,i)*final_state_rescaling
             appl_wR(2)=appl_wR(2)+wgt(2,i)*final_state_rescaling
@@ -1635,7 +1642,7 @@
       if (icontr.eq.0) return
       do i=1,icontr
          if (itype(i).eq.2 .or. itype(i).eq.3 .or. itype(i).eq.14 .or.
-     &        itype(i).eq.7) then
+     &        itype(i).eq.7 .or. itype(i).eq.15) then
             sig=sig+wgts(1,i)
          endif
       enddo
@@ -1654,7 +1661,7 @@
       if (icontr.eq.0) return
       do i=1,icontr
          if (itype(i).ne.2 .and. itype(i).ne.3 .and. itype(i).ne.14
-     &        .and. itype(i).ne.7) then
+     &        .and. itype(i).ne.7 .and. itype(i).ne.15) then
             sig=sig+wgts(1,i)
          endif
       enddo
@@ -2054,8 +2061,9 @@
             tmp_wgt=0d0
             do j=1,icontr_sum(0,i)
                ict=icontr_sum(j,i)
-               if (itype(ict).ne.2 .and. itype(ict).ne.3 .and.
-     $             itype(ict).ne.14) tmp_wgt=tmp_wgt+wgts(1,ict)
+               if ( itype(ict).ne.2  .and. itype(ict).ne.3 .and.
+     $              itype(ict).ne.14 .and. itype(ict).ne.15)
+     $                              tmp_wgt=tmp_wgt+wgts(1,ict)
             enddo
             n1body_wgt=n1body_wgt+abs(tmp_wgt)
          enddo
@@ -2162,7 +2170,8 @@
       include 'reweight0.inc'
       include 'genps.inc'
       include 'nFKSconfigs.inc'
-      integer i,ii,j,jj,ict,ipr,momenta_conf
+      include 'fks_info.inc'
+      integer k,i,ii,j,jj,ict,ipr,momenta_conf(2)
       logical momenta_equal,found
       double precision conv,momenta_str_l(0:3,nexternal,max_n_ctr)
       external momenta_equal
@@ -2183,24 +2192,32 @@
 c Check if the current set of momenta are already available in the
 c momenta_str_l array. If not, add it.
          found=.false.
-         do j=1,n_mom_conf
-            if (momenta_equal(momenta_str_l(0,1,j),momenta_m(0,1,ict)))
-     &           then
-               momenta_conf=j
-               found=.true.
-               exit
+         do k=1,2
+            do j=1,n_mom_conf
+               if (momenta_m(0,1,k,ict).le.0d0) then
+                  momenta_conf(k)=0
+                  cycle
+               endif
+               if (momenta_equal(momenta_str_l(0,1,j),
+     &                           momenta_m(0,1,k,ict))) then
+                  momenta_conf(k)=j
+                  found=.true.
+                  exit
+               endif
+            enddo
+            if (.not. found) then
+               n_mom_conf=n_mom_conf+1
+               do ii=1,nexternal
+                  do jj=0,3
+                     momenta_str(jj,ii,n_mom_conf)=
+     &                                      momenta_m(jj,ii,k,ict)
+                     momenta_str_l(jj,ii,n_mom_conf)=
+     &                                      momenta_m(jj,ii,k,ict)
+                  enddo
+               enddo
+               momenta_conf(k)=n_mom_conf
             endif
          enddo
-         if (.not. found) then
-            n_mom_conf=n_mom_conf+1
-            do ii=1,nexternal
-               do jj=0,3
-                  momenta_str(jj,ii,n_mom_conf)=momenta_m(jj,ii,ict)
-                  momenta_str_l(jj,ii,n_mom_conf)=momenta_m(jj,ii,ict)
-               enddo
-            enddo
-            momenta_conf=n_mom_conf
-         endif
          if (.not. Hevents) then
 c For S-events, be careful to take all the IPROC that contribute to the
 c iproc_picked:
@@ -2208,15 +2225,17 @@
             do ii=1,iproc_save(nFKS(ict))
                if (eto(ii,nFKS(ict)).ne.ipr) cycle
                n_ctr_found=n_ctr_found+1
+
                if (nincoming.eq.2) then
-                  write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
-     &                 (wgt(j,ict)*conv,j=1,3),
+                  write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
+     &                 (wgt(j,ict)*conv,j=1,3),(wgt_me_tree(j,ict),j=1,2),
      &                 nexternal
                else
-                  write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
-     &                 (wgt(j,ict),j=1,3),
+                  write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
+     &                 (wgt(j,ict),j=1,3),(wgt_me_tree(j,ict),j=1,2), 
      &                 nexternal
                endif
+
                procid=''
                do j=1,nexternal
                   write (str_temp,*) parton_pdg(j,ii,ict)
@@ -2226,13 +2245,18 @@
                n_ctr_str(n_ctr_found) =
      &              trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &              //trim(adjustl(procid))
-               write (str_temp,'(i2,5(1x,d14.8),3(1x,i2),1x,d18.12)')
+
+               write (str_temp,'(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12)')
      &              QCDpower(ict),
      &              (bjx(j,ict),j=1,2),
      &              (scales2(j,ict),j=1,3),
-     &              momenta_conf,
+     &              g_strong(ict),
+     &              (momenta_conf(j),j=1,2),
      &              itype(ict),
      &              nFKS(ict),
+     &              fks_i_d(nFKS(ict)),
+     &              fks_j_d(nFKS(ict)),
+     &              parton_pdg_uborn(fks_j_d(nFKS(ict)),ii,ict),
      &              parton_iproc(ii,ict)
                n_ctr_str(n_ctr_found) =
      &              trim(adjustl(n_ctr_str(n_ctr_found)))//' '
@@ -2242,15 +2266,17 @@
 c H-event
             ipr=iproc_picked
             n_ctr_found=n_ctr_found+1
+
             if (nincoming.eq.2) then
-               write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
-     &              (wgt(j,ict)*conv,j=1,3),
+               write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
+     &              (wgt(j,ict)*conv,j=1,3),(wgt_me_tree(j,ict),j=1,2),
      &              nexternal
             else
-               write (n_ctr_str(n_ctr_found),'(3(1x,d18.12),1x,i2)')
-     &              (wgt(j,ict),j=1,3),
+               write (n_ctr_str(n_ctr_found),'(5(1x,d18.12),1x,i2)')
+     &              (wgt(j,ict),j=1,3),(wgt_me_tree(j,ict),j=1,2),
      &              nexternal
             endif
+
             procid=''
             do j=1,nexternal
                write (str_temp,*) parton_pdg(j,ipr,ict)
@@ -2260,17 +2286,24 @@
             n_ctr_str(n_ctr_found) =
      &           trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &           //trim(adjustl(procid))
-            write (str_temp,'(i2,5(1x,d14.8),3(1x,i2),1x,d18.12)')
+
+            write (str_temp,'(i2,6(1x,d14.8),6(1x,i2),1x,i8,1x,d18.12)')
      &           QCDpower(ict),
      &           (bjx(j,ict),j=1,2),
      &           (scales2(j,ict),j=1,3),
-     &           momenta_conf,
-     &              itype(ict),
-     &              nFKS(ict),
+     &           g_strong(ict),
+     &           (momenta_conf(j),j=1,2),
+     &           itype(ict),
+     &           nFKS(ict),
+     &           fks_i_d(nFKS(ict)),
+     &           fks_j_d(nFKS(ict)),
+     &           parton_pdg_uborn(fks_j_d(nFKS(ict)),ipr,ict),
      &           parton_iproc(ipr,ict)
             n_ctr_str(n_ctr_found) =
      &           trim(adjustl(n_ctr_str(n_ctr_found)))//' '
      &           //trim(adjustl(str_temp))
+
+
          endif
          if (n_ctr_found.ge.max_n_ctr) then
             write (*,*) 'ERROR: too many contributions in <rwgt>'
@@ -2852,6 +2885,8 @@
       parameter (vtiny=1d-8)
       double complex ximag
       parameter (ximag=(0.d0,1.d0))
+      double precision iden_comp
+      common /c_iden_comp/iden_comp
 C  
       if(p_born(0,1).le.0.d0)then
 c Unphysical kinematics: set matrix elements equal to zero
@@ -2923,7 +2958,7 @@
          write(*,*) 'FATAL ERROR in sborncol_fsr',i_type,j_type,i_fks,j_fks
          stop
       endif
-      wgt=dble(wgt1(1)*ap+wgt1(2)*Q)
+      wgt=dble(wgt1(1)*ap+wgt1(2)*Q)*iden_comp
       return
       end
 
@@ -2971,6 +3006,8 @@
       parameter (vtiny=1d-8)
       double complex ximag
       parameter (ximag=(0.d0,1.d0))
+      double precision iden_comp
+      common /c_iden_comp/iden_comp
 C  
       if(p_born(0,1).le.0.d0)then
 c Unphysical kinematics: set matrix elements equal to zero
@@ -3047,7 +3084,7 @@
      #             wgt1(2) * dconjg(azifact)
          call Qterms_reduced_spacelike(m_type, i_type, t, z, Q)
       endif
-      wgt=dble(wgt1(1)*ap+wgt1(2)*Q)
+      wgt=dble(wgt1(1)*ap+wgt1(2)*Q)*iden_comp
       return
       end
 
@@ -3329,6 +3366,8 @@
 
       double precision zero,pmass(nexternal)
       parameter(zero=0d0)
+      double precision iden_comp
+      common /c_iden_comp/iden_comp
       include "pmass.inc"
 c
 c Call the Born to be sure that 'CalculatedBorn' is done correctly. This
@@ -3353,8 +3392,7 @@
             endif
          enddo
       enddo
-
-      wgt=softcontr
+      wgt=softcontr*iden_comp
 c Add minus sign to compensate the minus in the color factor
 c of the color-linked Borns (b_sf_0??.f)
 c Factor two to fix the limits.
@@ -3487,6 +3525,8 @@
       double precision one,pi
       parameter (one=1.d0)
       parameter (pi=3.1415926535897932385d0)
+      double precision iden_comp
+      common /c_iden_comp/iden_comp
 
       if(j_fks.gt.nincoming)then
 c Do not include this contribution for final-state branchings
@@ -3544,7 +3584,7 @@
 c The partonic flux 1/(2*s) is inserted in genps. Thus, an extra 
 c factor z (implicit in the flux of the reduced Born in FKS) 
 c has to be inserted here
-      xnorm=1.d0/z
+      xnorm=1.d0/z *iden_comp
 
       collrem_xi=oo2pi * born_wgt * collrem_xi * xnorm
       collrem_lxi=oo2pi * born_wgt * collrem_lxi * xnorm
@@ -4343,7 +4383,7 @@
       include "run.inc"
       include "fks_powers.inc"
       include 'reweight.inc'
-      double precision p(0:3,nexternal),bsv_wgt,born_wgt
+      double precision p(0:3,nexternal),bsv_wgt,born_wgt,avv_wgt
       double precision pp(0:3,nexternal)
       
       double complex wgt1(2)
@@ -4441,6 +4481,7 @@
          bsv_wgt=dble(wgt1(1))
          born_wgt=dble(wgt1(1))
          virt_wgt=0d0
+         avv_wgt=0d0 
 
          if (abrv.eq.'born' .or. abrv.eq.'grid') goto 549
          if (abrv.eq.'virt' .or. abrv.eq.'viSC' .or.
@@ -4594,7 +4635,7 @@
 c$$$            bsv_wgt=bsv_wgt+virt_wgt_save
          endif
          if (abrv(1:4).ne.'virt' .and. ickkw.ne.-1)
-     &        bsv_wgt=bsv_wgt+average_virtual*born_wgt*ao2pi
+     &        avv_wgt=average_virtual*born_wgt*ao2pi
 
 c eq.(MadFKS.C.13)
          if(abrv.eq.'viSA'.or.abrv.eq.'viSB')then
@@ -4642,9 +4683,11 @@
             wgtnstmp=bsv_wgt-born_wgt-
      #                wgtwnstmpmuf*log(q2fact(1)/QES2)-
      #                wgtwnstmpmur*log(scale**2/QES2)
+            wgtnstmp_avgvirt = avv_wgt
          else
             wgtnstmp=0d0
             wgtwnstmpmur=0.d0
+            wgtnstmp_avgvirt = 0d0
          endif
 
          if (abrv(1:2).eq.'vi') then
@@ -5376,11 +5419,15 @@
       common/numberofparticles/fkssymmetryfactor,fkssymmetryfactorBorn,
      &                         fkssymmetryfactorDeg,ngluons,nquarks
 
+      double precision iden_comp
+      common /c_iden_comp/iden_comp
+
       include 'coupl.inc'
       include 'genps.inc'
       include 'nexternal.inc'
       include 'fks_powers.inc'
       include 'nFKSconfigs.inc'
+      include 'c_weight.inc'
       integer fks_j_from_i(nexternal,0:nexternal)
      &     ,particle_type(nexternal),pdg_type(nexternal)
       common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
@@ -5416,13 +5463,14 @@
       character*1 integrate
       integer i_fks,j_fks
       common/fks_indices/i_fks,j_fks
-      integer fac_i,fac_j,i_fks_pdg,j_fks_pdg
+      integer fac_i,fac_j,i_fks_pdg,j_fks_pdg,iden(nexternal)
 
       integer fac_i_FKS(fks_configs),fac_j_FKS(fks_configs)
-     $     ,i_type_FKS(fks_configs),j_type_FKS(fks_configs)
-     $     ,m_type_FKS(fks_configs),ngluons_FKS(fks_configs)
+     &     ,i_type_FKS(fks_configs),j_type_FKS(fks_configs)
+     &     ,m_type_FKS(fks_configs),ngluons_FKS(fks_configs)
+     &     ,iden_real_FKS(fks_configs),iden_born_FKS(fks_configs)
       save fac_i_FKS,fac_j_FKS,i_type_FKS,j_type_FKS,m_type_FKS
-     $     ,ngluons_FKS
+     $     ,ngluons_FKS,iden_real_FKS,iden_born_FKS
 
       character*13 filename
 
@@ -5598,12 +5646,56 @@
          i_type_FKS(nFKSprocess)=i_type
          j_type_FKS(nFKSprocess)=j_type
          m_type_FKS(nFKSprocess)=m_type
+
+
+c Compute the identical particle symmetry factor that is in the
+c real-emission matrix elements.
+         iden_real_FKS(nFKSprocess)=1
+         do i=1,nexternal
+            iden(i)=1
+         enddo
+         do i=nincoming+2,nexternal
+            do j=nincoming+1,i-1
+               if (pdg_type(j).eq.pdg_type(i)) then
+                  iden(j)=iden(j)+1
+                  iden_real_FKS(nFKSprocess)=
+     &                 iden_real_FKS(nFKSprocess)*iden(j)
+                  exit
+               endif
+            enddo
+         enddo
+c Compute the identical particle symmetry factor that is in the
+c Born matrix elements.
+         iden_born_FKS(nFKSprocess)=1
+         call set_pdg(0,nFKSprocess)
+         do i=1,nexternal
+            iden(i)=1
+         enddo
+         do i=nincoming+2,nexternal-1
+            do j=nincoming+1,i-1
+               if (pdg_uborn(j,0).eq.pdg_uborn(i,0)) then
+                  iden(j)=iden(j)+1
+                  iden_born_FKS(nFKSprocess)=
+     &                 iden_born_FKS(nFKSprocess)*iden(j)
+                  exit
+               endif
+            enddo
+         enddo
       endif
 
       i_type=i_type_FKS(nFKSprocess)
       j_type=j_type_FKS(nFKSprocess)
       m_type=m_type_FKS(nFKSprocess)
 
+c Difference in identical particle factor in the Born and real emission
+c matrix elements. To define wgt_ME_tree for the Born, we need to
+c include this factor, because in the current Born the symmetry factor
+c for the real is used. THIS NEEDS TO BE CHANGED WHEN MERGING WITH THE
+c 'FKS_EW' STUFF
+      iden_comp=dble(iden_born_FKS(nFKSprocess))/
+     &          dble(iden_real_FKS(nFKSprocess))
+      
+      
 c Set matrices used by MC counterterms
       if (match_to_shower) call set_mc_matrices
 
@@ -5614,7 +5706,7 @@
       if (nbody.and.pdg_type(i_fks).eq.21) then
          fkssymmetryfactor=dble(ngluons)
          fkssymmetryfactorDeg=dble(ngluons)
-         fkssymmetryfactorBorn=dble(ngluons)
+         fkssymmetryfactorBorn=1d0
       elseif(pdg_type(i_fks).eq.-21) then
          fkssymmetryfactor=1d0
          fkssymmetryfactorDeg=1d0

=== modified file 'Template/NLO/SubProcesses/handling_lhe_events.f'
--- Template/NLO/SubProcesses/handling_lhe_events.f	2015-03-03 10:11:02 +0000
+++ Template/NLO/SubProcesses/handling_lhe_events.f	2016-03-16 17:59:03 +0000
@@ -661,20 +661,21 @@
               stop
            endif
            write(ifile,'(a)')'  </rwgt>'
-         elseif(jwgtinfo.eq.-5) then
-           write(ifile,'(a)')'  <rwgt>'
+        endif
+        if(jwgtinfo.eq.-5.or.jwgtinfo.eq.-9.or.jwgtinfo.eq.-8) then
+           write(ifile,'(a)')'  <mgrwgt>'
            write (ifile,'(1x,d16.10,3(1x,i4))') wgtref,n_ctr_found
-     &          ,n_mom_conf,nint(wgtcpower)
+     &          ,n_mom_conf, nint(wgtcpower)
            do i=1,n_mom_conf
               do j=1,mexternal
-                 write (ifile,'(4(1x,d16.10))')
+                 write (ifile,'(4(1x,d21.15))')
      &                (momenta_str(ii,j,i),ii=0,3)
               enddo
            enddo
            do i=1,n_ctr_found
               write (ifile,'(a)') trim(adjustl(n_ctr_str(i)))
            enddo
-           write(ifile,'(a)')'  </rwgt>'
+           write(ifile,'(a)')'  </mgrwgt>'
          elseif(jwgtinfo.eq.15) then
            write(ifile,'(a)')'  <unlops>'
            write(ifile,*)NUP_H
@@ -685,7 +686,8 @@
      $             VTIMUP_H(I),SPINUP_H(I)
            enddo
            write(ifile,'(a)')'  </unlops>'
-        elseif(jwgtinfo.eq.8)then
+        endif
+        if(jwgtinfo.eq.8.or.jwgtinfo.eq.-8)then
            write(ifile,'(a)') '  <rwgt>'
           write(ifile,406)wgtref,wgtxsecmu(1,1),numscales,numPDFpairs
           do i=1,numscales
@@ -697,8 +699,7 @@
             write(ifile,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
           enddo
           write(ifile,'(a)') '  </rwgt>'
-
-        elseif(jwgtinfo.eq.9)then
+        elseif(jwgtinfo.eq.9.or.jwgtinfo.eq.-9)then
            write(ifile,'(a)') '  <rwgt>'
            do i=1,numscales
               do j=1,numscales
@@ -904,7 +905,7 @@
               stop
            endif
            read(ifile,'(a)')string
-         elseif(jwgtinfo.eq.-5) then
+         elseif(jwgtinfo.eq.-5.or.jwgtinfo.eq.-8.or.jwgtinfo.eq.-9) then
            read(ifile,'(a)')string
            read(ifile,*) wgtref,n_ctr_found,n_mom_conf,wgtcpower
            do i=1,n_mom_conf
@@ -916,7 +917,8 @@
               read (ifile,'(a)') n_ctr_str(i)
            enddo
            read(ifile,'(a)')string
-         elseif(jwgtinfo.eq.15) then
+         endif
+         if(jwgtinfo.eq.15) then
            read(ifile,'(a)') string
            read(ifile,*)NUP_H
            do i=1,NUP_H
@@ -926,7 +928,7 @@
      $             VTIMUP_H(I),SPINUP_H(I)
            enddo
            read(ifile,'(a)') string
-        elseif(jwgtinfo.eq.8)then
+        elseif(jwgtinfo.eq.8.or.jwgtinfo.eq.-8)then
           read(ifile,'(a)')string
           read(ifile,406)wgtref,wgtxsecmu(1,1),numscales,numPDFpairs
           do i=1,numscales
@@ -938,7 +940,7 @@
             read(ifile,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
           enddo
           read(ifile,'(a)')string
-        elseif(jwgtinfo.eq.9)then
+        elseif(jwgtinfo.eq.9.or.jwgtinfo.eq.-9)then
            read(ifile,'(a)')string
            wgtref=XWGTUP
            do i=1,numscales
@@ -1156,7 +1158,7 @@
               stop
            endif
            read(ifile,'(a)')string
-         elseif(jwgtinfo.eq.-5) then
+         elseif(jwgtinfo.eq.-5.or.jwgtinfo.eq.-8.or.jwgtinfo.eq.-9) then
            read(ifile,'(a)')string
            read(ifile,*) wgtref,n_ctr_found,n_mom_conf,wgtcpower
            do i=1,n_mom_conf
@@ -1168,7 +1170,8 @@
               read (ifile,'(a)') n_ctr_str(i)
            enddo
            read(ifile,'(a)')string
-         elseif(jwgtinfo.eq.15) then
+         endif
+         if(jwgtinfo.eq.15) then
            read(ifile,'(a)') string
            read(ifile,*)NUP_H
            do i=1,NUP_H
@@ -1178,7 +1181,7 @@
      $             VTIMUP_H(I),SPINUP_H(I)
            enddo
            read(ifile,'(a)') string
-        elseif(jwgtinfo.eq.8)then
+        elseif(jwgtinfo.eq.8.or.jwgtinfo.eq.-8)then
           read(ifile,'(a)')string
           read(ifile,406)wgtref,wgtxsecmu(1,1),numscales,numPDFpairs
           do i=1,numscales
@@ -1190,7 +1193,7 @@
             read(ifile,404)wgtxsecPDF(nps),wgtxsecPDF(nng)
           enddo
           read(ifile,'(a)')string
-        elseif(jwgtinfo.eq.9)then
+        elseif(jwgtinfo.eq.9.or.jwgtinfo.eq.-9)then
            read(ifile,'(a)')string
            wgtref=XWGTUP
            do i=1,numscales

=== modified file 'Template/NLO/SubProcesses/montecarlocounter.f'
--- Template/NLO/SubProcesses/montecarlocounter.f	2015-09-14 18:56:28 +0000
+++ Template/NLO/SubProcesses/montecarlocounter.f	2016-03-16 17:59:03 +0000
@@ -1014,11 +1014,13 @@
 
       logical calculatedBorn
       common/ccalculatedBorn/calculatedBorn
+      double precision iden_comp
+      common /c_iden_comp/iden_comp
 
 c Particle types (=color) of i_fks, j_fks and fks_mother
       integer i_type,j_type,m_type
       common/cparticle_types/i_type,j_type,m_type
-
+      
 c
 c BORN
       call sborn(p_born,wgt1)
@@ -1149,7 +1151,7 @@
 c BARRED AMPLITUDES
       do i=1,max_bcol
          if (sumborn.ne.0d0) then
-            bornbars(i)=jamp2(i)/sumborn * born
+            bornbars(i)=jamp2(i)/sumborn * born *iden_comp
          elseif (born.eq.0d0 .or. jamp2(i).eq.0d0) then
             bornbars(i)=0d0
          else
@@ -1157,7 +1159,7 @@
             stop
          endif
          if (sumborn.ne.0d0) then
-            bornbarstilde(i)=jamp2(i)/sumborn * borntilde
+            bornbarstilde(i)=jamp2(i)/sumborn * borntilde *iden_comp
          elseif (borntilde.eq.0d0 .or. jamp2(i).eq.0d0) then
             bornbarstilde(i)=0d0
          else

=== modified file 'Template/NLO/SubProcesses/reweight.inc'
--- Template/NLO/SubProcesses/reweight.inc	2015-02-13 14:49:55 +0000
+++ Template/NLO/SubProcesses/reweight.inc	2016-03-16 17:59:03 +0000
@@ -30,8 +30,8 @@
       common/cwgtaux1/wgtdegrem_xi,wgtdegrem_lxi,wgtdegrem_muF
 
 c Temporary quantities relevant to n-body terms
-      double precision wgtnstmp,wgtwnstmpmuf,wgtwnstmpmur,wgtref_nbody,wgtvirttmp
-      common/cwgtaux2/wgtnstmp,wgtwnstmpmuf,wgtwnstmpmur,wgtref_nbody,wgtvirttmp
+      double precision wgtnstmp,wgtwnstmpmuf,wgtwnstmpmur,wgtref_nbody,wgtvirttmp,wgtnstmp_avgvirt
+      common/cwgtaux2/wgtnstmp,wgtwnstmpmuf,wgtwnstmpmur,wgtref_nbody,wgtvirttmp,wgtnstmp_avgvirt
 
 c Temporary weights for pure-NLO computations: event, counterevents+Born, 
 c and Born, with the same conventions as in fks_singular

=== modified file 'Template/NLO/SubProcesses/reweight_xsec_events.f'
--- Template/NLO/SubProcesses/reweight_xsec_events.f	2015-08-19 12:50:52 +0000
+++ Template/NLO/SubProcesses/reweight_xsec_events.f	2016-03-16 17:59:03 +0000
@@ -61,7 +61,11 @@
       write(*,*)'      0 otherwise'
       read(*,*)isave
       if(isave.eq.1)then
-        isave=9
+        if (keep_rwgt_info)then
+           isave = -9
+        else
+           isave=9
+        endif
       else
         isave=0
       endif
@@ -431,10 +435,9 @@
 
         write(buff,201)'#aMCatNLO',iSorH_lhe,ifks_lhe,jfks_lhe,
      $       fksfather_lhe,ipartner_lhe, scale1_lhe,scale2_lhe, isave
-     $       ,izero,izero, wgtcentral,wgtmumin,wgtmumax,wgtpdfmin
+     $       ,mexternal,iwgtnumpartn, wgtcentral,wgtmumin,wgtmumax,wgtpdfmin
      $       ,wgtpdfmax
 
-
 c renormalize all the scale & PDF weights to have the same normalization
 c as XWGTUP
         if(do_rwgt_scale)then
@@ -494,7 +497,7 @@
       endif
       close(34)
 
- 201  format(a9,1x,i1,4(1x,i2),2(1x,e14.8),1x,i1,2(1x,i2),5(1x,e14.8))
+ 201  format(a9,1x,i1,4(1x,i2),2(1x,e14.8),1x,i2,2(1x,i2),5(1x,e14.8))
 
       end
 
@@ -574,16 +577,25 @@
       include 'nexternal.inc'
       include 'c_weight.inc'
       include 'reweight0.inc'
-      integer i,idum,j,k,momenta_conf
+      integer i,idum,j,k,momenta_conf(2),ii
       icontr=n_ctr_found
       iwgt=1
       do i=1,icontr
-         read(n_ctr_str(i),*)(wgt(j,i),j=1,3),idum,(pdg(j,i),j=1
-     &        ,nexternal),QCDpower(i),(bjx(j,i),j=1,2),(scales2(j,i),j=1
-     &        ,3),momenta_conf,itype(i),nFKS(i),wgts(1,i)
-         do j=1,nexternal
-            do k=0,3
-               momenta(k,j,i)=momenta_str(k,j,momenta_conf)
+         read(n_ctr_str(i),*)(wgt(j,i),j=1,3),(wgt_ME_tree(j,i),j=1,2)
+     &        ,idum,(pdg(j,i),j=1,nexternal),QCDpower(i),(bjx(j,i),j=1
+     &        ,2),(scales2(j,i),j=1,3),g_strong(i),(momenta_conf(j),j=1
+     &        ,2),itype(i),nFKS(i),idum,idum,idum,wgts(1,i)
+         do ii=1,2
+            do j=1,nexternal
+               do k=0,3
+                  if (momenta_conf(ii).gt.0) then
+                     momenta_m(k,j,ii,i)=momenta_str(k,j
+     $                                               ,momenta_conf(ii))
+                  else
+                     momenta_m(k,j,ii,i)=-99d0
+                     exit
+                  endif
+               enddo
             enddo
          enddo
       enddo

=== modified file 'Template/NLO/bin/aMCatNLO'
--- Template/NLO/bin/aMCatNLO	2015-03-14 02:21:38 +0000
+++ Template/NLO/bin/aMCatNLO	2016-03-16 17:59:03 +0000
@@ -164,7 +164,10 @@
             cmd_line.run_cmd('quit')
     elif args:
         # a single command is provided
-        cmd_line = cmd_interface.aMCatNLOCmd()
+        if options.web:
+            cmd_line = cmd_interface.aMCatNLOCmd()
+        else:
+            cmd_line = cmd_interface.aMCatNLOCmdShell()
         if not hasattr(cmd_line, 'do_%s' % args[0]):
             if parser_error:
                 print parser_error

=== added directory 'Template/RWGTNLO'
=== added file 'Template/RWGTNLO/__init__.py'
=== added file 'Template/RWGTNLO/alfas.inc'
--- Template/RWGTNLO/alfas.inc	1970-01-01 00:00:00 +0000
+++ Template/RWGTNLO/alfas.inc	2016-03-16 17:59:03 +0000
@@ -0,0 +1,11 @@
+c***********************************************************************
+c     this files contains the common blocks for the 
+c     the alpha_s settings
+c
+c     asmz   = alpha_s(Mz) is set based on the pdf chosen in setcuts.f
+c     nloop  = order of the running of alpha_s based on the pdf chosen  
+c***********************************************************************
+      integer nloop
+      double precision asmz
+      common/a_block/asmz,nloop
+

=== added file 'Template/RWGTNLO/alfas_functions_lhapdf.f'
--- Template/RWGTNLO/alfas_functions_lhapdf.f	1970-01-01 00:00:00 +0000
+++ Template/RWGTNLO/alfas_functions_lhapdf.f	2016-03-16 17:59:03 +0000
@@ -0,0 +1,167 @@
+C
+C-----------------------------------------------------------------------------
+C
+      double precision function alfa(alfa0,qsq )
+C
+C-----------------------------------------------------------------------------
+C
+C	This function returns the 1-loop value of alpha.
+C
+C	INPUT: 
+C		qsq   = Q^2
+C
+C-----------------------------------------------------------------------------
+C
+      implicit none
+      double precision  qsq,alfa0
+c
+c constants
+c
+      double precision  One, Three, Pi,zmass
+      parameter( One = 1.0d0, Three = 3.0d0 )
+      parameter( Pi = 3.14159265358979323846d0 )
+      parameter( zmass = 91.188d0 )
+cc
+      alfa = alfa0 / ( 1.0d0 - alfa0*dlog( qsq/zmass**2 ) /Three /Pi )
+ccc
+      return
+      end
+
+C
+C-----------------------------------------------------------------------------
+C
+      double precision function alfaw(alfaw0,qsq,nh )
+C
+C-----------------------------------------------------------------------------
+C
+C	This function returns the 1-loop value of alpha_w.
+C
+C	INPUT: 
+C		qsq = Q^2
+C               nh  = # of Higgs doublets
+C
+C-----------------------------------------------------------------------------
+C
+      implicit none
+      double precision  qsq, alphaw, dum,alfaw0
+      integer  nh, nq
+c
+c	  include
+c
+	  
+c
+c constants
+c
+      double precision  Two, Four, Pi, Twpi, zmass,tmass
+      parameter( Two = 2.0d0, Four = 4.0d0 )
+      parameter( Pi = 3.14159265358979323846d0 )
+      parameter( Twpi = 3.0d0*Four*Pi )
+      parameter( zmass = 91.188d0,tmass=174d0 )
+cc
+      if ( qsq.ge.tmass**2 ) then
+         nq = 6
+      else
+         nq = 5
+      end if
+      dum = ( 22.0d0 - Four*nq - nh/Two ) / Twpi
+      alfaw = alfaw0 / ( 1.0d0 + dum*alfaw0*dlog( qsq/zmass**2 ) )
+ccc
+      return
+      end
+
+C-----------------------------------------------------------------------------
+C
+      DOUBLE PRECISION FUNCTION ALPHAS(Q)
+C     wrapper to the lhapdf alphaS
+C-----------------------------------------------------------------------------
+      IMPLICIT NONE
+c
+      include 'alfas.inc'
+      REAL*8 Q,alphasPDF
+      external alphasPDF
+
+c timing statistics
+c      include "timing_variables.inc"
+
+c     This function takes 20 micro seconds to run, so it is ok to profile
+c     it with the 0.3 ms function cpu_time.
+c      call cpu_time(tbefore)
+      ALPHAS=alphasPDF(Q)
+c      call cpu_time(tAfter)
+      
+c      tPDF = tPDF + (tAfter-tBefore)
+
+      RETURN
+      END
+
+C-----------------------------------------------------------------------------
+C
+      double precision function mfrun(mf,scale,asmz,nloop)
+C
+C-----------------------------------------------------------------------------
+C
+C	This function returns the 2-loop value of a MSbar fermion mass
+C       at a given scale.
+C
+C	INPUT: mf    = MSbar mass of fermion at MSbar fermion mass scale 
+C	       scale = scale at which the running mass is evaluated
+C	       asmz  = AS(MZ) : this is passed to alphas(scale,asmz,nloop)
+C              nloop = # of loops in the evolution
+C       
+C
+C
+C	EXTERNAL:      double precision alphas(scale,asmz,nloop)
+C                      
+C-----------------------------------------------------------------------------
+C
+      implicit none
+C
+C     ARGUMENTS
+C
+      double precision  mf,scale,asmz
+      integer           nloop
+C
+C     LOCAL
+C
+      double precision  beta0, beta1,gamma0,gamma1
+      double precision  A1,as,asmf,l2
+      integer  nf
+C
+C     EXTERNAL
+C
+      double precision  alphas
+      external          alphas
+c
+c     CONSTANTS
+c
+      double precision  One, Two, Three, Pi
+      parameter( One = 1.0d0, Two = 2.0d0, Three = 3.0d0 )
+      parameter( Pi = 3.14159265358979323846d0) 
+      double precision tmass
+      parameter(tmass=174d0)
+cc
+C
+C
+      if ( mf.gt.tmass ) then
+         nf = 6
+      else
+         nf = 5
+      end if
+
+      beta0 = ( 11.0d0 - Two/Three *nf )/4d0
+      beta1 = ( 102d0  - 38d0/Three*nf )/16d0
+      gamma0= 1d0
+      gamma1= ( 202d0/3d0  - 20d0/9d0*nf )/16d0
+      A1    = -beta1*gamma0/beta0**2+gamma1/beta0
+      as    = alphas(scale)
+      asmf  = alphas(mf)
+      l2    = (1+ A1*as/Pi)/(1+ A1*asmf/Pi)
+      
+      
+      mfrun = mf * (as/asmf)**(gamma0/beta0)
+
+      if(nloop.eq.2) mfrun =mfrun*l2
+ccc
+      return
+      end
+

=== added file 'Template/RWGTNLO/makefile'
--- Template/RWGTNLO/makefile	1970-01-01 00:00:00 +0000
+++ Template/RWGTNLO/makefile	2016-03-16 17:59:03 +0000
@@ -0,0 +1,19 @@
+include ../Source/make_opts
+
+LIBDIR = ../lib/
+LINKLIBS =  -L../lib/ -lpdf $(llhapdf)
+LIBS = $(LIBDIR)libpdf.$(libext) 
+
+#note only support lhapdf=True
+RWGT=rwgt.f $(alphas_functions)
+
+rwgt2py.so: rwgt.f makefile $(LIBS)
+	touch __init__.py	
+	touch ../__init__.py
+	$(F2PY) $(LDFLAGS) $(LINKLIBS) -c  $(RWGT) -m rwgt2py --fcompiler=$(FC)
+
+$(LIBS):
+	cd PDF; make
+
+
+

=== added file 'Template/RWGTNLO/rwgt.f'
--- Template/RWGTNLO/rwgt.f	1970-01-01 00:00:00 +0000
+++ Template/RWGTNLO/rwgt.f	2016-03-16 17:59:03 +0000
@@ -0,0 +1,114 @@
+      subroutine initialise(lpp, pdfid)
+      implicit none
+      integer lpp(2)
+      integer pdfid
+CF2PY INTENT(IN) :: LPP(2)
+CF2PY INTENT(IN) :: PDF
+
+      integer save_lpp(2)
+      common/beam_type/save_lpp
+      double precision value(20)
+      character*20 parm(20)
+c       Save the beam type information        
+        save_lpp(1) = lpp(1)
+        save_lpp(2) = lpp(2)
+
+c       Initialise the pdf
+        if (abs(lpp(1)).eq.1.or.abs(lpp(2)).eq.1)then
+           value(1)=pdfid
+           parm(1)='DEFAULT'
+           call pdfset(parm,value)
+        endif
+        return
+        end
+
+      subroutine test_pdf()
+      implicit none 
+      integer lpp(2)
+      double precision pdf, PDG2PDF
+      lpp(1)=1
+      lpp(2) =1
+      call initialise(lpp, 230000)
+      pdf = PDG2PDF(1,1,0.1,100)
+      write(*,*) pdf
+      end
+     
+
+
+      subroutine get_wgt(scales2, pdg, bjx, wgt, g, qcdpower, ymur, 
+     &                   ymuf, n_ctr, new_wgt, all_wgt)
+      implicit none
+      double precision scales2(3,n_ctr)
+      integer pdg(2, n_ctr)
+      double precision bjx(2,n_ctr)
+      double precision g(n_ctr)
+      integer qcdpower(n_ctr)
+      double precision wgt(3, n_ctr)
+      double precision ymur, ymuf
+      integer n_ctr
+      double precision new_wgt
+      double precision all_wgt(n_ctr)
+
+CF2PY INTENT(IN) :: scales2(3, n_ctr)
+CF2PY INTENT(IN) :: pdg(2, n_ctr)
+CF2PY INTENT(IN) :: bjx(2, n_ctr)
+CF2PY INTENT(IN) :: wgt(3, n_ctr)
+CF2PY INTENT(IN) :: g(n_ctr)
+CF2PY INTENT(IN) :: qcdpower(n_ctr)
+CF2PY INTENT(IN) :: ymur
+CF2PY INTENT(IN) :: ymuf
+CF2PY depend(n_ctr) :: scales2, pdg, bjx, qcdpower
+CF2PY INTENT(OUT) :: new_wgt
+CF2PY INTENT(OUT) :: all_wgt(n_ctr)
+
+
+      integer save_lpp(2)
+      common/beam_type/save_lpp
+
+c     LOCAL      
+      integer i,kr
+      double precision mu2_q
+      double precision mu2_r, mu2_f
+      double precision xlum
+      integer lp, pd
+      double precision pi
+      parameter (pi=3.14159265358979323846d0)
+c     debug
+      double precision gs,add_wgt
+c     external
+      double precision alphas, pdg2pdf
+      external alphas, pdg2pdf
+      new_wgt = 0d0
+
+      do i=1, n_ctr
+         mu2_q=scales2(1,i)
+         mu2_r=scales2(2,i)*ymuR**2
+c Update the strong coupling
+         mu2_f=scales2(3,i)*ymuF**2
+c         gs = sqrt(4d0*pi*alphas(sqrt(mu2_r)))
+c call the PDFs
+         xlum=1d0
+         LP=SIGN(1,save_LPP(1))
+         pd=pdg(1,i)
+         if (pd.eq.21) pd=0
+         xlum=xlum*PDG2PDF(ABS(save_LPP(1)),pd*LP,bjx(1,i)
+     &           ,DSQRT(mu2_f))
+         LP=SIGN(1,save_LPP(2))
+         pd=pdg(2,i)
+         if (pd.eq.21) pd=0
+         xlum=xlum*PDG2PDF(ABS(save_LPP(2)),pd*LP,bjx(2,i)
+     &           ,DSQRT(mu2_f))
+c         write(*,*) "xlum", xlum
+c add the weights to the array
+         add_wgt =  xlum * (wgt(1,i)+wgt(2,i)*log(mu2_r
+     &              /mu2_q)+wgt(3,i)*log(mu2_f/mu2_q))*g(i)
+     &              **QCDpower(i)
+c         write(*,*) i, add_wgt
+         all_wgt(i) = add_wgt
+         new_wgt = new_wgt + add_wgt
+c         wgts(iwgt,i)=wgts(iwgt,i)
+c     &              *rwgt_muR_dep_fac(sqrt(mu2_r(kr)))
+      enddo
+      return
+      end
+

=== added file 'Template/RWGTNLO/setrun.f'
--- Template/RWGTNLO/setrun.f	1970-01-01 00:00:00 +0000
+++ Template/RWGTNLO/setrun.f	2016-03-16 17:59:03 +0000
@@ -0,0 +1,88 @@
+C-------------------------------------------------
+C   GET_PDFUP
+C   Convert MadEvent pdf name to LHAPDF number
+C-------------------------------------------------
+
+      subroutine get_pdfup(pdfin,pdfgup,pdfsup,lhaid)
+      implicit none
+
+      character*(*) pdfin
+      integer mpdf
+      integer npdfs,i,pdfgup(2),pdfsup(2),lhaid
+
+      parameter (npdfs=16)
+      character*7 pdflabs(npdfs)
+      data pdflabs/
+     $   'none',
+     $   'mrs02nl',
+     $   'mrs02nn',
+     $   'cteq4_m',
+     $   'cteq4_l',
+     $   'cteq4_d',
+     $   'cteq5_m',
+     $   'cteq5_d',
+     $   'cteq5_l',
+     $   'cteq5m1',
+     $   'cteq6_m',
+     $   'cteq6_l',
+     $   'cteq6l1',     
+     $   'nn23lo',
+     $   'nn23lo1',
+     $   'nn23nlo'/
+      integer numspdf(npdfs)
+      data numspdf/
+     $   00000,
+     $   20250,
+     $   20270,
+     $   19150,
+     $   19170,
+     $   19160,
+     $   19050,
+     $   19060,
+     $   19070,
+     $   19051,
+     $   10000,
+     $   10041,
+     $   10042,
+     $   246800,
+     $   247000,
+     $   244600/
+
+
+      if(pdfin.eq."lhapdf") then
+        write(*,*)'using LHAPDF'
+        do i=1,2
+           pdfgup(i)=-1
+           pdfsup(i)=lhaid
+        enddo
+        return
+      endif
+
+      
+      mpdf=-1
+      do i=1,npdfs
+        if(pdfin(1:len_trim(pdfin)) .eq. pdflabs(i))then
+          mpdf=numspdf(i)
+        endif
+      enddo
+
+      if(mpdf.eq.-1) then
+        write(*,*)'ERROR: pdf ',pdfin,' not implemented in get_pdfup.'
+        write(*,*)'known pdfs are'
+        write(*,*) pdflabs
+        open(unit=26,file='../../error',status='unknown')
+        write(26,*)'ERROR: pdf ',pdfin,' not implemented in get_pdfup.'
+        write(26,*)'known pdfs are'
+        write(26,*) pdflabs
+        stop 1
+c$$$        write(*,*)'using ',pdflabs(12)
+c$$$        mpdf=numspdf(12)
+      endif
+
+      do i=1,2
+        pdfgup(i)=-1
+        pdfsup(i)=mpdf
+      enddo
+
+      return
+      end

=== added file 'Template/RWGTNLO/timing_variables.inc'
--- Template/RWGTNLO/timing_variables.inc	1970-01-01 00:00:00 +0000
+++ Template/RWGTNLO/timing_variables.inc	2016-03-16 17:59:03 +0000
@@ -0,0 +1,9 @@
+*     -*-fortran-*-
+
+      real tBefore, tAfter
+      real tBorn,tIS,tReal,tCount,tFxFx,tf_nb,tf_all,t_as,tr_s,tr_pdf
+     &     ,t_plot,t_cuts,t_MC_subt,t_isum,t_p_unw,t_write
+      common/timings/tBorn,tIS,tReal,tCount,tFxFx,tf_nb,tf_all,t_as
+     &     ,tr_s,tr_pdf,t_plot,t_cuts,t_MC_subt,t_isum,t_p_unw,t_write
+      real tOLP, tFastJet, tPDF, tGenPS, tDSigI, tDSigR
+      common/timings_old/tOLP, tFastJet, tPDF, tGenPS, tDSigI, tDSigR

=== modified file 'Template/loop_material/StandAlone/SubProcesses/makefile'
--- Template/loop_material/StandAlone/SubProcesses/makefile	2016-02-12 01:48:16 +0000
+++ Template/loop_material/StandAlone/SubProcesses/makefile	2016-03-16 17:59:03 +0000
@@ -8,14 +8,15 @@
 else
     ROOT = ..
 endif
+LIBDIR = $(ROOT)/lib/
+
+PROG   = check
+all : $(PROG)
 
 include $(ROOT)/Source/make_opts
 include $(ROOT)/SubProcesses/MadLoop_makefile_definitions
 SHELL = /bin/bash
 
-LIBDIR = $(ROOT)/lib/
-PROG   = check
-all : $(PROG)
 OLP    = OLP
 STABCHECKDRIVER = StabilityCheckDriver
 CHECK_SA_BORN_SPLITORDERS = check_sa_born_splitOrders
@@ -80,19 +81,22 @@
 	ar rcs libMadLoop.$(libext) $(OLP_PROCESS)
 	mv libMadLoop.$(libext) $(MADLOOP_LIB)
 
-clean:
-	@rm -f *.o *.so
-
 ../$(OLP):
-	cd ..; make $(OLP); cd -
-	rm -f libMadLoop.dylib
-	ln -s ../libMadLoop.dylib
-
-libMadLoop.dylib: ../$(OLP)
+	rm -f libMadLoop.$(dylibext)
+	ln -s ../libMadLoop.$(dylibext)
+	cd $(ROOT)/SubProcesses; make $(OLP)
+
+../$(OLP)_static:
+	cd $(ROOT)/SubProcesses; make $(OLP)_static
+
+libMadLoop.$(dylibext): ../$(OLP)
 
 ifeq ($(origin MENUM),undefined)
   MENUM=2
 endif
-matrix$(MENUM)py.so: ../$(OLP)
+matrix$(MENUM)py.so: ../$(OLP)_static f2py_wrapper.f 
 	touch __init__.py
-	$(F2PY) -L. -lMadLoop -m matrix$(MENUM)py -c check_sa.f --fcompiler=$(FC)
+	$(F2PY) $(STDLIB) $(MADLOOP_LIB) -m matrix$(MENUM)py -c f2py_wrapper.f --fcompiler=$(FC)   -L../../lib/ -ldhelas -lmodel $(LINK_LOOP_LIBS)
+
+clean:
+	@rm -f *.o *.so *.$(libext) *.$(dylibext)

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2016-03-03 10:41:29 +0000
+++ UpdateNotes.txt	2016-03-16 17:59:03 +0000
@@ -1,7 +1,8 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
-2.4.0(XX/XX/XX)
-        MZ: updated FJcore to version 3.1.3 (was 3.0.5)
+
+2.4.0 (XX/XX/XX)
+        OM: Allowing the proper NLO reweighting for NLO sample
         MZ: new NLO generation mode. It is more efficient from the memory and CPU point of 
             view, in particular for high-multiplicity processes. 
             Many thanks to Josh Bendavid for his fundamental contribution for this.
@@ -24,6 +25,8 @@
         MZ: Fix problem with slurm cluster
         OM: Improve scan functionalities
         PT: New way of handling Pythia8 decays
+        OM: Allow to keep the reweight information in the final lhe file for future computation 
+        MZ: updated FJcore to version 3.1.3 (was 3.0.5)
 
 2.3.3 (15/10/15)
 	OM: Having two mode for "output pythia8" one (default) for pythia8.2 and one for pythia8.1 (with --version=8.1)

=== modified file 'VERSION'
--- VERSION	2015-10-25 15:53:09 +0000
+++ VERSION	2016-03-16 17:59:03 +0000
@@ -1,5 +1,5 @@
-version = 2.3.3
-date = 2015-10-25
+version = 2.4.0.beta
+date = 2016-01-25
 
 
 

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2016-02-18 20:27:51 +0000
+++ madgraph/core/base_objects.py	2016-03-16 17:59:03 +0000
@@ -2949,13 +2949,16 @@
                 'legs_with_decays', 'perturbation_couplings', 'has_born', 
                 'NLO_mode','split_orders']
 
-    def nice_string(self, indent=0, print_weighted = True):
+    def nice_string(self, indent=0, print_weighted = True, prefix=True):
         """Returns a nicely formated string about current process
         content. Since the WEIGHTED order is automatically set and added to 
         the user-defined list of orders, it can be ommitted for some info
         displays."""
 
-        mystr = " " * indent + "Process: "
+        if prefix:
+            mystr = " " * indent + "Process: "
+        else:
+            mystr = ""
         prevleg = None
         for leg in self['legs']:
             mypart = self['model'].get('particle_dict')[leg['id']]
@@ -3651,11 +3654,42 @@
 
         return max_order_now, particles, hierarchy
 
-    def nice_string(self, indent=0, print_weighted=False):
+    def __iter__(self):
+        """basic way to loop over all the process definition. 
+           not used by MG which used some smarter version (use by ML)"""
+        
+        isids = [leg['ids'] for leg in self['legs'] \
+                 if leg['state'] == False]
+        fsids = [leg['ids'] for leg in self['legs'] \
+                 if leg['state'] == True]
+
+        red_isidlist = []
+        # Generate all combinations for the initial state
+        for prod in itertools.product(*isids):
+            islegs = [Leg({'id':id, 'state': False}) for id in prod]  
+            if tuple(sorted(prod)) in red_isidlist:
+                    continue        
+            red_isidlist.append(tuple(sorted(prod)))      
+            red_fsidlist = []
+            for prod in itertools.product(*fsids):
+                # Remove double counting between final states
+                if tuple(sorted(prod)) in red_fsidlist:
+                    continue        
+                red_fsidlist.append(tuple(sorted(prod)))
+                leg_list = [copy.copy(leg) for leg in islegs]
+                leg_list.extend([Leg({'id':id, 'state': True}) for id in prod])
+                legs = LegList(leg_list)
+                process = self.get_process_with_legs(legs)
+                yield process
+
+    def nice_string(self, indent=0, print_weighted=False, prefix=True):
         """Returns a nicely formated string about current process
         content"""
 
-        mystr = " " * indent + "Process: "
+        if prefix:
+            mystr = " " * indent + "Process: "
+        else:
+            mystr=""
         prevleg = None
         for leg in self['legs']:
             myparts = \

=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2016-02-23 22:55:27 +0000
+++ madgraph/core/diagram_generation.py	2016-03-16 17:59:03 +0000
@@ -847,9 +847,6 @@
         for name, (value, operator) in self['process'].get('constrained_orders').items():
             res.filter_constrained_orders(name, value, operator)
             
-        
-        
-        
         # 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).

=== modified file 'madgraph/interface/amcatnlo_interface.py'
--- madgraph/interface/amcatnlo_interface.py	2016-02-24 13:44:29 +0000
+++ madgraph/interface/amcatnlo_interface.py	2016-03-16 17:59:03 +0000
@@ -737,7 +737,7 @@
             subproc_path = os.path.join(path, os.path.pardir, 'SubProcesses', \
                                      'initial_states_map.dat')
             self._curr_exporter.write_init_map(subproc_path,
-                                               self._curr_matrix_elements.get('initial_states'))
+                                self._curr_matrix_elements.get('initial_states'))
             
         cpu_time1 = time.time()
 

=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2016-03-14 10:17:10 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2016-03-16 17:59:03 +0000
@@ -938,16 +938,8 @@
         self.nb_core = 0
         self.prompt = "%s>"%os.path.basename(pjoin(self.me_dir))
 
-        # load the current status of the directory
-        if os.path.exists(pjoin(self.me_dir,'HTML','results.pkl')):
-            self.results = save_load_object.load_from_file(pjoin(self.me_dir,'HTML','results.pkl'))
-            self.results.resetall(self.me_dir)
-            self.last_mode = self.results[self.results.lastrun][-1]['run_mode']
-        else:
-            model = self.find_model_name()
-            process = self.process # define in find_model_name
-            self.results = gen_crossxhtml.AllResultsNLO(model, process, self.me_dir)
-            self.last_mode = ''
+
+        self.load_results_db()
         self.results.def_web_mode(self.web)
         # check that compiler is gfortran 4.6 or later if virtuals have been exported
         proc_card = open(pjoin(self.me_dir, 'Cards', 'proc_card_mg5.dat')).read()
@@ -1236,7 +1228,7 @@
 you have to remove some events after showering 'by hand'.
 Please read http://amcatnlo.cern.ch/FxFx_merging.htm for more details.""")
 
-
+        self.store_result()
         #check if the param_card defines a scan.
         if self.param_card_iterator:
             param_card_iterator = self.param_card_iterator
@@ -3324,6 +3316,17 @@
 
         if not self.to_store:
             return 
+
+        if 'event' in self.to_store:
+            if os.path.exists(pjoin(self.me_dir,'Events', self.run_name, 'events.lhe')):
+                if not  os.path.exists(pjoin(self.me_dir,'Events', self.run_name, 'events.lhe.gz')):
+                    self.update_status('gzipping output file: events.lhe', level='parton', error=True)
+                    misc.gzip(pjoin(self.me_dir,'Events', self.run_name, 'events.lhe'))
+                else:
+                    os.remove(pjoin(self.me_dir,'Events', self.run_name, 'events.lhe'))
+            if os.path.exists(pjoin(self.me_dir,'Events','reweight.lhe')):
+                os.remove(pjoin(self.me_dir,'Events', 'reweight.lhe'))
+                
         
         tag = self.run_card['run_tag']
         
@@ -4333,10 +4336,6 @@
                                 'madspin': default_switch,
                                 'reweight': default_switch}
 
-            
-            
-        
-        
         description = {'order':  'Perturbative order of the calculation:',
                        'fixed_order': 'Fixed order (no event generation and no MC@[N]LO matching):',
                        'shower': 'Shower the generated events:',
@@ -4539,9 +4538,17 @@
         if mode =='onlyshower':
             cards = ['shower_card.dat']
         
+        
+        # automatically switch to keep_wgt option
+        first_cmd = [] # force to change some switch
+        #if switch['reweight'] == 'ON':
+        #    logger.info("Automatically switch \"keep_rwgt_info\" in the run_card to allow NLO reweighting", '$MG:color:BLACK')
+        #    first_cmd.append('set keep_rwgt_info T')
+        
         if not options['force'] and not self.force:
-            self.ask_edit_cards(cards, plot=False)
+            self.ask_edit_cards(cards, plot=False, first_cmd=first_cmd)
 
+        
         self.banner = banner_mod.Banner()
 
         # store the cards in the banner
@@ -4590,7 +4597,6 @@
             analyse_card_path = pjoin(self.me_dir, 'Cards','FO_analyse_card.dat')
             self.analyse_card = self.banner.charge_card('FO_analyse_card')
 
-        
         return mode
 
 
@@ -4710,4 +4716,100 @@
 _shower_parser.add_option("-f", "--force", default=False, action='store_true',
                                 help="Use the shower_card present in the directory for the launch, without editing")
 
+if '__main__' == __name__:
+    # Launch the interface without any check if one code is already running.
+    # This can ONLY run a single command !!
+    import sys
+    if not sys.version_info[0] == 2 or sys.version_info[1] < 6:
+        sys.exit('MadGraph/MadEvent 5 works only with python 2.6 or later (but not python 3.X).\n'+\
+               'Please upgrate your version of python.')
+
+    import os
+    import optparse
+    # Get the directory of the script real path (bin)                                                                                                                                                           
+    # and add it to the current PYTHONPATH                                                                                                                                                                      
+    root_path = os.path.dirname(os.path.dirname(os.path.realpath( __file__ )))
+    sys.path.insert(0, root_path)
+
+    class MyOptParser(optparse.OptionParser):    
+        class InvalidOption(Exception): pass
+        def error(self, msg=''):
+            raise MyOptParser.InvalidOption(msg)
+    # Write out nice usage message if called with -h or --help                                                                                                                                                  
+    usage = "usage: %prog [options] [FILE] "
+    parser = MyOptParser(usage=usage)
+    parser.add_option("-l", "--logging", default='INFO',
+                      help="logging level (DEBUG|INFO|WARNING|ERROR|CRITICAL) [%default]")
+    parser.add_option("","--web", action="store_true", default=False, dest='web', \
+                     help='force toce to be in secure mode')
+    parser.add_option("","--debug", action="store_true", default=False, dest='debug', \
+                     help='force to launch debug mode')
+    parser_error = ''
+    done = False
+    
+    for i in range(len(sys.argv)-1):
+        try:
+            (options, args) = parser.parse_args(sys.argv[1:len(sys.argv)-i])
+            done = True
+        except MyOptParser.InvalidOption, error:
+            pass
+        else:
+            args += sys.argv[len(sys.argv)-i:]
+    if not done:
+        # raise correct error:                                                                                                                                                                                  
+        try:
+            (options, args) = parser.parse_args()
+        except MyOptParser.InvalidOption, error:
+            print error
+            sys.exit(2)
+
+    if len(args) == 0:
+        args = ''
+
+    import subprocess
+    import logging
+    import logging.config
+    # Set logging level according to the logging level given by options                                                                                                                                         
+    #logging.basicConfig(level=vars(logging)[options.logging])                                                                                                                                                  
+    import internal.coloring_logging
+    try:
+        if __debug__ and options.logging == 'INFO':
+            options.logging = 'DEBUG'
+        if options.logging.isdigit():
+            level = int(options.logging)
+        else:
+            level = eval('logging.' + options.logging)
+        print os.path.join(root_path, 'internal', 'me5_logging.conf')
+        logging.config.fileConfig(os.path.join(root_path, 'internal', 'me5_logging.conf'))
+        logging.root.setLevel(level)
+        logging.getLogger('madgraph').setLevel(level)
+    except:
+        raise
+        pass
+
+    # Call the cmd interface main loop                                                                                                                                                                          
+    try:
+        if args:
+            # a single command is provided   
+            if '--web' in args:
+                i = args.index('--web') 
+                args.pop(i)                                                                                                                                                                     
+                cmd_line =  aMCatNLOCmd(force_run=True)
+            else:
+                cmd_line =  aMCatNLOCmdShell(force_run=True)
+
+            if not hasattr(cmd_line, 'do_%s' % args[0]):
+                if parser_error:
+                    print parser_error
+                    print 'and %s  can not be interpreted as a valid command.' % args[0]
+                else:
+                    print 'ERROR: %s  not a valid command. Please retry' % args[0]
+            else:
+                cmd_line.use_rawinput = False
+                cmd_line.run_cmd(' '.join(args))
+                cmd_line.run_cmd('quit')
+
+    except KeyboardInterrupt:
+        print 'quit on KeyboardInterrupt'
+        pass
 

=== modified file 'madgraph/interface/coloring_logging.py'
--- madgraph/interface/coloring_logging.py	2015-10-01 16:00:08 +0000
+++ madgraph/interface/coloring_logging.py	2016-03-16 17:59:03 +0000
@@ -1,4 +1,7 @@
 import logging
+# method to add color to a logging.info add a second argument:
+# '$MG:color:BLACK'
+
 
 BLACK, RED, GREEN, YELLOW, BLUE, MAGENTA, CYAN, WHITE = range(8)
 

=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2016-02-21 00:08:42 +0000
+++ madgraph/interface/common_run_interface.py	2016-03-16 17:59:03 +0000
@@ -65,6 +65,8 @@
     import internal.cluster as cluster
     import internal.check_param_card as check_param_card
     import internal.files as files
+    import internal.save_load_object as save_load_object
+    import internal.gen_crossxhtml as gen_crossxhtml
     from internal import InvalidCmd, MadGraph5Error
     MADEVENT=True    
 else:
@@ -75,7 +77,10 @@
     import madgraph.various.misc as misc
     import madgraph.iolibs.files as files
     import madgraph.various.cluster as cluster
+    import madgraph.iolibs.save_load_object as save_load_object
+    import madgraph.madevent.gen_crossxhtml as gen_crossxhtml
     import models.check_param_card as check_param_card
+    
     from madgraph import InvalidCmd, MadGraph5Error, MG5DIR
     MADEVENT=False
 
@@ -391,7 +396,7 @@
 
     def check_decay_events(self,args):
         """Check the argument for decay_events command
-        syntax: decay_events [NAME]
+        syntax is "decay_events [NAME]"
         Note that other option are already remove at this point
         """
 
@@ -418,7 +423,7 @@
 
     def check_check_events(self,args):
         """Check the argument for decay_events command
-        syntax: decay_events [NAME]
+        syntax is "decay_events [NAME]"
         Note that other option are already remove at this point
         """
 
@@ -441,9 +446,7 @@
 
 
     def get_events_path(self, run_name):
-        """Check the argument for decay_events command
-        syntax: decay_events [NAME]
-        Note that other option are already remove at this point
+        """return the path to the output events
         """
 
 
@@ -523,6 +526,11 @@
     def __init__(self, me_dir, options, *args, **opts):
         """common"""
 
+        self.force_run = False # this flag force the run even if RunWeb is present
+        if 'force_run' in opts and opts['force_run']:
+            self.force_run = True
+            del opts['force_run']
+
         cmd.Cmd.__init__(self, *args, **opts)
         # Define current MadEvent directory
         if me_dir is None and MADEVENT:
@@ -538,21 +546,22 @@
         self.error =  pjoin(self.me_dir, 'error')
         self.dirbin = pjoin(self.me_dir, 'bin', 'internal')
 
-        # Check that the directory is not currently running
-        if os.path.exists(pjoin(me_dir,'RunWeb')): 
-            message = '''Another instance of the program is currently running.
-            (for this exact same directory) Please wait that this is instance is 
-            closed. If no instance is running, you can delete the file
-            %s and try again.''' % pjoin(me_dir,'RunWeb')
-            raise AlreadyRunning, message
-        else:
-            pid = os.getpid()
-            fsock = open(pjoin(me_dir,'RunWeb'),'w')
-            fsock.write(`pid`)
-            fsock.close()
-
-            misc.Popen([os.path.relpath(pjoin(self.dirbin, 'gen_cardhtml-pl'), me_dir)],
-                        cwd=me_dir)
+        # Check that the directory is not currently running_in_idle
+        if not self.force_run:
+            if os.path.exists(pjoin(me_dir,'RunWeb')): 
+                message = '''Another instance of the program is currently running.
+                (for this exact same directory) Please wait that this is instance is 
+                closed. If no instance is running, you can delete the file
+                %s and try again.''' % pjoin(me_dir,'RunWeb')
+                raise AlreadyRunning, message
+            else:
+                pid = os.getpid()
+                fsock = open(pjoin(me_dir,'RunWeb'),'w')
+                fsock.write(`pid`)
+                fsock.close()
+    
+                misc.Popen([os.path.relpath(pjoin(self.dirbin, 'gen_cardhtml-pl'), me_dir)],
+                            cwd=me_dir)
 
         self.to_store = []
         self.run_name = None
@@ -609,6 +618,35 @@
 
         return args
 
+    @misc.multiple_try(nb_try=5, sleep=2)
+    def load_results_db(self):
+        """load the current results status"""
+        
+        # load the current status of the directory
+        if os.path.exists(pjoin(self.me_dir,'HTML','results.pkl')):
+            try:
+                self.results = save_load_object.load_from_file(pjoin(self.me_dir,'HTML','results.pkl'))
+            except Exception:
+                #the pickle fail -> need to recreate the library
+                model = self.find_model_name()
+                process = self.process # define in find_model_name
+                self.results = gen_crossxhtml.AllResults(model, process, self.me_dir)
+                self.results.resetall(self.me_dir)
+            else:                                
+                self.results.resetall(self.me_dir)
+            try:
+                self.last_mode = self.results[self.results.lastrun][-1]['run_mode']
+            except:
+                self.results.resetall(self.me_dir)
+                self.last_mode = ''
+        else:
+            model = self.find_model_name()
+            process = self.process # define in find_model_name
+            self.results = gen_crossxhtml.AllResults(model, process, self.me_dir)
+            self.results.resetall(self.me_dir)
+            self.last_mode=''
+        return self.results
+
     ############################################################################
     def do_treatcards(self, line, amcatnlo=False):
         """Advanced commands: create .inc files from param_card.dat/run_card.dat"""
@@ -696,13 +734,13 @@
             param_card.write_inc_file(outfile, ident_card, default)
 
 
-    def ask_edit_cards(self, cards, mode='fixed', plot=True):
+    def ask_edit_cards(self, cards, mode='fixed', plot=True, first_cmd=None):
         """ """
         if not self.options['madanalysis_path']:
             plot = False
 
         self.ask_edit_card_static(cards, mode, plot, self.options['timeout'],
-                                  self.ask)
+                                  self.ask, first_cmd=first_cmd)
 
     @staticmethod
     def ask_edit_card_static(cards, mode='fixed', plot=True,
@@ -1024,10 +1062,13 @@
         """Dummy routine, to be overwritten by daughter classes"""
 
         pass
+    
     ############################################################################
     def do_reweight(self, line):
-        """ Allow to reweight the events generated with a new choices of model
-            parameter.
+        """ syntax is "reweight RUN_NAME"
+            Allow to reweight the events generated with a new choices of model
+            parameter. Description of the methods are available here
+            cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Reweight
         """
         
         if '-from_cards' in line and not os.path.exists(pjoin(self.me_dir, 'Cards', 'reweight_card.dat')):
@@ -1045,28 +1086,91 @@
             raise self.ConfigurationError, '''Can\'t load Reweight module.
             The variable mg5_path might not be correctly configured.'''
         
-        self.to_store.append('event')
+
+                        
         if not '-from_cards' in line:
             self.keep_cards(['reweight_card.dat'])
             self.ask_edit_cards(['reweight_card.dat'], 'fixed', plot=False)        
 
+        # load the name of the event file
+        args = self.split_arg(line) 
+
+        if not self.force_run:
+            # forbid this function to create an empty item in results.
+            if self.run_name and self.results.current and  self.results.current['cross'] == 0:
+                self.results.delete_run(self.run_name, self.run_tag)
+            self.results.save()
+            # we want to run this in a separate shell to avoid hard f2py crash
+            command =  [sys.executable]
+            if os.path.exists(pjoin(self.me_dir, 'bin', 'madevent')):
+                command.append(pjoin(self.me_dir, 'bin', 'internal','madevent_interface.py'))
+            else:
+                command.append(pjoin(self.me_dir, 'bin', 'internal', 'amcatnlo_run_interface.py'))
+            if not isinstance(self, cmd.CmdShell):
+                command.append('--web')
+            command.append('reweight')
+            if self.run_name:
+                command.append(self.run_name)
+            else:
+                command += args
+            if '-from_cards' not in command:
+                command.append('-from_cards')
+            p = misc.Popen(command, stdout = subprocess.PIPE, stderr = subprocess.STDOUT, cwd=self.me_dir)
+            while p.poll() is None:
+                line = p.stdout.readline()
+                if any(t in line for t in ['INFO:', 'WARNING:', 'CRITICAL:', 'ERROR:', 'root:','KEEP:']) and \
+                   not '***********' in line:
+                        print line[:-1].replace('INFO', 'REWEIGHT').replace('KEEP:','')
+                elif __debug__ and line:
+                    logger.debug(line[:-1])
+            if p.returncode !=0:
+                logger.error("Reweighting failed")
+                return
+            self.results = self.load_results_db()
+            # forbid this function to create an empty item in results.
+            try:
+                if self.results[self.run_name][-2]['cross']==0:
+                    self.results.delete_run(self.run_name,self.results[self.run_name][-2]['tag'])
+            except:
+                pass
+            try:
+                if self.results.current['cross'] == 0 and self.run_name:
+                    self.results.delete_run(self.run_name, self.run_tag)
+            except:
+                pass                    
+            # re-define current run     
+            try:
+                self.results.def_current(self.run_name, self.run_tag)
+            except Exception:
+                pass
+            return              
+
+        self.to_store.append('event')
         # forbid this function to create an empty item in results.
         if self.results.current['cross'] == 0 and self.run_name:
             self.results.delete_run(self.run_name, self.run_tag)
 
-        # load the name of the event file
-        args = self.split_arg(line) 
         self.check_decay_events(args) 
         # args now alway content the path to the valid files
-        reweight_cmd = reweight_interface.ReweightInterface(args[0])
-        reweight_cmd.mother = self
-        self.update_status('Running Reweight', level='madspin')
-        
+        reweight_cmd = reweight_interface.ReweightInterface(args[0], mother=self)
+        #reweight_cmd.use_rawinput = False
+        #reweight_cmd.mother = self
+        wgt_names = reweight_cmd.get_weight_names()
+        if wgt_names == [''] and reweight_cmd.has_nlo:
+            self.update_status('Running Reweighting (LO approximate)', level='madspin')
+        else:
+            self.update_status('Running Reweighting', level='madspin')
         
         path = pjoin(self.me_dir, 'Cards', 'reweight_card.dat')
+        reweight_cmd.raw_input=False
         reweight_cmd.me_dir = self.me_dir
         reweight_cmd.import_command_file(path)
         reweight_cmd.do_quit('')
+            
+        logger.info("quit rwgt")
+        
+        
+        
         # re-define current run
         try:
             self.results.def_current(self.run_name, self.run_tag)
@@ -1422,44 +1526,7 @@
                 else:
                     self.pdffile = " "
                 return self.pdffile
-                
-    def do_quit(self, line):
-        """Not in help: exit """
-  
-  
-        try:
-            os.remove(pjoin(self.me_dir,'RunWeb'))
-        except Exception, error:
-            pass
-        
-        try:
-            self.store_result()
-        except Exception:
-            # If nothing runs they they are no result to update
-            pass
-        
-        try:
-            self.update_status('', level=None)
-        except Exception, error:        
-            pass
-        try:
-            devnull = open(os.devnull, 'w') 
-            misc.call(['./bin/internal/gen_cardhtml-pl'], cwd=self.me_dir,
-                        stdout=devnull, stderr=devnull)
-        except Exception:
-            pass
-        try:
-            devnull.close()
-        except Exception:
-            pass
-        
-        return super(CommonRunCmd, self).do_quit(line)
-
-    
-    # Aliases
-    do_EOF = do_quit
-    do_exit = do_quit
-      
+                      
     ############################################################################
     def do_open(self, line):
         """Open a text file/ eps file / html file"""
@@ -1677,10 +1744,11 @@
     def do_quit(self, line):
         """Not in help: exit """
 
-        try:
-            os.remove(pjoin(self.me_dir,'RunWeb'))
-        except Exception:
-            pass
+        if not self.force_run:
+            try:
+                os.remove(pjoin(self.me_dir,'RunWeb'))
+            except Exception:
+                pass
         try:
             self.store_result()
         except Exception:
@@ -2039,6 +2107,20 @@
                                         os.path.join('.',*[a for a in args \
                                                     if a.endswith(os.path.sep)]))
 
+    def complete_reweight(self,text, line, begidx, endidx):
+        "Complete the pythia command"
+        args = self.split_arg(line[0:begidx], error=False)
+
+        #return valid run_name
+        data = glob.glob(pjoin(self.me_dir, 'Events', '*','*events.lhe*'))
+        data = [n.rsplit('/',2)[1] for n in data]
+        if not '-f' in args:
+            data.append('-f')
+        tmp1 =  self.list_completion(text, data)
+        return tmp1 
+
+
+
     def complete_compute_widths(self, text, line, begidx, endidx):
         "Complete the compute_widths command"
 
@@ -2066,20 +2148,57 @@
         """update the make_opts file writing the environmental variables
         stored in make_opts_var"""
         make_opts = os.path.join(self.me_dir, 'Source', 'make_opts')
+        
+        return self.update_make_opts_full(make_opts, self.make_opts_var)
+        
+
+    @staticmethod
+    def update_make_opts_full(path, def_variables, keep_old=True):
+        """update the make_opts file writing the environmental variables
+        of def_variables.
+        if a value of the dictionary is None then it is not written.
+        """
+        make_opts = path
+        pattern = re.compile(r'^(\w+)\s*=\s*(.*)$',re.DOTALL)
+        diff = False # set to True if one varible need to be updated 
+                     #if on False the file is not modify
+        
         tag = '#end_of_make_opts_variables\n'
-        content = open(make_opts).read()
-
-        # if this is not the first time that the file is updated, there
-        # should be a line #end_of_make_opts_variables
-        if tag in content:
-            content = content.split(tag)[1]
-
-        variables = '\n'.join('%s=%s' % (k,v) for k, v in self.make_opts_var.items())
-        variables += '\n%s' % tag
-
-        open(make_opts, 'w').write(variables + content)
-        return
-
+        make_opts_variable = True # flag to say if we are in edition area or not
+        content = []
+        variables = dict(def_variables)
+        need_keys = variables.keys()
+        for line in open(make_opts):
+            line = line.strip()
+            if make_opts_variable: 
+                if line.startswith('#') or not line:
+                    if line.startswith('#end_of_make_opts_variables'):
+                        make_opts_variable = False
+                    continue
+                elif pattern.search(line):
+                    key, value = pattern.search(line).groups()
+                    if key not in variables:
+                        variables[key] = value
+                    elif value !=  variables[key]:
+                        diff=True
+                    else:
+                        need_keys.remove(key)
+                else: 
+                    misc.sprint("end on line", line)
+                    make_opts_variable = False
+                    content.append(line)
+            else:                  
+                content.append(line)
+                     
+        if need_keys:
+            diff=True #This means that new definition are added to the file. 
+
+        content_variables = '\n'.join('%s=%s' % (k,v) for k, v in variables.items() if v is not None)
+        content_variables += '\n%s' % tag
+
+        if diff:
+            open(make_opts, 'w').write(content_variables + '\n'.join(content))
+        return       
 
 
 # lhapdf-related functions
@@ -2128,12 +2247,18 @@
         """copy (if needed) the lhapdf set corresponding to the lhaid in lhaid_list 
         into lib/PDFsets"""
 
+        if not hasattr(self, 'lhapdf_pdfsets'):
+            self.lhapdf_pdfsets = self.get_lhapdf_pdfsets_list(pdfsets_dir)
+
+
         pdfsetname = ''
         for lhaid in lhaid_list:
             try:
                 if not pdfsetname:
                     if lhaid in self.lhapdf_pdfsets:
                         pdfsetname = self.lhapdf_pdfsets[lhaid]['filename']
+                    elif isinstance(lhaid, str):
+                        pdfsetname = lhaid
                     else:
                         raise MadGraph5Error('lhaid %s not valid input number for the current lhapdf' % lhaid )
                 # just check the other ids refer to the same pdfsetname
@@ -2159,7 +2284,7 @@
                 os.mkdir(pdfsets_dir)
             except OSError:
                 pdfsets_dir = pjoin(self.me_dir, 'lib', 'PDFsets')
-        else:
+        elif os.path.exists(pjoin(self.me_dir, 'lib', 'PDFsets')):
             #clean previous set of pdf used
             for name in os.listdir(pjoin(self.me_dir, 'lib', 'PDFsets')):
                 if name != pdfsetname:
@@ -2213,18 +2338,49 @@
     def install_lhapdf_pdfset(self, pdfsets_dir, filename):
         """idownloads and install the pdfset filename in the pdfsets_dir"""
         lhapdf_version = self.get_lhapdf_version()
+        local_path = pjoin(self.me_dir, 'lib', 'PDFsets')
+        return self.install_lhapdf_pdfset_static(self.options['lhapdf'],
+                                             pdfsets_dir, filename,
+                                             lhapdf_version=lhapdf_version,
+                                             alternate_path=local_path)
+                                                 
+
+    @staticmethod
+    def install_lhapdf_pdfset_static(lhapdf_config, pdfsets_dir, filename, 
+                                        lhapdf_version=None, alternate_path=None):
+        """idownloads and install the pdfset filename in the pdfsets_dir.
+        Version which can be used independently of the class.
+        local path is used if the global installation fails.
+        """
+           
+        if not lhapdf_version:
+            lhapdf_version = subprocess.Popen([lhapdf_config, '--version'], 
+                        stdout = subprocess.PIPE).stdout.read().strip()
+        if not pdfsets_dir:
+            pdfsets_dir = subprocess.Popen([lhapdf_config, '--datadir'], 
+                        stdout = subprocess.PIPE).stdout.read().strip()
+                                
+        if isinstance(filename, int):
+            pdf_info = CommonRunCmd.get_lhapdf_pdfsets_list_static(pdfsets_dir, lhapdf_version)
+            filename = pdf_info[filename]['filename']
+        
+        if os.path.exists(pjoin(pdfsets_dir, filename)):
+            logger.debug('%s is already present in %s', (filename, pdfsets_dir))
+            return
+             
         logger.info('Trying to download %s' % filename)
 
         if lhapdf_version.startswith('5.'):
+
             # use the lhapdf-getdata command, which is in the same path as
             # lhapdf-config
-            getdata = self.options['lhapdf'].replace('lhapdf-config', ('lhapdf-getdata'))
+            getdata = lhapdf_config.replace('lhapdf-config', ('lhapdf-getdata'))
             misc.call([getdata, filename], cwd = pdfsets_dir)
 
         elif lhapdf_version.startswith('6.'):
             # use the "lhapdf install xxx" command, which is in the same path as
             # lhapdf-config
-            getdata = self.options['lhapdf'].replace('lhapdf-config', ('lhapdf'))
+            getdata = lhapdf_config.replace('lhapdf-config', ('lhapdf'))
 
             misc.call([getdata, 'install', filename], cwd = pdfsets_dir)
 
@@ -2240,18 +2396,24 @@
         elif lhapdf_version.startswith('5.'):
             logger.warning('Could not download %s into %s. Trying to save it locally' \
                     % (filename, pdfsets_dir))
-            self.install_lhapdf_pdfset(pjoin(self.me_dir, 'lib', 'PDFsets'), filename)
+            CommonRunCmd.install_lhapdf_pdfset_static(lhapdf_config, alternate_path, filename,
+                                                      lhapdf_version=lhapdf_version)
         else:
             raise MadGraph5Error, \
                 'Could not download %s into %s. Please try to install it manually.' \
                     % (filename, pdfsets_dir)
 
 
+
     def get_lhapdf_pdfsets_list(self, pdfsets_dir):
         """read the PDFsets.index file, which should be located in the same
         place as pdfsets_dir, and return a list of dictionaries with the information
         about each pdf set"""
         lhapdf_version = self.get_lhapdf_version()
+        return self.get_lhapdf_pdfsets_list_static(pdfsets_dir, lhapdf_version)
+
+    @staticmethod
+    def get_lhapdf_pdfsets_list_static(pdfsets_dir, lhapdf_version):
 
         if lhapdf_version.startswith('5.'):
             if os.path.exists('%s.index' % pdfsets_dir):
@@ -2319,6 +2481,19 @@
 
         return datadir
 
+    def get_lhapdf_libdir(self):
+        lhapdf_version = self.get_lhapdf_version()
+
+        if lhapdf_version.startswith('5.'):
+            libdir = subprocess.Popen([self.options['lhapdf-config'], '--libdir'],
+                         stdout = subprocess.PIPE).stdout.read().strip()
+
+        elif lhapdf_version.startswith('6.'):
+            libdir = subprocess.Popen([self.options['lhapdf'], '--libs'],
+                         stdout = subprocess.PIPE).stdout.read().strip()
+
+        return libdir
+
 class AskforEditCard(cmd.OneLinePathCompletion):
     """A class for asking a question where in addition you can have the
     set command define and modifying the param_card/run_card correctly"""
@@ -3082,12 +3257,39 @@
                     (block, lhaid, value), '$MG:color:BLACK')
         self.param_card[block].param_dict[lhaid].value = value
     
+    def check_card_consistency(self):
+        """This is run on quitting the class. Apply here all the self-consistency
+        rule that you want. Do the modification via the set command."""
+
+        # if NLO reweighting is ON: ensure that we keep the rwgt information
+        if 'reweight' in self.allow_arg and 'run' in self.allow_arg and \
+            isinstance(self.run_card,banner_mod.RunCardNLO) and \
+            not self.run_card['keep_rwgt_info']:
+            #check if a NLO reweighting is required
+                re_pattern = re.compile(r'''^\s*change\s*mode\s* (LO\+NLO|LO|NLO)\s*(?:#|$)''', re.M+re.I)
+                text = open(pjoin(self.me_dir,'Cards','reweight_card.dat')).read()
+                options = re_pattern.findall(text)
+                if any(o in ['NLO', 'LO+NLO'] for o in options):
+                    logger.info('NLO reweighting is on ON. Automatically set keep_rwgt_info to True', '$MG:color:BLACK' )
+                    self.do_set('run_card keep_rwgt_info True')
+    
+    
     def reask(self, *args, **opt):
         
         cmd.OneLinePathCompletion.reask(self,*args, **opt)
         if self.has_mw and not os.path.exists(pjoin(self.me_dir,'Cards','transfer_card.dat')):
             logger.warning('No transfer function currently define. Please use the change_tf command to define one.')
-            
+    
+    def postcmd(self, stop, line):
+        
+        ending_question = cmd.OneLinePathCompletion.postcmd(self,stop,line)
+        if ending_question:
+            self.check_card_consistency()
+            return ending_question
+    
+    def check_answer_consistency(self):
+        """function called if the code reads a file"""
+        self.check_card_consistency() 
       
     def help_set(self):
         '''help message for set'''

=== modified file 'madgraph/interface/extended_cmd.py'
--- madgraph/interface/extended_cmd.py	2016-03-02 12:15:58 +0000
+++ madgraph/interface/extended_cmd.py	2016-03-16 17:59:03 +0000
@@ -692,8 +692,11 @@
                                                    mother_interface=self, **opt)
         
         if first_cmd:
-            question_instance.onecmd(first_cmd)
-        
+            if isinstance(first_cmd, str):
+                question_instance.onecmd(first_cmd)
+            else:
+                for line in first_cmd:
+                    question_instance.onecmd(line)
         if not self.haspiping:
             if hasattr(obj, "haspiping"):
                 obj.haspiping = self.haspiping
@@ -707,13 +710,13 @@
                 answer = alias[answer]
             if ask_class:
                 answer = question_instance.default(answer)
+            if hasattr(question_instance, 'check_answer_consistency'):
+                question_instance.check_answer_consistency()
             return answer
         
         question = question_instance.question
         value =   Cmd.timed_input(question, default, timeout=timeout,
                                  fct=question_instance, fct_timeout=fct_timeout)
-
-
         
         try:
             if value in alias:
@@ -765,8 +768,8 @@
                     self.store_line(line)
                     return None # print the question and use the pipe
                 logger.info(question_instance.question)
-                logger.warning('The answer to the previous question is not set in your input file')
-                logger.warning('Use %s value' % default)
+                logger.info('The answer to the previous question is not set in your input file', '$MG:color:BLACK')
+                logger.info('Use %s value' % default, '$MG:color:BLACK')
                 return str(default)
         
         line = line.replace('\n','').strip()

=== modified file 'madgraph/interface/loop_interface.py'
--- madgraph/interface/loop_interface.py	2016-03-04 16:16:23 +0000
+++ madgraph/interface/loop_interface.py	2016-03-16 17:59:03 +0000
@@ -539,24 +539,16 @@
                             me, self._curr_fortran_model)
             # If all ME's do not share the same maximum loop vertex rank and the
             # same loop maximum wavefunction size, we need to set the maximum
-            # in coef_specs.inc of the HELAS Source and warn the user that this
-            # might be a problem
+            # in coef_specs.inc of the HELAS Source. The SubProcesses/P* directory
+            # all link this file, so it should be properly propagated
             if self.options['loop_optimized_output'] and len(matrix_elements)>1:
                 max_lwfspins = [m.get_max_loop_particle_spin() for m in \
                                                                 matrix_elements]
-                try:
-                    max_loop_vert_ranks = [me.get_max_loop_vertex_rank() for me in \
+                max_loop_vert_ranks = [me.get_max_loop_vertex_rank() for me in \
                                                                 matrix_elements]
-                except MadGraph5Error:
-                    pass
-                else:
-                    if len(set(max_lwfspins))>1 or len(set(max_loop_vert_ranks))>1:
-                        self._curr_exporter.fix_coef_specs(max(max_lwfspins),\
+                if len(set(max_lwfspins))>1 or len(set(max_loop_vert_ranks))>1:
+                    self._curr_exporter.fix_coef_specs(max(max_lwfspins),\
                                                        max(max_loop_vert_ranks))
-                        logger.warning('ML5 has just output processes which do not'+\
-                      ' share the same maximum loop wavefunction size or the '+\
-                      ' same maximum loop vertex rank. This is potentially '+\
-                      ' dangerous. Please prefer to output them separately.')
 
         # Just the matrix.f files
         if self._export_format == 'matrix':
@@ -719,6 +711,7 @@
                 if not isinstance(self, extended_cmd.CmdShell):
                     raise InvalidCmd, "loop_filter is not allowed in web mode"
             args = [a for a in args if not a.startswith('--loop_filter=')]
+
             # Rejoin line
             line = ' '.join(args[1:])
             
@@ -731,6 +724,22 @@
             
         # Extract process from process definition
         myprocdef = self.extract_process(line)
+        # hack for multiprocess:
+        if myprocdef.has_multiparticle_label():
+            # split it in a loop
+            succes, failed = 0, 0
+            for base_proc in myprocdef:
+                try:
+                    self.exec_cmd("add process %s" % base_proc.nice_string(prefix=False, print_weighted=True))
+                    succes += 1
+                except Exception:
+                    failed +=1
+            logger.info("%s/%s processes succeeded" % (succes, failed+succes))
+            if succes == 0:
+                raise
+            else:
+                return
+             
              
         # If it is a process for MadLoop standalone, make sure it has a 
         # unique ID. It is important for building a BLHA library which

=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2016-03-02 17:48:27 +0000
+++ madgraph/interface/madevent_interface.py	2016-03-16 17:59:03 +0000
@@ -938,7 +938,7 @@
     
     def check_pythia(self, args):
         """Check the argument for pythia command
-        syntax: pythia [NAME] 
+        syntax is  "pythia [NAME]" 
         Note that other option are already remove at this point
         """
         
@@ -1124,7 +1124,7 @@
     
     def check_pgs(self, arg, no_default=False):
         """Check the argument for pythia command
-        syntax: pgs [NAME] 
+        syntax is  "pgs [NAME]" 
         Note that other option are already remove at this point
         """
         
@@ -1182,7 +1182,7 @@
 
     def check_delphes(self, arg):
         """Check the argument for pythia command
-        syntax: delphes [NAME] 
+        syntax is "delphes [NAME]" 
         Note that other option are already remove at this point
         """
         
@@ -1239,7 +1239,7 @@
 
     def check_display(self, args):
         """check the validity of line
-        syntax: display XXXXX
+        syntax is "display XXXXX"
         """
             
         if len(args) < 1 or args[0] not in self._display_opts:
@@ -1656,24 +1656,7 @@
         if self.web:
             os.system('touch %s' % pjoin(self.me_dir,'Online'))
 
-        
-        # load the current status of the directory
-        if os.path.exists(pjoin(self.me_dir,'HTML','results.pkl')):
-            try:
-                self.results = save_load_object.load_from_file(pjoin(self.me_dir,'HTML','results.pkl'))
-            except Exception:
-                #the pickle fail -> need to recreate the library
-                model = self.find_model_name()
-                process = self.process # define in find_model_name
-                self.results = gen_crossxhtml.AllResults(model, process, self.me_dir)
-                self.results.resetall(self.me_dir)
-            else:                                
-                self.results.resetall(self.me_dir)
-        else:
-            model = self.find_model_name()
-            process = self.process # define in find_model_name
-            self.results = gen_crossxhtml.AllResults(model, process, self.me_dir)
-            self.results.resetall(self.me_dir)
+        self.load_results_db()        
         self.results.def_web_mode(self.web)
         
         self.prompt = "%s>"%os.path.basename(pjoin(self.me_dir))
@@ -3540,7 +3523,17 @@
         
         self.ask_edit_cards(['run_card'], args)
         self.run_card = banner_mod.RunCard(pjoin(self.me_dir, 'Cards', 'run_card.dat'))
-                
+        
+        # Check that all pdfset are correctly installed
+        lhaid = [self.run_card.get_lhapdf_id()]
+        sys_pdf = self.run_card['sys_pdf'].split('&&')
+        lhaid += [l.split()[0] for l in sys_pdf]
+        pdfsets_dir = self.get_lhapdf_pdfsetsdir()
+        # Copy all the relevant PDF sets
+        [self.copy_lhapdf_set([onelha], pdfsets_dir) for onelha in lhaid]
+        
+        
+            
         if any([arg in ['all','parton'] for arg in args]):
             filename = pjoin(self.me_dir, 'Events', self.run_name, 'unweighted_events.lhe')
             if os.path.exists(filename+'.gz'):
@@ -3593,7 +3586,10 @@
         if 'event' in self.to_store:
             if not os.path.exists(pjoin(self.me_dir, 'Events',self.run_name, 'unweighted_events.lhe.gz')) and\
                os.path.exists(pjoin(self.me_dir, 'Events',self.run_name, 'unweighted_events.lhe')):
+                logger.info("gzipping output file: unweighted_events.lhe")
                 misc.gzip(pjoin(self.me_dir,'Events',self.run_name,"unweighted_events.lhe"))
+            if os.path.exists(pjoin(self.me_dir,'Events','reweight.lhe')):
+                os.remove(pjoin(self.me_dir,'Events', 'reweight.lhe'))
         
         if 'pythia' in self.to_store:
             self.update_status('Storing Pythia files of previous run', level='pythia', error=True)
@@ -5098,3 +5094,111 @@
 
 AskforEditCard = common_run.AskforEditCard
 
+
+if '__main__' == __name__:
+    # Launch the interface without any check if one code is already running.
+    # This can ONLY run a single command !!
+    import sys
+    if not sys.version_info[0] == 2 or sys.version_info[1] < 6:
+        sys.exit('MadGraph/MadEvent 5 works only with python 2.6 or later (but not python 3.X).\n'+\
+               'Please upgrate your version of python.')
+
+    import os
+    import optparse
+    # Get the directory of the script real path (bin)                                                                                                                                                           
+    # and add it to the current PYTHONPATH                                                                                                                                                                      
+    root_path = os.path.dirname(os.path.dirname(os.path.realpath( __file__ )))
+    sys.path.insert(0, root_path)
+
+    class MyOptParser(optparse.OptionParser):    
+        class InvalidOption(Exception): pass
+        def error(self, msg=''):
+            raise MyOptParser.InvalidOption(msg)
+    # Write out nice usage message if called with -h or --help                                                                                                                                                  
+    usage = "usage: %prog [options] [FILE] "
+    parser = MyOptParser(usage=usage)
+    parser.add_option("-l", "--logging", default='INFO',
+                      help="logging level (DEBUG|INFO|WARNING|ERROR|CRITICAL) [%default]")
+    parser.add_option("","--web", action="store_true", default=False, dest='web', \
+                     help='force toce to be in secure mode')
+    parser.add_option("","--debug", action="store_true", default=False, dest='debug', \
+                     help='force to launch debug mode')
+    parser_error = ''
+    done = False
+    
+    for i in range(len(sys.argv)-1):
+        try:
+            (options, args) = parser.parse_args(sys.argv[1:len(sys.argv)-i])
+            done = True
+        except MyOptParser.InvalidOption, error:
+            pass
+        else:
+            args += sys.argv[len(sys.argv)-i:]
+    if not done:
+        # raise correct error:                                                                                                                                                                                  
+        try:
+            (options, args) = parser.parse_args()
+        except MyOptParser.InvalidOption, error:
+            print error
+            sys.exit(2)
+
+    if len(args) == 0:
+        args = ''
+
+    import subprocess
+    import logging
+    import logging.config
+    # Set logging level according to the logging level given by options                                                                                                                                         
+    #logging.basicConfig(level=vars(logging)[options.logging])                                                                                                                                                  
+    import internal.coloring_logging
+    try:
+        if __debug__ and options.logging == 'INFO':
+            options.logging = 'DEBUG'
+        if options.logging.isdigit():
+            level = int(options.logging)
+        else:
+            level = eval('logging.' + options.logging)
+        print os.path.join(root_path, 'internal', 'me5_logging.conf')
+        logging.config.fileConfig(os.path.join(root_path, 'internal', 'me5_logging.conf'))
+        logging.root.setLevel(level)
+        logging.getLogger('madgraph').setLevel(level)
+    except:
+        raise
+        pass
+
+    # Call the cmd interface main loop                                                                                                                                                                          
+    try:
+        if args:
+            # a single command is provided
+            if '--web' in args:
+                i = args.index('--web') 
+                args.pop(i)                                                                                                                                                                     
+                cmd_line = MadEventCmd(force_run=True)
+            else:
+                cmd_line = MadEventCmdShell(force_run=True)
+            if not hasattr(cmd_line, 'do_%s' % args[0]):
+                if parser_error:
+                    print parser_error
+                    print 'and %s  can not be interpreted as a valid command.' % args[0]
+                else:
+                    print 'ERROR: %s  not a valid command. Please retry' % args[0]
+            else:
+                cmd_line.use_rawinput = False
+                cmd_line.run_cmd(' '.join(args))
+                cmd_line.run_cmd('quit')
+
+    except KeyboardInterrupt:
+        print 'quit on KeyboardInterrupt'
+        pass
+
+
+    
+    
+    
+    
+    
+    
+    
+    
+
+

=== modified file 'madgraph/interface/reweight_interface.py'
--- madgraph/interface/reweight_interface.py	2016-01-12 14:46:58 +0000
+++ madgraph/interface/reweight_interface.py	2016-03-16 17:59:03 +0000
@@ -33,6 +33,7 @@
 import madgraph.interface.madgraph_interface as mg_interface
 import madgraph.interface.master_interface as master_interface
 import madgraph.interface.common_run_interface as common_run_interface
+import madgraph.interface.madevent_interface as madevent_interface
 import madgraph.iolibs.files as files
 import MadSpin.interface_madspin as madspin_interface
 import madgraph.various.misc as misc
@@ -67,7 +68,7 @@
     debug_output = 'Reweight_debug'
     
     @misc.mute_logger()
-    def __init__(self, event_path=None, allow_madspin=False, *completekey, **stdin):
+    def __init__(self, event_path=None, allow_madspin=False, mother=None, *completekey, **stdin):
         """initialize the interface with potentially an event_path"""
         
         if not event_path:
@@ -80,8 +81,11 @@
         
         self.model = None
         self.has_standalone_dir = False
+        self.mother= mother # calling interface
+            
         
-        self.options = {'curr_dir': os.path.realpath(os.getcwd())}
+        self.options = {'curr_dir': os.path.realpath(os.getcwd()),
+                        'rwgt_name':None}
         
         self.events_file = None
         self.processes = {}
@@ -91,8 +95,14 @@
         self.seed = None
         self.output_type = "default"
         self.helicity_reweighting = True
+        self.rwgt_mode = None # can be LO or NLO, None is default 
+        self.has_nlo = False
         self.rwgt_dir = None
         self.exitted = False # Flag to know if do_quit was already called.
+        
+        if 'lhapdf' in self.mother.options and not mother.options['lhapdf']:
+            logger.warning('NLO accurate reweighting requires lhapdf to be installed. Pass in approximate LO mode.')
+            self.rwgt_mode = 'LO'
         if event_path:
             logger.info("Extracting the banner ...")
             self.do_import(event_path, allow_madspin=allow_madspin)
@@ -155,7 +165,6 @@
         
         # Check the validity of the banner:
         if 'slha' not in self.banner:
-            misc.sprint(self.banner)
             self.events_file = None
             raise self.InvalidCmd('Event file does not contain model information')
         elif 'mg5proccard' not in self.banner:
@@ -169,8 +178,16 @@
         # load information
         process = self.banner.get_detail('proc_card', 'generate')
         if '[' in process:
-            logger.warning("Remember that the reweighting is performed at Leading Order. NLO precision is not guarantee.")
-
+            if not self.banner.get_detail('run_card', 'keep_rwgt_info'):
+                logger.warning("The information to perform a proper NLO reweighting is not present in the event file.")
+                logger.warning("       We will perform a LO reweighting instead. This does not guarantee NLO precision.")
+                self.rwgt_mode = 'LO'
+            
+            if 'OLP' in self.mother.options:
+                if self.mother.options['OLP'].lower() != 'madloop':
+                    logger.warning("Accurate NLO mode only works for OLP=MadLoop not for OLP=%s. An approximate (LO) reweighting will be performed instead")
+                    self.rwgt_mode = 'LO'
+                
         if not process:
             msg = 'Invalid proc_card information in the file (no generate line):\n %s' % self.banner['mg5proccard']
             raise Exception, msg
@@ -191,7 +208,7 @@
         if not order:
             #NO NLO tag => nothing to do actually return input
             return proc
-        elif not order.startswith(('virt=','loonly=','noborn=')):
+        elif not order.startswith(('virt','loonly','noborn')):
             # OK this a standard NLO process
             if '=' in order:
                 # get the type NLO QCD/QED/...
@@ -212,7 +229,7 @@
                         r='QCD=%i' % (int(ior[1])+1)
                     process=process+r+' '
             #handle special tag $ | / @
-            result = re.split('([/$@]|\w+=\w+)', process, 1)                    
+            result = re.split('([/$@]|\w+(?:^2)?(?:=|<=|>)?\w+)', process, 1)                    
             if len(result) ==3:
                 process, split, rest = result
                 commandline+="add process %s pert_%s %s%s %s --no_warning=duplicate;" % (process, order.replace(' ','') ,split, rest, final)
@@ -220,11 +237,11 @@
                 commandline +='add process %s pert_%s %s --no_warning=duplicate;' % (process,order.replace(' ',''), final)
         elif order.startswith(('noborn=')):
             # pass in sqrvirt=
-            return "add process %s " % proc.replace('noborn=', 'sqrvirt=')
+            return "add process %s ;" % proc.replace('noborn=', 'sqrvirt=')
             
         else:
             #just return the input. since this Madloop.
-            return "add process %s " % proc                                       
+            return "add process %s ;" % proc                                       
         return commandline
 
 
@@ -319,6 +336,21 @@
                 self.output_type = args[1]
         elif args[0] == "helicity":
             self.helicity_reweighting = banner.ConfigFile.format_variable(args[1], bool, "helicity")
+        elif args[0] == "mode":
+            if args[1] != 'LO':
+                if 'OLP' in self.mother.options and self.mother.options['OLP'].lower() != 'madloop':
+                    logger.warning("Only LO reweighting is allowed for OLP!=MadLoop. Keeping the mode to LO.")
+                    self.rwgt_mode = 'LO'
+                elif not self.banner.get_detail('run_card','keep_rwgt_info', default=False):
+                    logger.warning("Missing information for NLO type of reweighting. Keeping the mode to LO.")
+                    self.rwgt_mode = 'LO'
+                elif 'lhapdf' in self.mother.options and not self.mother.options['lhapdf']:
+                    logger.warning('NLO accurate reweighting requires lhapdf to be installed. Pass in approximate LO mode.')
+                    self.rwgt_mode = 'LO'
+                else:
+                    self.rwgt_mode = args[1]
+            else:
+                self.rwgt_mode = args[1]
         elif args[0] == "rwgt_dir":
             self.rwgt_dir = args[1]
             if not os.path.exists(self.rwgt_dir):
@@ -335,6 +367,14 @@
                 self.lhe_input = lhe_parser.EventFile(self.lhe_input.name)
             else:
                 raise self.InvalidCmd("No events files defined.")
+            
+        opts = {'rwgt_name':None}
+        if any(a.startswith('--') for a in args):
+            for a in args[:]:
+                if a.startswith('--') and '=' in a:
+                    key,value = a[2:].split('=')
+                    opts[key] = value .replace("'","") .replace('"','')
+        return opts
 
     def help_launch(self):
         """help for the launch command"""
@@ -345,12 +385,30 @@
         All scale are kept fix for this re-weighting.''')
 
 
-    #@misc.mute_logger()
+    def get_weight_names(self):
+        """ return the various name for the computed weights """
+        
+        if self.rwgt_mode == 'LO':
+            return ['']
+        elif self.rwgt_mode == 'NLO':
+            return ['_nlo']
+        elif self.rwgt_mode == 'LO+NLO':
+            return ['_lo', '_nlo']
+        elif self.rwgt_mode == 'NLO_tree':
+            return ['_tree']        
+        elif not self.rwgt_mode and self.has_nlo :
+            return ['_nlo']
+        else:
+            return ['']
+
+    @misc.mute_logger()
     def do_launch(self, line):
         """end of the configuration launched the code"""
         
         args = self.split_arg(line)
-        self.check_launch(args)
+        opts = self.check_launch(args)
+        if opts['rwgt_name']:
+            self.options['rwgt_name'] = opts['rwgt_name']
 
         model_line = self.banner.get('proc_card', 'full_model_line')
 
@@ -370,28 +428,43 @@
             rw_dir = pjoin(path_me, 'rw_me_second')
         else:
             rw_dir = pjoin(path_me, 'rw_me')
-        
+
         if not '--keep_card' in args:
             ff = open(pjoin(rw_dir,'Cards', 'param_card.dat'), 'w')
             ff.write(self.banner['slha'])
             ff.close()
+            if self.has_nlo and self.rwgt_mode != "LO":
+                rwdir_virt = rw_dir.replace('rw_me', 'rw_mevirt')
+                files.ln(ff.name, starting_dir=pjoin(rwdir_virt, 'Cards')) 
             ff = open(pjoin(path_me, 'rw_me','Cards', 'param_card_orig.dat'), 'w')
             ff.write(self.banner['slha'])
             ff.close()      
+            if self.has_nlo and self.rwgt_mode != "LO":
+                files.ln(ff.name, starting_dir=pjoin(path_me, 'rw_mevirt', 'Cards'))
             cmd = common_run_interface.CommonRunCmd.ask_edit_card_static(cards=['param_card.dat'],
                                    ask=self.ask, pwd=rw_dir, first_cmd=self.stored_line)
             self.stored_line = None
         
+        # get the names of type of reweighting requested
+        type_rwgt = self.get_weight_names()
+
         # check for potential scan in the new card 
         new_card = open(pjoin(rw_dir, 'Cards', 'param_card.dat')).read()
         pattern_scan = re.compile(r'''^[\s\d]*scan''', re.I+re.M) 
         param_card_iterator = []
         if pattern_scan.search(new_card):
+            try:
+                import internal.extended_cmd as extended_internal
+                Shell_internal = extended_internal.CmdShell
+            except:
+                Shell_internal = None
             import madgraph.interface.extended_cmd as extended_cmd
-            if not isinstance(self.mother, extended_cmd.CmdShell): 
+            if not isinstance(self.mother, (extended_cmd.CmdShell, Shell_internal)): 
                 raise Exception, "scan are not allowed on the Web"
             # at least one scan parameter found. create an iterator to go trough the cards
             main_card = check_param_card.ParamCardIterator(new_card)
+            if self.options['rwgt_name']:
+                self.options['rwgt_name'] = '%s_0' % self.options['rwgt_name']
 
             param_card_iterator = main_card
             first_card = param_card_iterator.next(autostart=True)
@@ -402,20 +475,37 @@
             self.mother.check_param_card(pjoin(rw_dir, 'Cards', 'param_card.dat'))
             new_card = open(pjoin(rw_dir, 'Cards', 'param_card.dat')).read()
 
+
         # Find new tag in the banner and add information if needed
         if 'initrwgt' in self.banner:
-            if 'type=\'mg_reweighting\'' in self.banner['initrwgt']:
-                blockpat = re.compile(r'''<weightgroup type=\'mg_reweighting\'\s*>(?P<text>.*?)</weightgroup>''', re.I+re.M+re.S)
+            if 'name=\'mg_reweighting\'' in self.banner['initrwgt']:
+                blockpat = re.compile(r'''<weightgroup name=\'mg_reweighting\'\s*>(?P<text>.*?)</weightgroup>''', re.I+re.M+re.S)
                 before, content, after = blockpat.split(self.banner['initrwgt'])
                 header_rwgt_other = before + after
-                pattern = re.compile('<weight id=\'mg_reweight_(?P<id>\d+)\'>(?P<info>[^<]*)</weight>', re.S+re.I+re.M)
+                pattern = re.compile('<weight id=\'(?:rwgt_(?P<id>\d+)|(?P<id2>[_\w]+))(?P<rwgttype>\s*|_\w+)\'>(?P<info>.*?)</weight>', re.S+re.I+re.M)
                 mg_rwgt_info = pattern.findall(content)
+                
                 maxid = 0
-                for i, diff in mg_rwgt_info:
-                    if int(i) > maxid:
-                        maxid = int(i)
+                for k,(i, fulltag, nlotype, diff) in enumerate(mg_rwgt_info):
+                    if i:
+                        if int(i) > maxid:
+                            maxid = int(i)
+                        mg_rwgt_info[k] = (i, nlotype, diff) # remove the pointless fulltag tag
+                    else:
+                        mg_rwgt_info[k] = (fulltag, nlotype, diff) # remove the pointless id tag
+                         
                 maxid += 1
                 rewgtid = maxid
+                if self.options['rwgt_name']:
+                    #ensure that the entry is not already define if so overwrites it
+                    misc.sprint(mg_rwgt_info)
+                    for (i, nlotype, diff) in mg_rwgt_info[:]:
+                        for flag in type_rwgt:
+                            if 'rwgt_%s' % i == '%s%s' %(self.options['rwgt_name'],flag) or \
+                                i == '%s%s' % (self.options['rwgt_name'], flag):
+                                    logger.warning("tag %s%s already defines, will replace it", self.options['rwgt_name'],flag)
+                                    mg_rwgt_info.remove((i, nlotype, diff))
+                                                
             else:
                 header_rwgt_other = self.banner['initrwgt'] 
                 mg_rwgt_info = []
@@ -425,18 +515,28 @@
             header_rwgt_other = ''
             mg_rwgt_info = []
             rewgtid = 1
-                    
+
         # add the reweighting in the banner information:
         #starts by computing the difference in the cards.
         s_orig = self.banner['slha']
         s_new = new_card
+        
+        #define tag for the run
+        if self.options['rwgt_name']:
+            tag = self.options['rwgt_name']
+        else:
+            tag = str(rewgtid)
+        
         if not self.second_model:
             old_param = check_param_card.ParamCard(s_orig.splitlines())
             new_param =  check_param_card.ParamCard(s_new.splitlines())
             card_diff = old_param.create_diff(new_param)
             if card_diff == '' and not self.second_process:
-                logger.warning(' REWEIGHTING: original card and new card are identical. Bypass this run')
-                return
+                if not __debug__:
+                    logger.warning(' REWEIGHTING: original card and new card are identical. Bypass this run')
+                    return
+                else:
+                    logger.warning(' REWEIGHTING: original card and new card are identical. Run it due to debug mode')
                 #raise self.InvalidCmd, 'original card and new card are identical'
             try:
                 if old_param['sminputs'].get(3)- new_param['sminputs'].get(3) > 1e-3 * new_param['sminputs'].get(3):
@@ -445,35 +545,53 @@
                 logger.debug("error in check of alphas: %s" % str(error))
                 pass #this is a security                
             if not self.second_process:
-                mg_rwgt_info.append((str(rewgtid), card_diff))
+                for name in type_rwgt:
+                    mg_rwgt_info.append((tag, name, card_diff))
             else:
                 str_proc = "\n change process  ".join([""]+self.second_process)
-                mg_rwgt_info.append((str(rewgtid), str_proc + '\n'+ card_diff))
+                for name in type_rwgt:
+                    mg_rwgt_info.append((tag, name, str_proc + '\n'+ card_diff))
         else:
             str_info = "change model %s" % self.second_model
             if self.second_process:
                 str_info += "\n change process  ".join([""]+self.second_process)
             card_diff = str_info
             str_info += '\n' + s_new
-            mg_rwgt_info.append((str(rewgtid), str_info))
-
+            for name in type_rwgt:
+                mg_rwgt_info.append((tag, name, str_info))
         # re-create the banner.
         self.banner['initrwgt'] = header_rwgt_other
-        self.banner['initrwgt'] += '\n<weightgroup type=\'mg_reweighting\'>\n'
-        for tag, diff in mg_rwgt_info:
-            self.banner['initrwgt'] += '<weight id=\'mg_reweight_%s\'>%s</weight>\n' % \
-                                       (tag, diff)
+        self.banner['initrwgt'] += '\n<weightgroup name=\'mg_reweighting\'>\n'
+        for tag, rwgttype, diff in mg_rwgt_info:
+            if tag.isdigit():
+                self.banner['initrwgt'] += '<weight id=\'rwgt_%s%s\'>%s</weight>\n' % \
+                                       (tag, rwgttype, diff)
+            else:
+                self.banner['initrwgt'] += '<weight id=\'%s%s\'>%s</weight>\n' % \
+                                       (tag, rwgttype, diff)
         self.banner['initrwgt'] += '\n</weightgroup>\n'
         self.banner['initrwgt'] = self.banner['initrwgt'].replace('\n\n', '\n')
-            
-        output = open( self.lhe_input.name +'rw', 'w')
-        
+
+
+        start = time.time()
+        cross, ratio, ratio_square,error = {},{},{}, {}
+        for name in type_rwgt + ['orig']:
+            cross[name], error[name] = 0.,0.
+            ratio[name],ratio_square[name] = 0., 0.# to compute the variance and associate error
+
+        if self.output_type == "default":
+            output = open( self.lhe_input.name +'rw', 'w')
+            #write the banner to the output file
+            self.banner.write(output, close_tag=False)
+        else:
+            output = {}
+            for name in type_rwgt:
+                output[name] = open( self.lhe_input.name +'rw'+name, 'w')
+                #write the banner to the output file
+                self.banner.write(output[name], close_tag=False)
         
         logger.info('starts to compute weight for events with the following modification to the param_card:')
-        logger.info(card_diff)
-        
-        #write the banner to the output file
-        self.banner.write(output, close_tag=False)
+        logger.info(card_diff.replace('\n','\nKEEP:'))
         # prepare the output file for the weight plot
         if self.mother:
             out_path = pjoin(self.mother.me_dir, 'Events', 'reweight.lhe')
@@ -495,14 +613,13 @@
             ff.close()
         
         # Loop over all events
-        tag_name = 'mg_reweight_%s' % rewgtid
-        start = time.time()
-        cross = 0
-        ratio, ratio_square = 0, 0 # to compute the variance and associate error
-        
+        if self.options['rwgt_name']:
+            tag_name = self.options['rwgt_name']
+        else:
+            tag_name = 'rwgt_%s' % rewgtid
+                
         os.environ['GFORTRAN_UNBUFFERED_ALL'] = 'y'
         if self.lhe_input.closed:
-            misc.sprint("using", self.lhe_input.name)
             self.lhe_input = lhe_parser.EventFile(self.lhe_input.name)
 
 #        Multicore option not really stable -> not use it
@@ -529,38 +646,67 @@
                 continue
             else:
                 weight = self.calculate_weight(event)
-                cross += weight
-                ratio += weight/event.wgt
-                ratio_square += (weight/event.wgt)**2
+                if not isinstance(weight, dict):
+                    weight = {'':weight}
+                
+                for name in weight:
+                    cross[name] += weight[name]
+                    ratio[name] += weight[name]/event.wgt
+                    ratio_square[name] += (weight[name]/event.wgt)**2
+                    
+                # ensure to have a consistent order of the weights. new one are put 
+                # at the back, remove old position if already defines
+                for tag in type_rwgt:
+                    try:
+                        event.reweight_order.remove('%s%s'  % (tag_name,tag))
+                    except ValueError:
+                        continue
+                
+                event.reweight_order += ['%s%s' % (tag_name,name) for name in type_rwgt]  
                 if self.output_type == "default":
-                    event.reweight_data[tag_name] = weight
-                    #write this event with weight
+                    for name in weight:
+                        if 'orig' in name:
+                            continue             
+                        event.reweight_data['%s%s' % (tag_name,name)] = weight[name]
+                        #write this event with weight
                     output.write(str(event))
                     if self.mother:
-                        event.wgt = weight
+                        event.wgt = weight[type_rwgt[0]]
                         event.reweight_data = {}
                         output2.write(str(event))
-    
                 else:
-                    event.wgt = weight
-                    event.reweight_data = {}
-                    if self.mother:
-                        output2.write(str(event))
-                    else:
-                        output.write(str(event))
+                    for i,name in enumerate(weight):
+                        event.wgt = weight[name]
+                        event.reweight_data = {}
+                        if self.mother and len(weight)==1:
+                            output2.write(str(event))
+                        elif self.mother and i == 0:
+                            output[name].write(str(event))
+                            output2.write(str(event))
+                        else:
+                            output[name].write(str(event))
                 
         # check normalisation of the events:
         if 'event_norm' in self.run_card:
             if self.run_card['event_norm'] == 'average':
-                cross /= event_nb+1
+                for key, value in cross.items():
+                    cross[key] = value / (event_nb+1)
 
                 
         running_time = misc.format_timer(time.time()-start)
         logger.info('All event done  (nb_event: %s) %s' % (event_nb+1, running_time))        
         
-        output.write('</LesHouchesEvents>\n')
-        output.close()
+        
+        if self.output_type == "default":
+            output.write('</LesHouchesEvents>\n')
+            output.close()
+        else:
+            for key in output:
+                output[key].write('</LesHouchesEvents>\n')
+                output.close()
+            
         os.environ['GFORTRAN_UNBUFFERED_ALL'] = 'n'
+        
         if self.mother:
             output2.write('</LesHouchesEvents>\n')
             output2.close()        
@@ -570,14 +716,17 @@
                 results = self.mother.results
                 results.add_run(run_name, self.run_card, current=True)
                 results.add_detail('nb_event', event_nb+1)
-                results.add_detail('cross', cross)
+                name = type_rwgt[0]
+                results.add_detail('cross', cross[name])
                 event_nb +=1
-                variance = ratio_square/event_nb - (ratio/event_nb)**2
-                orig_cross, orig_error = self.orig_cross
-                error = variance/math.sqrt(event_nb) * orig_cross + ratio/event_nb * orig_error
-                
-                results.add_detail('error', error)
-                self.mother.create_plot(mode='reweight', event_path=output2.name,
+                for name in type_rwgt:
+                    variance = ratio_square[name]/event_nb - (ratio[name]/event_nb)**2
+                    orig_cross, orig_error = self.orig_cross
+                    error[name] = variance/math.sqrt(event_nb) * orig_cross + ratio[name]/event_nb * orig_error
+                results.add_detail('error', error[type_rwgt[0]])
+                import madgraph.interface.madevent_interface as ME_interface
+                if isinstance(self.mother, ME_interface.MadEventCmd):
+                    self.mother.create_plot(mode='reweight', event_path=output2.name,
                                         tag=self.run_card['run_tag'])
                 #modify the html output to add the original run
                 if 'plot' in results.current.reweight:
@@ -610,6 +759,7 @@
         
         if self.output_type == "default":
             files.mv(output.name, target)
+            logger.info('Event %s have now the additional weight' % self.lhe_input.name)
         elif self.output_type == "unweight":
             output2.close()
             lhe = lhe_parser.EventFile(output2.name)
@@ -618,25 +768,38 @@
                 results = self.mother.results
                 results.add_detail('nb_event', nb_event)
                 results.current.parton.append('lhe')
-                
+            logger.info('Event %s is now unweighted under the new theory' % output2.name)                
         else:
             files.mv(output2.name, self.lhe_input.name)     
             if self.mother and  hasattr(self.mother, 'results'):
                 results = self.mother.results
-                results.current.parton.append('lhe')   
-
-        logger.info('Event %s have now the additional weight' % self.lhe_input.name)
-        logger.info('new cross-section is : %g pb (indicative error: %g pb)' % (cross,error))
+                results.current.parton.append('lhe')       
+            logger.info('Event %s is now created with new central weight' % output2.name)
+        
+        for name in cross:
+            if name == 'orig':
+                continue
+            logger.info('new cross-section is %s: %g pb (indicative error: %g pb)' %\
+                        ('(%s)' %name if name else '',cross[name], error[name]))
+            
         self.terminate_fortran_executables(new_card_only=True)
         #store result
-        self.all_cross_section[rewgtid] = (cross, error)
+        for name in cross:
+            if name == 'orig':
+                self.all_cross_section[name] = (cross[name], error[name])
+            else:
+                self.all_cross_section[(tag_name,name)] = (cross[name], error[name])
         
+
         # perform the scanning
         if param_card_iterator:
-            for card in param_card_iterator:
+            for i,card in enumerate(param_card_iterator):
+                if self.options['rwgt_name']:
+                    self.options['rwgt_name'] = '%s_%s' % (self.options['rwgt_name'].rsplit('_',1)[0], i+1)
                 card.write(pjoin(rw_dir, 'Cards', 'param_card.dat'))
                 self.exec_cmd("launch --keep_card", printcmd=False, precmd=True)
         
+        self.options['rwgt_name'] = None
     
     def do_set(self, line):
         "Not in help"
@@ -691,20 +854,137 @@
     def calculate_weight(self, event, space=None):
         """space defines where to find the calculator (in multicore)"""
         
+        if self.has_nlo and self.rwgt_mode != "LO":
+            return self.calculate_nlo_weight(event, space)
+        
         if not space: 
             space = self
-        
-        event.parse_reweight()
-        
+        event.parse_reweight()                    
+           
+        # LO reweighting    
         w_orig = self.calculate_matrix_element(event, 0, space)
         w_new =  self.calculate_matrix_element(event, 1, space)
+        
         if w_orig == 0:
             tag, order = event.get_tag_and_order()
             orig_order, Pdir, hel_dict = self.id_to_path[tag]
             misc.sprint(w_orig, w_new)
             misc.sprint(event)
             raise Exception, "Invalid matrix element for original computation (weight=0)"
-        return w_new/w_orig*event.wgt
+
+        return {'orig': event.wgt, '': w_new/w_orig*event.wgt}
+     
+    def calculate_nlo_weight(self, event, space=None):
+
+
+        type_nlo = self.get_weight_names()
+        final_weight = {'orig': event.wgt}
+        
+        if not space: 
+            space = self #for multicore: not use so far
+            
+        event.parse_reweight()
+        event.parse_nlo_weight() 
+
+        #initialise the input to the function which recompute the weight
+        scales2 = []
+        pdg = []
+        bjx = []
+        wgt_tree = [] # reweight for loop-improved type
+        wgt_virt  = [] #reweight b+v together
+        base_wgt = []
+        gs=[]
+        qcdpower = []
+        ref_wgts = [] #for debugging
+        
+        orig_wgt = 0
+        for cevent in event.nloweight.cevents:
+            #check if we need to compute the virtual for that cevent
+            need_V = False # the real is nothing else than the born for a N+1 config
+            all_ctype = [w.type for w in cevent.wgts]
+            if '_nlo' in type_nlo and any(c in all_ctype for c in [2,14,15]):
+                need_V =True
+                
+            w_orig = self.calculate_matrix_element(cevent, 0, space)
+            w_new =  self.calculate_matrix_element(cevent, 1, space)
+            ratio_T = w_new/w_orig
+            if need_V:
+                scale2 = cevent.wgts[0].scales2[0]
+                #for scale2 in set(c.scales2[1] for c in cevent.wgts): 
+                w_origV = self.calculate_matrix_element(cevent, 'V0', space, scale2=scale2)
+                w_newV =  self.calculate_matrix_element(cevent, 'V1', space, scale2=scale2)                    
+                ratio_BV = (w_newV + w_new) / (w_origV + w_orig)
+                ratio_V = w_newV/w_origV
+            else:
+                ratio_V = "should not be used"
+                ratio_BV = "should not be used"
+                
+            for c_wgt in cevent.wgts:
+                orig_wgt += c_wgt.ref_wgt
+                #add the information to the input
+                scales2.append(c_wgt.scales2)
+                pdg.append(c_wgt.pdgs[:2])
+
+                bjx.append(c_wgt.bjks)
+                qcdpower.append(c_wgt.qcdpower)
+                gs.append(c_wgt.gs)
+                ref_wgts.append(c_wgt.ref_wgt)
+                
+                if '_nlo' in type_nlo:
+                    if c_wgt.type in  [2,14,15]:
+                        R = ratio_BV
+                    else:
+                        R = ratio_T
+                    
+                    new_wgt = [c_wgt.pwgt[0] * R,
+                               c_wgt.pwgt[1] * ratio_T,
+                               c_wgt.pwgt[2] * ratio_T]
+                    wgt_virt.append(new_wgt)
+                    
+                if '_tree' in type_nlo:
+                    new_wgt = [c_wgt.pwgt[0] * ratio_T,
+                               c_wgt.pwgt[1] * ratio_T,
+                               c_wgt.pwgt[2] * ratio_T]
+                    wgt_tree.append(new_wgt)
+                base_wgt.append(c_wgt.pwgt[:3])
+        
+        #change the ordering to the fortran one:
+        scales2 = self.invert_momenta(scales2)
+        pdg = self.invert_momenta(pdg)
+        bjx = self.invert_momenta(bjx)
+        # re-compute original weight to reduce numerical inacurracy
+        base_wgt = self.invert_momenta(base_wgt)
+        
+        orig_wgt_check, partial_check = self.combine_wgt(scales2, pdg, bjx, base_wgt, gs, qcdpower, 1., 1.)
+        
+        if '_nlo' in type_nlo:
+            wgt = self.invert_momenta(wgt_virt)
+            with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                new_out, partial = self.combine_wgt(scales2, pdg, bjx, wgt, gs, qcdpower, 1., 1.)
+            # try to correct for precision issue
+            avg = [partial_check[i]/ref_wgts[i] for i in range(len(ref_wgts))]
+            out = sum(partial[i]/avg[i] if 0.85<avg[i]<1.15 else 0 \
+                          for i in range(len(avg)))
+            final_weight['_nlo'] = out/orig_wgt*event.wgt
+
+            
+        if '_tree' in type_nlo:
+            wgt = self.invert_momenta(wgt_tree)
+            with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                out, partial = self.combine_wgt(scales2, pdg, bjx, wgt, gs, qcdpower, 1., 1.)
+            # try to correct for precision issue
+            avg = [partial_check[i]/ref_wgts[i] for i in range(len(ref_wgts))]
+            new_out = sum(partial[i]/avg[i] if 0.85<avg[i]<1.15 else partial[i] \
+                          for i in range(len(avg)))
+            final_weight['_tree'] = new_out/orig_wgt*event.wgt            
+             
+        if '_lo' in type_nlo:
+            w_orig = self.calculate_matrix_element(event, 0, space)
+            w_new =  self.calculate_matrix_element(event, 1, space)            
+            final_weight['_lo'] = w_new/w_orig*event.wgt
+            
+        return final_weight 
+        
      
     @staticmethod   
     def invert_momenta(p):
@@ -726,26 +1006,38 @@
             open(pjoin(Pdir, 'matrix%spy.so' % tag),'w').write(open(pjoin(Pdir, 'matrix2py.so')
                                         ).read().replace('matrix2py', 'matrix%spy' % tag))
     
-    def calculate_matrix_element(self, event, hypp_id, space):
+    def calculate_matrix_element(self, event, hypp_id, space, scale2=0):
         """routine to return the matrix element"""
-
+        
         tag, order = event.get_tag_and_order()
+        if isinstance(hypp_id, str) and hypp_id.startswith('V'):
+            tag = (tag,'V')
+            hypp_id = int(hypp_id[1:])
+            base = "rw_mevirt"
+        else:
+            base = "rw_me"
         
         if (not self.second_model and not self.second_process) or hypp_id==0:
             orig_order, Pdir, hel_dict = self.id_to_path[tag]
         else:
-            orig_order, Pdir, hel_dict = self.id_to_path_second[tag]
+            orig_order, Pdir, hel_dict = self.id_to_path_second[tag] 
+            
 
         run_id = (tag, hypp_id)
 
-
         assert space == self
         start = False
         if run_id in space.calculator:
             external = space.calculator[run_id]
+            #    mod = space.calculator[(run_id,'module')]
+            #    #with misc.chdir(Pdir):
+            #    #    if hypp_id==0:
+            #    #        mod.initialise('param_card_orig.dat')
+            #    #    else:
+            #    #        mod.initialise('param_card.dat')
         elif (not self.second_model and not self.second_process) or hypp_id==0:
             # create the executable for this param_card  
-            subdir = pjoin(self.me_dir,'rw_me ', 'SubProcesses')
+            subdir = pjoin(self.me_dir, base, 'SubProcesses')
             if self.me_dir not in sys.path:
                 sys.path.insert(0,self.me_dir)
             if self.rwgt_dir and self.rwgt_dir not in sys.path:
@@ -760,45 +1052,57 @@
                     if old_module != nb_f2py_module:
                         metag += 1
                         dir_to_f2py_free_mod[(Pdir,0)] = (metag, nb_f2py_module)
-                
                 os.environ['MENUM'] = '2'
                 if not self.rwgt_dir or not os.path.exists(pjoin(Pdir, 'matrix2py.so')):
-                    misc.compile(['matrix2py.so'], cwd=Pdir)
-                    self.rename_f2py_lib(Pdir, 2*metag)
-                    
-                mymod = __import__('rw_me.SubProcesses.%s.matrix%spy' % (Pname, 2*metag), globals(), locals(), [],-1)
+                    misc.compile(['matrix2py.so'], cwd=Pdir)                
+                self.rename_f2py_lib(Pdir, 2*metag)
+                try:
+                    mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, 2*metag), globals(), locals(), [],-1)
+                except:
+                    import platform
+                    if platform.system() == 'Darwin':
+                        os.system('install_name_tool -change libMadLoop.dylib %s/libMadLoop.dylib matrix%spy.so' % (Pdir,2*metag))
+                        mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, 2*metag), globals(), locals(), [],-1)
+                    else:
+                        misc.sprint("fail compilation")
+                        raise
                 S = mymod.SubProcesses
                 P = getattr(S, Pname)
                 mymod = getattr(P, 'matrix%spy' % (2*metag))
                 with misc.chdir(Pdir):
-                    mymod.initialise('param_card_orig.dat')
-                    
+                    with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                        mymod.initialise('param_card_orig.dat')
+                 
+                     
             if hypp_id == 1:
                 #incorrect line
                 metag = dir_to_f2py_free_mod[(Pdir,0)][0]
                 newtag = 2*metag+1
                 self.rename_f2py_lib(Pdir, newtag)
                 try:
-                    mymod = __import__('rw_me.SubProcesses.%s.matrix%spy' % (Pname, newtag), globals(), locals(), [],-1)
+                    mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, newtag), globals(), locals(), [],-1)
                 except Exception, error:
                     os.remove(pjoin(Pdir, 'matrix%spy.so' % newtag))
                     newtag  = "L%s" % newtag
                     os.environ['MENUM'] = newtag
                     misc.compile(['matrix%spy.so' % newtag], cwd=Pdir)
-                    mymod = __import__('rw_me.SubProcesses.%s.matrix%spy' % (Pname, newtag), globals(), locals(), [],-1)
-                
+                    mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, newtag), globals(), locals(), [],-1)
+                 
                 S = mymod.SubProcesses
                 P = getattr(S, Pname)
                 mymod = getattr(P, 'matrix%spy' % newtag)
                 with misc.chdir(Pdir):
-                    mymod.initialise('param_card.dat') 
+                    with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                        mymod.initialise('param_card.dat') 
+                     
             space.calculator[run_id] = mymod.get_me
+            space.calculator[(run_id,'module')] = mymod
             external = space.calculator[run_id]                      
         else:
-            subdir = pjoin(self.me_dir,'rw_me_second', 'SubProcesses')
+            subdir = pjoin(self.me_dir,'%s_second' % base, 'SubProcesses')
             if self.me_dir not in sys.path:
                 sys.path.append(self.me_dir)
-
+ 
             assert hypp_id == 1
             Pname = os.path.basename(Pdir)
             os.environ['MENUM'] = '2'
@@ -813,26 +1117,29 @@
                     metag += 1
                     dir_to_f2py_free_mod[(Pdir,1)] = (metag, nb_f2py_module)
             self.rename_f2py_lib(Pdir, metag)
-            
             try:
-                mymod = __import__("rw_me_second.SubProcesses.%s.matrix%spy" % (Pname, metag))
+                mymod = __import__("%s_second.SubProcesses.%s.matrix%spy" % (base, Pname, metag))
             except ImportError:
                 os.remove(pjoin(Pdir, 'matrix%spy.so' % metag ))
-                metag = "L%s" % metag
+                metag = "L" % metag
                 os.environ['MENUM'] = str(metag)
                 misc.compile(['matrix%spy.so' % metag], cwd=pjoin(subdir, Pdir))
-                mymod = __import__("rw_me_second.SubProcesses.%s.matrix%spy" % (Pname, metag))
+                mymod = __import__("%s_second.SubProcesses.%s.matrix%spy" % (base, Pname, metag))
+                     
             reload(mymod)
             S = mymod.SubProcesses
             P = getattr(S, Pname)
-            mymod = getattr(P, 'matrix%spy' % metag)                
-            with misc.chdir(Pdir):      
-                mymod.initialise('param_card.dat')              
+            mymod = getattr(P, 'matrix%spy' % metag)
+            with misc.chdir(Pdir):   
+                with misc.stdchannel_redirected(sys.stdout, os.devnull):   
+                    mymod.initialise('param_card.dat')              
             space.calculator[run_id] = mymod.get_me
+            space.calculator[(run_id,'module')] = mymod
             external = space.calculator[run_id]                
         
-        
         p = self.invert_momenta(event.get_momenta(orig_order))
+        
+
         # add helicity information
         
         hel_order = event.get_helicity(orig_order)
@@ -843,7 +1150,15 @@
 
         with misc.chdir(Pdir):
             with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                me_value = external(p,event.aqcd, nhel)
+                if 'V' in tag or \
+                   (hypp_id ==1  and self.second_process and any('sqrvirt' in l for l in self.second_process)):
+                    me_value = external(p,event.aqcd, math.sqrt(scale2), nhel)
+                else:
+                    try:
+                        me_value = external(p,event.aqcd, nhel)
+                    except TypeError:
+                        me_value = external(p,event.aqcd, math.sqrt(scale2), nhel)                    
+
         # for NLO we have also the stability status code
         if isinstance(me_value, tuple):
             me_value, code = me_value
@@ -875,17 +1190,28 @@
                 split = line.split()
                 if len(split) == 4:
                     cross, error = float(split[0]), float(split[1])
+        if 'orig' not in self.all_cross_section:
             logger.info('Original cross-section: %s +- %s pb' % (cross, error))
-        
+        else: 
+            logger.info('Original cross-section: %s +- %s pb (cross-section from sum of weights: %s)' % (cross, error, self.all_cross_section['orig'][0]))
         logger.info('Computed cross-section:')
         keys = self.all_cross_section.keys()
         keys.sort()
         for key in keys:
-            logger.info('%s : %s +- %s pb' % (key,self.all_cross_section[key][0],self.all_cross_section[key][1] ))  
+            if key == 'orig':
+                continue
+            logger.info('%s : %s +- %s pb' % (key[0] if not key[1] else '%s%s' % key,
+                self.all_cross_section[key][0],self.all_cross_section[key][1] ))  
         self.terminate_fortran_executables()
     
         if self.rwgt_dir:
             self.save_to_pickle()
+        
+        with misc.stdchannel_redirected(sys.stdout, os.devnull):
+            for run_id in self.calculator:
+                del self.calculator[run_id]
+            del self.calculator
+        
             
     def __del__(self):
         self.do_quit('')
@@ -895,165 +1221,89 @@
         """Adding one element to the list based on the matrix element"""
         
 
-
-
-        
-    
     @misc.mute_logger()
     def create_standalone_directory(self, second=False):
         """generate the various directory for the weight evaluation"""
         
+        data={}
+        if not second:
+            data['paths'] = ['rw_me', 'rw_mevirt']
+            # model
+            info = self.banner.get('proc_card', 'full_model_line')
+            if '-modelname' in info:
+                data['mg_names'] = False
+            else:
+                data['mg_names'] = True
+            data['model_name'] = self.banner.get('proc_card', 'model')
+            #processes
+            data['processes'] = [line[9:].strip() for line in self.banner.proc_card
+                     if line.startswith('generate')]
+            data['processes'] += [' '.join(line.split()[2:]) for line in self.banner.proc_card
+                      if re.search('^\s*add\s+process', line)]  
+            #object_collector
+            self.id_to_path = {}
+            data['id2path'] = self.id_to_path
+        else:
+            data['paths'] = ['rw_me_second', 'rw_mevirt_second']
+            # model
+            if self.second_model:
+                data['mg_names'] = True
+                if ' ' in self.second_model:
+                    args = self.second_model.split()
+                    if '--modelname' in args:
+                        data['mg_names'] = False
+                    data['model_name'] = args[0]
+                else:
+                    data['model_name'] = self.second_model
+            else:
+                data['model_name'] = None
+            #processes
+            if self.second_process:
+                data['processes'] = self.second_process
+            else:
+                data['processes'] = [line[9:].strip() for line in self.banner.proc_card
+                                 if line.startswith('generate')]
+                data['processes'] += [' '.join(line.split()[2:]) 
+                                      for line in self.banner.proc_card
+                                      if re.search('^\s*add\s+process', line)]
+            #object_collector
+            self.id_to_path_second = {}   
+            data['id2path'] = self.id_to_path_second 
+        
         # 0. clean previous run ------------------------------------------------
         if not self.rwgt_dir:
             path_me = self.me_dir
         else:
             path_me = self.rwgt_dir
         try:
-            shutil.rmtree(pjoin(path_me,'rw_me'))
+            shutil.rmtree(pjoin(path_me,data['paths'][0]))
+        except Exception: 
+            pass
+        try:
+            shutil.rmtree(pjoin(path_me, data['paths'][1]))
         except Exception: 
             pass
 
-        # 1. Load model---------------------------------------------------------  
+        # 1. prepare the interface----------------------------------------------
+        mgcmd = self.mg5cmd
         complex_mass = False   
         has_cms = re.compile(r'''set\s+complex_mass_scheme\s*(True|T|1|true|$|;)''')
         for line in self.banner.proc_card:
             if line.startswith('set'):
-                self.mg5cmd.exec_cmd(line, printcmd=False, precmd=False, postcmd=False)
+                mgcmd.exec_cmd(line, printcmd=False, precmd=False, postcmd=False)
                 if has_cms.search(line):
                     complex_mass = True
             elif line.startswith('define'):
                 try:
-                    self.mg5cmd.exec_cmd(line, printcmd=False, precmd=False, postcmd=False)
+                    mgcmd.exec_cmd(line, printcmd=False, precmd=False, postcmd=False)
                 except Exception:
                     pass 
                           
-        info = self.banner.get('proc_card', 'full_model_line')
-        if '-modelname' in info:
-            mg_names = False
-        else:
-            mg_names = True
-        model_name = self.banner.get('proc_card', 'model')
-        if model_name:
-            self.load_model(model_name, mg_names, complex_mass)
-        else:
+        # 1. Load model---------------------------------------------------------  
+        if  not data['model_name'] and not second:
             raise self.InvalidCmd('Only UFO model can be loaded in this module.')
-        
-        mgcmd = self.mg5cmd
-        modelpath = self.model.get('modelpath')
-        if os.path.basename(modelpath) != mgcmd._curr_model['name']:
-            name, restrict = mgcmd._curr_model['name'].rsplit('-',1)
-            if os.path.exists(pjoin(os.path.dirname(modelpath),name, 'restrict_%s.dat' % restrict)):
-                modelpath = pjoin(os.path.dirname(modelpath), mgcmd._curr_model['name'])
-            
-        commandline="import model %s " % modelpath
-        mgcmd.exec_cmd(commandline)
-        
-        #multiparticles
-        for name, content in self.banner.get('proc_card', 'multiparticles'):
-            mgcmd.exec_cmd("define %s = %s" % (name, content))
-        
-        # 2. compute the production matrix element -----------------------------
-        processes = [line[9:].strip() for line in self.banner.proc_card
-                     if line.startswith('generate')]
-        processes += [' '.join(line.split()[2:]) for line in self.banner.proc_card
-                      if re.search('^\s*add\s+process', line)]   
-        mgcmd.exec_cmd("set group_subprocesses False")
-
-        logger.info('generating the square matrix element for reweighting')
-        start = time.time()
-        commandline=''
-        for proc in processes:
-            if '[' not in proc:
-                commandline += "add process %s ;" % proc
-            else:
-                commandline += self.get_LO_definition_from_NLO(proc)
-        
-        commandline = commandline.replace('add process', 'generate',1)
-        logger.info(commandline)
-        try:
-            mgcmd.exec_cmd(commandline, precmd=True, errorhandling=False)
-        except diagram_generation.NoDiagramException:
-            commandline=''
-            for proc in processes:
-                if '[' not in proc:
-                    raise
-                commandline += "add process %s ;" % proc
-            commandline = commandline.replace('add process', 'generate',1)
-            logger.info("RETRY with %s", commandline)
-            mgcmd.exec_cmd(commandline, precmd=True)
-        except Exception, error:
-            raise
-        
-        commandline = 'output standalone_rw %s' % pjoin(path_me,'rw_me')
-        mgcmd.exec_cmd(commandline, precmd=True)        
-        logger.info('Done %.4g' % (time.time()-start))
-        self.has_standalone_dir = True
-        
-
-        # 3. Store id to directory information ---------------------------------
-        matrix_elements = mgcmd._curr_matrix_elements.get_matrix_elements()
-        
-        self.id_to_path = {}
-        for me in matrix_elements:
-            for proc in me.get('processes'):
-                initial = []    #filled in the next line
-                final = [l.get('id') for l in proc.get('legs')\
-                      if l.get('state') or initial.append(l.get('id'))]
-                order = (initial, final)
-                tag = proc.get_initial_final_ids()
-                decay_finals = proc.get_final_ids_after_decay()
-
-                if tag[1] != decay_finals:
-                    order = (initial, list(decay_finals))
-                    decay_finals.sort()
-                    tag = (tag[0], tuple(decay_finals))
-                Pdir = pjoin(path_me, 'rw_me', 'SubProcesses', 
-                                  'P%s' % me.get('processes')[0].shell_string())
-                assert os.path.exists(Pdir), "Pdir %s do not exists" % Pdir                        
-                if tag in self.id_to_path:
-                    if not Pdir == self.id_to_path[tag][1]:
-                        misc.sprint(tag, Pdir, self.id_to_path[tag][1])
-                        raise self.InvalidCmd, '2 different process have the same final states. This module can not handle such situation'
-                    else:
-                        continue
-                # build the helicity dictionary
-                hel_nb = 0
-                hel_dict = {9:0} # unknown helicity -> use full ME
-                for helicities in me.get_helicity_matrix():
-                    hel_nb +=1 #fortran starts at 1
-                    hel_dict[tuple(helicities)] = hel_nb
-
-                self.id_to_path[tag] = [order, Pdir, hel_dict]        
-        
-        # 3. If we need a new model/process-------------------------------------
-        if self.second_model or self.second_process:
-            self.create_second_standalone_directory()    
-
-    @misc.mute_logger()
-    def create_second_standalone_directory(self, second=False):
-        """generate the various directory for the weight evaluation"""
-        
-        # 0. clean previous run ------------------------------------------------
-        path_me = self.me_dir
-        try:
-            shutil.rmtree(pjoin(path_me,'rw_me_second'))
-        except Exception: 
-            pass
-        
-        mgcmd = self.mg5cmd
-        # 1. Load model---------------------------------------------------------  
-        if self.second_model:
-            use_mgdefault= True
-            complex_mass = False
-            if ' ' in self.second_model:
-                args = self.second_model.split()
-                if '--modelname' in args:
-                    use_mgdefault = False
-                model_name = args[0]
-            else:
-                model_name = self.second_model
-            self.load_model(model_name, use_mgdefault, complex_mass)
-                
+        elif data['model_name']:
+            self.load_model(data['model_name'], data['mg_names'], complex_mass)
             modelpath = self.model.get('modelpath')
             if os.path.basename(modelpath) != mgcmd._curr_model['name']:
                 name, restrict = mgcmd._curr_model['name'].rsplit('-',1)
@@ -1061,46 +1311,58 @@
                     modelpath = pjoin(os.path.dirname(modelpath), mgcmd._curr_model['name'])
                 
             commandline="import model %s " % modelpath
-            if not use_mgdefault:
+            if not data['mg_names']:
                 commandline += ' -modelname '
             mgcmd.exec_cmd(commandline)
             
+            #multiparticles
+            for name, content in self.banner.get('proc_card', 'multiparticles'):
+                mgcmd.exec_cmd("define %s = %s" % (name, content))
+        
         # 2. compute the production matrix element -----------------------------
-        if self.second_process:
-            processes = self.second_process
+        has_nlo = False  
+        mgcmd.exec_cmd("set group_subprocesses False")
+
+        if not second:
+            logger.info('generating the square matrix element for reweighting')
         else:
-            processes = [line[9:].strip() for line in self.banner.proc_card
-                         if line.startswith('generate')]
-            processes += [' '.join(line.split()[2:]) for line in self.banner.proc_card
-                          if re.search('^\s*add\s+process', line)]   
-        mgcmd.exec_cmd("set group_subprocesses False")
-
-        logger.info('generating the square matrix element for reweighting')
+            logger.info('generating the square matrix element for reweighting (second model and/or processes)')
         start = time.time()
         commandline=''
-        for proc in processes:
+        for proc in data['processes']:
             if '[' not in proc:
-                commandline+="add process %s ;" % proc
-            elif 'sqrvirt' in proc:
-                commandline+="add process %s ;" % proc
+                commandline += "add process %s ;" % proc
             else:
-                raise self.InvalidCmd('NLO processes can\'t be reweight (for Loop induced reweighting use [sqrvirt =])')
+                has_nlo = True
+                commandline += self.get_LO_definition_from_NLO(proc)
         
         commandline = commandline.replace('add process', 'generate',1)
         logger.info(commandline)
-        mgcmd.exec_cmd(commandline, precmd=True)
-        
-        matrix_elements = mgcmd._curr_matrix_elements.get_matrix_elements()
-        
-        commandline = 'output standalone_rw %s' % pjoin(path_me,'rw_me_second')
+        try:
+            mgcmd.exec_cmd(commandline, precmd=True, errorhandling=False)
+        except diagram_generation.NoDiagramException:
+            commandline=''
+            for proc in data['processes']:
+                if '[' not in proc:
+                    raise
+                commandline += "add process %s ;" % proc
+            commandline = commandline.replace('add process', 'generate',1)
+            logger.info("RETRY with %s", commandline)
+            mgcmd.exec_cmd(commandline, precmd=True)
+            has_nlo = False
+        except Exception, error:
+            raise
+        
+        commandline = 'output standalone_rw %s' % pjoin(path_me,data['paths'][0])
         mgcmd.exec_cmd(commandline, precmd=True)        
         logger.info('Done %.4g' % (time.time()-start))
+        self.has_standalone_dir = True
         
+
         # 3. Store id to directory information ---------------------------------
         matrix_elements = mgcmd._curr_matrix_elements.get_matrix_elements()
         
-        self.id_to_path_second = {}
-        to_check = [] # list of tag that do not have a Pdir at creation time. 
+        to_check = [] # list of tag that do not have a Pdir at creation time.
         for me in matrix_elements:
             for proc in me.get('processes'):
                 initial = []    #filled in the next line
@@ -1114,43 +1376,196 @@
                     order = (initial, list(decay_finals))
                     decay_finals.sort()
                     tag = (tag[0], tuple(decay_finals))
-                Pdir = pjoin(path_me, 'rw_me_second', 'SubProcesses', 
+                Pdir = pjoin(path_me, data['paths'][0], 'SubProcesses', 
                                   'P%s' % me.get('processes')[0].shell_string())
+
                 if not os.path.exists(Pdir):
                     to_check.append(tag)
-                    continue
-                if tag in self.id_to_path_second:
-                    if not Pdir == self.id_to_path_second[tag][1]:
-                        misc.sprint(tag, Pdir, self.id_to_path_second[tag][1])
+                    continue                        
+                if tag in data['id2path']:
+                    if not Pdir == data['id2path'][tag][1]:
+                        misc.sprint(tag, Pdir, data['id2path'][tag][1])
                         raise self.InvalidCmd, '2 different process have the same final states. This module can not handle such situation'
                     else:
                         continue
-                
                 # build the helicity dictionary
                 hel_nb = 0
                 hel_dict = {9:0} # unknown helicity -> use full ME
                 for helicities in me.get_helicity_matrix():
                     hel_nb +=1 #fortran starts at 1
-                    hel_dict[tuple(helicities)] = hel_nb    
-                self.id_to_path_second[tag] = [order, Pdir, hel_dict]
+                    hel_dict[tuple(helicities)] = hel_nb
 
+                data['id2path'][tag] = [order, Pdir, hel_dict]        
+ 
         for tag in to_check:
             if tag not in self.id_to_path:
                 logger.warning("no valid path for %s" % (tag,))
                 #raise self.InvalidCmd, "no valid path for %s" % (tag,)
         
-        # 4. Check MadLoopParam
-        if os.path.exists(pjoin(path_me, 'rw_me_second', 'Cards', 'MadLoopParams.dat')):
-            MLCard = banner.MadLoopParam(pjoin(path_me, 'rw_me_second', 'Cards', 'MadLoopParams.dat'))
+        # 4. Check MadLoopParam for Loop induced
+        if os.path.exists(pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat')):
+            MLCard = banner.MadLoopParam(pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat'))
             MLCard.set('WriteOutFilters', False)
             MLCard.set('UseLoopFilter', False)
             MLCard.set("DoubleCheckHelicityFilter", False)
             MLCard.set("HelicityFilterLevel", 0)
-            MLCard.write(pjoin(path_me, 'rw_me_second', 'SubProcesses', 'MadLoopParams.dat'),
-                         pjoin(path_me, 'rw_me_second', 'Cards', 'MadLoopParams.dat'), 
+            MLCard.write(pjoin(path_me, data['paths'][0], 'SubProcesses', 'MadLoopParams.dat'),
+                         pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat'), 
                          commentdefault=False)
             
-            
+        # 5. create the virtual for NLO reweighting  ---------------------------
+        if has_nlo and 'NLO' in self.rwgt_mode:
+            # Do not pass here for LO/NLO_tree
+            start = time.time()
+            commandline=''
+            for proc in data['processes']:
+                if '[' not in proc:
+                    pass
+                else:
+                    proc = proc.replace('[', '[ virt=')
+                    commandline += "add process %s ;" % proc
+            # deactivate golem since it creates troubles
+            old_options = dict(mgcmd.options)
+            if mgcmd.options['golem'] or mgcmd.options['pjfry']:
+                logger.info(" When doing NLO reweighting, MG5aMC cannot use the loop reduction algorithms Golem and/or PJFry++")
+            mgcmd.options['golem'] = None            
+            mgcmd.options['pjfry'] = None 
+            commandline = commandline.replace('add process', 'generate',1)
+            logger.info(commandline)
+            mgcmd.exec_cmd(commandline, precmd=True)
+            commandline = 'output standalone_rw %s -f' % pjoin(path_me, data['paths'][1])
+            mgcmd.exec_cmd(commandline, precmd=True) 
+            
+            #put back golem to original value
+            mgcmd.options['golem'] = old_options['golem']
+            mgcmd.options['pjfry'] = old_options['pjfry']
+            # update make_opts
+            m_opts = {}
+            if mgcmd.options['lhapdf']:
+                #lhapdfversion = subprocess.Popen([mgcmd.options['lhapdf'], '--version'], 
+                #        stdout = subprocess.PIPE).stdout.read().strip()[0]
+                m_opts['lhapdf'] = True
+                m_opts['f2pymode'] = True
+                m_opts['lhapdfversion'] = 5 # 6 always fail on my computer since 5 is compatible but slower always use 5
+                m_opts['llhapdf'] = self.mother.get_lhapdf_libdir()                       
+            else:
+                raise Exception, "NLO reweighting requires LHAPDF to work correctly"
+ 
+            path = pjoin(path_me,data['paths'][1], 'Source', 'make_opts')             
+            common_run_interface.CommonRunCmd.update_make_opts_full(path, m_opts)      
+            logger.info('Done %.4g' % (time.time()-start))
+
+
+            # Download LHAPDF SET
+            common_run_interface.CommonRunCmd.install_lhapdf_pdfset_static(\
+                mgcmd.options['lhapdf'], None, self.banner.run_card.get_lhapdf_id())
+            
+            # now store the id information             
+            matrix_elements = mgcmd._curr_matrix_elements.get_matrix_elements()            
+            for me in matrix_elements:
+                for proc in me.get('processes'):
+                    initial = []    #filled in the next line
+                    final = [l.get('id') for l in proc.get('legs')\
+                          if l.get('state') or initial.append(l.get('id'))]
+                    order = (initial, final)
+                    tag = proc.get_initial_final_ids()
+                    decay_finals = proc.get_final_ids_after_decay()
+    
+                    if tag[1] != decay_finals:
+                        order = (initial, list(decay_finals))
+                        decay_finals.sort()
+                        tag = (tag[0], tuple(decay_finals))
+                    Pdir = pjoin(path_me, data['paths'][1], 'SubProcesses', 
+                                      'P%s' % me.get('processes')[0].shell_string())
+                    assert os.path.exists(Pdir), "Pdir %s do not exists" % Pdir                        
+                    if (tag,'V') in data['id2path']:
+                        if not Pdir == data['id2path'][(tag,'V')][1]:
+                            misc.sprint(tag, Pdir, self.id_to_path[(tag,'V')][1])
+                            raise self.InvalidCmd, '2 different process have the same final states. This module can not handle such situation'
+                        else:
+                            continue
+                    # build the helicity dictionary
+                    hel_nb = 0
+                    hel_dict = {9:0} # unknown helicity -> use full ME
+                    for helicities in me.get_helicity_matrix():
+                        hel_nb +=1 #fortran starts at 1
+                        hel_dict[tuple(helicities)] = hel_nb
+    
+                    data['id2path'][(tag,'V')] = [order, Pdir, hel_dict] 
+                
+            #compile the module to combine the weight
+            misc.compile(cwd=pjoin(path_me, data['paths'][1], 'Source'))
+            #link it 
+            with misc.chdir(pjoin(path_me)):
+                if path_me not in sys.path:
+                    sys.path.insert(0, path_me)
+                mymod = __import__('%s.Source.rwgt2py' % data['paths'][1], globals(), locals(), [],-1)
+                mymod =  mymod.Source.rwgt2py
+                with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                    mymod.initialise([self.banner.run_card['lpp1'], 
+                                  self.banner.run_card['lpp2']],
+                                 self.banner.run_card.get_lhapdf_id())
+                self.combine_wgt = mymod.get_wgt
+                
+        elif has_nlo and not second and self.rwgt_mode == ['NLO_tree']:
+            # We do not have any virtual reweighting to do but we still have to
+            #combine the weights.
+            #Idea:create a fake directory.
+            start = time.time()
+            commandline='import model loop_sm;generate g g > e+ ve [virt=QCD]'
+            # deactivate golem since it creates troubles
+            old_options = dict(mgcmd.options)
+            mgcmd.options['golem'] = None            
+            mgcmd.options['pjfry'] = None 
+            commandline = commandline.replace('add process', 'generate',1)
+            logger.info(commandline)
+            mgcmd.exec_cmd(commandline, precmd=True)
+            commandline = 'output standalone_rw %s -f' % pjoin(path_me, data['paths'][1])
+            mgcmd.exec_cmd(commandline, precmd=True)    
+            #put back golem to original value
+            mgcmd.options['golem'] = old_options['golem']
+            mgcmd.options['pjfry'] = old_options['pjfry']
+            # update make_opts
+            m_opts = {}
+            if mgcmd.options['lhapdf']:
+                #lhapdfversion = subprocess.Popen([mgcmd.options['lhapdf'], '--version'], 
+                #        stdout = subprocess.PIPE).stdout.read().strip()[0]
+                m_opts['lhapdf'] = True
+                m_opts['f2pymode'] = True
+                m_opts['lhapdfversion'] = 5 # 6 always fail on my computer since 5 is compatible but slower always use 5
+                m_opts['llhapdf'] = self.mother.get_lhapdf_libdir()                        
+            else:
+                raise Exception, "NLO_tree reweighting requires LHAPDF to work correctly"
+ 
+            path = pjoin(path_me,data['paths'][1], 'Source', 'make_opts')             
+            common_run_interface.CommonRunCmd.update_make_opts_full(path, m_opts)      
+            logger.info('Done %.4g' % (time.time()-start))
+
+            # Download LHAPDF SET
+            common_run_interface.CommonRunCmd.install_lhapdf_pdfset_static(\
+                mgcmd.options['lhapdf'], None, self.banner.run_card.get_lhapdf_id())
+            
+            #compile the module to combine the weight
+            misc.compile(cwd=pjoin(path_me, data['paths'][1], 'Source'))
+            #link it 
+            with misc.chdir(pjoin(path_me)):
+                if path_me not in sys.path:
+                    sys.path.insert(0, path_me)
+                mymod = __import__('%s.Source.rwgt2py' % data['paths'][1], globals(), locals(), [],-1)
+                mymod =  mymod.Source.rwgt2py
+                with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                    mymod.initialise([self.banner.run_card['lpp1'], 
+                                  self.banner.run_card['lpp2']],
+                                 self.banner.run_card.get_lhapdf_id())
+                self.combine_wgt = mymod.get_wgt
+            
+             
+        # 6. If we need a new model/process-------------------------------------
+        if (self.second_model or self.second_process) and not second:
+            self.create_standalone_directory(second=True)    
+
+        if not second:
+            self.has_nlo = has_nlo
             
         
         
@@ -1192,6 +1607,8 @@
         else:
             to_save['second_model'] = None
         to_save['rwgt_dir'] = self.rwgt_dir
+        to_save['has_nlo'] = self.has_nlo
+        to_save['rwgt_mode'] = self.rwgt_mode
 
         name = pjoin(self.rwgt_dir, 'rw_me', 'rwgt.pkl')
         save_load_object.save_to_file(name, to_save)
@@ -1202,7 +1619,8 @@
         obj = save_load_object.load_from_file( pjoin(self.rwgt_dir, 'rw_me', 'rwgt.pkl'))
         
         self.has_standalone_dir = True
-        self.options = {'curr_dir': os.path.realpath(os.getcwd())}
+        self.options = {'curr_dir': os.path.realpath(os.getcwd()),
+                        'rewgt_name': None}
         
         old_rwgt = obj['rwgt_dir']
            
@@ -1222,8 +1640,20 @@
         self.processes = obj['processes']
         self.second_process = obj['second_process']
         self.second_model = obj['second_model']
-
-        
+        self.has_nlo = obj['has_nlo']
+        if not self.rwgt_mode:
+            self.rwgt_mode = obj['rwgt_mode']
+            logger.info("mode set to %s" % self.rwgt_mode)
+        if self.has_nlo:
+            path = pjoin(obj['rwgt_dir'], 'rw_mevirt','Source')
+            sys.path.insert(0, path)
+            mymod = __import__('rwgt2py', globals(), locals())
+            with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                mymod.initialise([self.banner.run_card['lpp1'], 
+                              self.banner.run_card['lpp2']],
+                             self.banner.run_card.get_lhapdf_id())
+            self.combine_wgt = mymod.get_wgt
+                    
         
         
         

=== modified file 'madgraph/iolibs/export_fks.py'
--- madgraph/iolibs/export_fks.py	2016-02-24 13:54:17 +0000
+++ madgraph/iolibs/export_fks.py	2016-03-16 17:59:03 +0000
@@ -1588,6 +1588,11 @@
         for file in linkfiles:
             ln('../../../lib/%s' % file)
 
+        linkfiles = ['coef_specs.inc']
+
+        for file in linkfiles:        
+            ln('../../../Source/DHELAS/%s' % file)
+
         # Return to original PWD
         os.chdir(cwd)
 
@@ -2632,9 +2637,7 @@
             lines.append('INTEGER IDEN_VALUES(%d)' % len(info_list))
             lines.append('DATA IDEN_VALUES /' + \
                          ', '.join(['%d' % ( 
-                         fks_born.born_matrix_element.get_denominator_factor() / \
-                         fks_born.born_matrix_element['identical_particle_factor'] * \
-                         fks_born.real_processes[info['n_me'] - 1].matrix_element['identical_particle_factor'] ) \
+                         fks_born.born_matrix_element.get_denominator_factor() ) \
                          for info in info_list]) + '/')
         else:
             # otherwise use the born
@@ -3286,6 +3289,11 @@
         for file in linkfiles:
             ln('../../../lib/%s' % file)
 
+        linkfiles = ['coef_specs.inc']
+
+        for file in linkfiles:        
+            ln('../../../Source/DHELAS/%s' % file)
+
         # Return to original PWD
         os.chdir(cwd)
 
@@ -3303,12 +3311,10 @@
         filename = os.path.join(self.dir_path, 'Source', 'DHELAS', 'coef_specs.inc')
 
         replace_dict = {}
-        replace_dict['max_lwf_size'] = 4 
-
+        replace_dict['max_lwf_size'] = 4
         replace_dict['vertex_max_coefs'] = max(\
                 [q_polynomial.get_number_of_coefs_for_rank(n) 
                     for n in max_loop_vertex_ranks])
-
         IncWriter=writers.FortranWriter(filename,'w')
         IncWriter.writelines("""INTEGER MAXLWFSIZE
                            PARAMETER (MAXLWFSIZE=%(max_lwf_size)d)

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2016-03-04 16:16:23 +0000
+++ madgraph/iolibs/export_v4.py	2016-03-16 17:59:03 +0000
@@ -42,6 +42,7 @@
 import madgraph.iolibs.template_files as template_files
 import madgraph.iolibs.ufo_expression_parsers as parsers
 import madgraph.iolibs.helas_call_writers as helas_call_writers
+import madgraph.interface.common_run_interface as common_run_interface
 import madgraph.various.diagram_symmetry as diagram_symmetry
 import madgraph.various.misc as misc
 import madgraph.various.banner as banner_mod
@@ -1687,35 +1688,16 @@
         f2py_compiler = compilers['f2py']
         if not f2py_compiler:
             f2py_compiler = 'f2py'
-
-            
+        for_update= {'DEFAULT_F_COMPILER':compiler,
+                     'DEFAULT_F2PY_COMPILER':f2py_compiler}
         make_opts = pjoin(root_dir, 'Source', 'make_opts')
-        lines = open(make_opts).read().split('\n')
-        FC_re = re.compile('^(\s*)(FC|F2PY)\s*=\s*(.+)\s*$')
-        for iline, line in enumerate(lines):
 
-            FC_result = FC_re.match(line)
-            if FC_result:
-                if 'FC' == FC_result.group(2):
-                    if compiler != FC_result.group(3):
-                        mod = True
-                    lines[iline] = FC_result.group(1) + "FC=" + compiler
-                elif 'F2PY' ==  FC_result.group(2):
-                    if f2py_compiler != FC_result.group(3):
-                        mod = True
-                    lines[iline] = FC_result.group(1) + "F2PY=" + f2py_compiler
-                                    
-        if not mod:
-            return
-        
         try:
-            outfile = open(make_opts, 'w')
+            common_run_interface.CommonRunCmd.update_make_opts_full(
+                            make_opts, for_update)
         except IOError:
             if root_dir == self.dir_path:
-                logger.info('Fail to set compiler. Trying to continue anyway.')
-            return
-        outfile.write('\n'.join(lines))
-
+                logger.info('Fail to set compiler. Trying to continue anyway.')            
 
     def replace_make_opt_c_compiler(self, compiler, root_dir = ""):
         """Set CXX=compiler in Source/make_opts.
@@ -1738,56 +1720,26 @@
                     is_lc = True
                 else:
                     is_lc = False
-    
-        mod = False #avoid to rewrite the file if not needed
+        
+        # list of the variable to set in the make_opts file
+        for_update= {'DEFAULT_CPP_COMPILER':compiler,
+                     'MACFLAG':'-mmacosx-version-min=10.7' if is_clang and is_lc else '',
+                     'STDLIB': '-lc++' if is_lc else '-lstdc++',
+                     'STDLIB_FLAG': '-stdlib=libc++' if is_lc and is_clang else ''
+                     }
+        
         if not root_dir:
             root_dir = self.dir_path
         make_opts = pjoin(root_dir, 'Source', 'make_opts')
-        lines = open(make_opts).read().split('\n')
-        CC_re = re.compile('^(\s*)CXX\s*=\s*(.+)\s*$')
-        for iline, line in enumerate(lines):
-            CC_result = CC_re.match(line)
-            if CC_result:
-                if compiler != CC_result.group(2):
-                    mod = True
-                lines[iline] = CC_result.group(1) + "CXX=" + compiler
-            if 'LDFLAGS=-lc++' in line:
-                if not is_lc:
-                    lines[iline] = 'LDFLAGS=-lstdc++'
-                    mod = True
-                elif is_clang and not 'mmacosx-version' in line:
-                    lines[iline] += " -mmacosx-version-min=10.7"
-                    mod=True
-                elif not is_clang and 'mmacosx-version' in line:
-                    lines[iline] = lines[iline].replace('-mmacosx-version-min=10.7', '')
-                    mod = True
-            elif 'LDFLAGS=-lstdc++' in line:
-                if is_lc and is_clang:
-                    lines[iline] = 'LDFLAGS=-lc++  -mmacosx-version-min=10.7'
-                    mod = True
 
-        if is_clang and is_lc:
-            CFLAGS_re=re.compile('^(\s*)CFLAGS\s*=\s*(.+)\s*$')
-            CXXFLAGS_re=re.compile('^(\s*)CXXFLAGS\s*=\s*(.+)\s*$')
-            flags= '-O -stdlib=libc++ -mmacosx-version-min=10.7'
-            for iline, line in enumerate(lines):
-                CF_result = CFLAGS_re.match(line)
-                CXXF_result = CXXFLAGS_re.match(line)
-                if CF_result and CF_result.group(2) != flags:
-                    lines[iline] = CF_result.group(1) + "CFLAGS= " + flags
-                    mod=True
-                if CXXF_result and CXXF_result.group(1)!= flags:
-                    lines[iline] = CXXF_result.group(1) + "CXXFLAGS= " + flags
-                    mod=True                    
-        if not mod:
-            return
         try:
-            outfile = open(make_opts, 'w')
+            common_run_interface.CommonRunCmd.update_make_opts_full(
+                            make_opts, for_update)
         except IOError:
             if root_dir == self.dir_path:
-                logger.info('Fail to set compiler. Trying to continue anyway.')
-            return
-        outfile.write('\n'.join(lines))
+                logger.info('Fail to set compiler. Trying to continue anyway.')  
+    
+        return
 
 #===============================================================================
 # ProcessExporterFortranSA
@@ -1918,6 +1870,15 @@
         open(pjoin(self.dir_path,'__init__.py'),'w')
         open(pjoin(self.dir_path,'SubProcesses','__init__.py'),'w')
 
+        if 'mode' in self.opt and self.opt['mode'] == "reweight":
+            #add the module to hande the NLO weight
+            files.copytree(pjoin(MG5DIR, 'Template', 'RWGTNLO'),
+                          pjoin(self.dir_path, 'Source'))
+            files.copytree(pjoin(MG5DIR, 'Template', 'NLO', 'Source', 'PDF'),
+                           pjoin(self.dir_path, 'Source', 'PDF'))
+            self.write_pdf_opendata()
+
+
 
     def compiler_choice(self, compiler):
         """ Different daughter classes might want different compilers.
@@ -5751,7 +5712,7 @@
         fsock.writelines("""include \'input.inc\'
                             include \'coupl.inc\'""")
         fsock.writelines("""
-                            MU_R = mu_r2
+                            if (mu_r2.gt.0d0) MU_R = mu_r2
                             G = SQRT(4.0d0*PI*AS2) 
                             AS = as2
 
@@ -6342,8 +6303,9 @@
 own and set the path to its library in the MG5aMC option 'ninja'.""")
                     cmd.exec_cmd('set ninja None')
                     cmd.exec_cmd('save options')  
+
+
     # ==========================================================================
-
     # First treat the MadLoop5 standalone case       
     MadLoop_SA_options = {'clean': not noclean, 
       'complex_mass':cmd.options['complex_mass_scheme'],
@@ -6361,10 +6323,11 @@
       'output_dependencies':cmd.options['output_dependencies'],
       'SubProc_prefix':'P',
       'compute_color_flows':cmd.options['loop_color_flows'],
-      'mode': 'reweight' if cmd._export_format == "standalone_rw" else ''
+      'mode': 'reweight' if cmd._export_format == "standalone_rw" else '',
+      'cluster_local_path': cmd.options['cluster_local_path']
       }
 
-    if output_type.startswith('madloop'):
+    if output_type.startswith('madloop'):        
         import madgraph.loop.loop_exporters as loop_exporters
         if os.path.isdir(os.path.join(cmd._mgme_dir, 'Template/loop_material')):
             ExporterClass=None

=== modified file 'madgraph/iolibs/file_writers.py'
--- madgraph/iolibs/file_writers.py	2015-10-01 16:00:08 +0000
+++ madgraph/iolibs/file_writers.py	2016-03-16 17:59:03 +0000
@@ -185,7 +185,7 @@
     comment_char = 'c'
     downcase = False
     line_length = 71
-    max_split = 10
+    max_split = 20
     split_characters = "+-*/,) "
     comment_split_characters = " "
 
@@ -332,20 +332,46 @@
         res_lines = [line]