Merge lp:~maddevelopers/mg5amcnlo/python_syscalc into lp:~maddevelopers/mg5amcnlo/2.5.0

[Olivier Mattelaer]

This version allows to compute systematics error band on both LO and NLO sample.

For LO code, this is run automatically if 'use_syst' is set on True in the run_card.
For NLO code, this need to be called manually via the command 'systematics RUN_XX OPTIONS'
The sample need to be generated with the options 'store_rwgt_info=True'

The computation can also be called outside of madgraph via the code
madgraph/various/systematics.py INPUT.LHE[.GZ] OUTPUT.LHE[.GZ] OPTIONS

Here is the help of the systematics function such that you can see the options/make comment on this help if this is not clear/ comment on the default value.

INFO: syntax: systematics RUN_NAME [OUTPUT] [options]
INFO: -- Run the systematics run on the run_name run.
INFO: RUN_NAME can be a path to a lhef file.
INFO: OUTPUT can be the path to the output lhe file. we overwritte the input file otherwise
INFO:
INFO: options: (value written are default)
INFO:
INFO: --mur=0.5,1,2 # specify the values for renormalisation scale variation
INFO: --muf=0.5,1,2 # specify the values for factorisation scale variation
INFO: --alps=1 # specify the values for MLM emission scale variation
INFO: --dyn=-1,1,2,3,4 # specify the dynamical schemes to use.
INFO: # -1 is the one used by the sample.
INFO: # > 0 correspond to options of dynamical_scale_choice of the run_card.
INFO: --pdf=errorset # specify the pdfs to use for pdf variation. (see below)
INFO: --together=mur,muf,dyn # which parameter should we vary together to have all uncorrelated combination of the variations
INFO: --from_card # use the information from the run_card (LO only).
INFO:
INFO: Allowed value for the pdf options:
INFO: central : Do not perform any pdf variation
INFO: errorset : runs over the errorset
INFO: 244800 : runs over the associated set and his errorset
INFO: 244800@0 : runs over the associated set ONLY
INFO: CT10 : runs over the associated set and his errorset
INFO: CT10@0 : runs over the associated set ONLY
INFO: CT10@X : runs over the Xth set of the associated error set
INFO: XX,YY,ZZ : runs over the sets for XX,YY,ZZ (those three follows above syntax)

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/python_syscalc]

294. By Olivier Mattelaer on 2016-08-12

making the code runnning correctly on cluster mode

295. By Olivier Mattelaer on 2016-08-13

forget to push the last minut changes in the pdf syntax

296. By Olivier Mattelaer on 2016-08-17

merge with latest 2.5.0 version

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier,

a)

I tried running everything out of the box for p p > e+ ve and python itself
crashed (so the error was not handled), because I guess one cannot use your
py-syscalc without lhapdf. I got

terminate called after throwing an instance of 'LHAPDF::ReadError'
  what(): Info file not found for PDF set 'NNPDF23_lo_as_0130_qed'

So the case of lhapdf not being used should be better handled.

Also, in the run_card, the commentary of the entry

NNPDF23_nlo_as_0116_mc = sys_pdf # matching scales

is not correct. And it is not clear what the user is supposed to or can put
for this field. It seems that putting the LHAID works as well as putting
directly the name of the pdf. Is this also true when using the old C++
syscalc?
The commentary of this entry should clarify this point.

Also, when putting a PDF set in 'sys_pdf' that I don't have, the codes
doesn't try to download it (unless of course this is also the PDF set that
I want to use for my central value).
It would be a nice feature if it did, because as of now the code simply
crashes when doing the systematics.
There should at least be a check that the PDF is available.

This comment applies both for LO and NLO.

b)

The help of the command 'systematics' is nice. I would maybe rephrase the
message (which is tricky for users not used to that)

--together=mur,muf,dyn # which parameter should we vary together to have
all uncorrelated combination of the variations
to

'lists the parameter that must be varied simultaneously so as to compute
the weights for all combinations of their variations'.

but up to you.

c)

At NLO, I tried

systematics run_02 --pdf=21100 --together=mur,muf

on p p > e+ ve [QCD] and it crashed with:

WARNING: program
/opt/local/Library/Frameworks/Python.framework/Versions/2.7/Resources/Python.app/Contents/MacOS/Python
/Users/valentin/Documents/Work/MG5/python_syscalc/PROCNLO_loop_sm_0/bin/internal/systematics.py
events.lhe.gz ./tmp_1_events.lhe.gz --pdf=21100 --together=mur,muf
--start_event=2500 --stop_event=5000 --result=./log_sys_1.txt
--lhapdf_config=/Users/valentin/Documents/HEP_softs/LHAPDF/LHAPDF-6.1.5_installation/bin/lhapdf-config
launch ends with non zero status: 1. Stop all computation

There seems to be a 'launch' too many in the command no?

I tried to run it by hand without the 'launch', but I got:

    pdfset = lhapdf.mkPDF(int(data)).set()
  File "lhapdf.pyx", line 469, in lhapdf.mkPDF (lhapdf.cpp:7253)
  File "lhapdf.pyx", line 454, in lhapdf.mkPDF_lhaid (lhapdf.cpp:7044)
RuntimeError: Empty PDF file name given to Info::load

On Thu, Aug 11, 2016 at 4:55 PM, Olivier Mattelaer <
<email address hidden>> wrote:

> Olivier Mattelaer has proposed merging lp:~maddevelopers/mg5amcnlo/python_syscalc
> into lp:~maddevelopers/mg5amcnlo/2.5.0.
>
> Requested reviews:
> Eleni Vryonidou (evryonidou)
> Rikkert Frederix (frederix)
> Valentin Hirschi (valentin-hirschi)
>
> For more details, see:
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/
> python_syscalc/+merge/302735
>
> This version allows to compute systematics error band on both LO and NLO
> sample.
>
> For LO code, this is run automatically if 'use_syst' is set on True in the
> run_car...

[lp:~maddevelopers/mg5amcnlo/python_syscalc]

297. By Olivier Mattelaer on 2016-08-19

Fix related to the review of Valentin + fixing reweighting for FxFx

[Olivier Mattelaer (olivier-mattelaer)]

Hi Valentin,

Thanks for this review.

> So the case of lhapdf not being used should be better handled.

We actually do not need to run with lhapdf, just to have the lhapdf (version6) available.

> Also, in the run_card, the commentary of the entry
> NNPDF23_nlo_as_0116_mc = sys_pdf # matching scales
> is not correct

fixed

> It seems that putting the LHAID works as well as putting
> directly the name of the pdf. Is this also true when using the old C++
> syscalc?
> The commentary of this entry should clarify this point.

This is not a valid input for SysCalc. So I will not mention it here.

> Also, when putting a PDF set in 'sys_pdf' that I don't have, the codes
> doesn't try to download it (unless of course this is also the PDF set that
> I want to use for my central value).
> It would be a nice feature if it did, because as of now the code simply
> crashes when doing the systematics.

Done. I also automatically load the central if the run was not performed with lhapdf.
(knowing that nnpdf is likely to be slightly off due to their bug in lhapdf.)

> The help of the command 'systematics' is nice. I would maybe rephrase the
> message

Done

> pdfset = lhapdf.mkPDF(int(data)).set()
> File "lhapdf.pyx", line 469, in lhapdf.mkPDF (lhapdf.cpp:7253)
> File "lhapdf.pyx", line 454, in lhapdf.mkPDF_lhaid (lhapdf.cpp:7044)
> RuntimeError: Empty PDF file name given to Info::load

I do not reproduce this. Looks like either some PDF sets are missing or (worse) corrupted
I have face as well some problem with some lhapdf set so I have use another method
to load them and it fixes the problem that i faced.

Cheers,

Olivier

> On Aug 19, 2016, at 07:21, Valentin Hirschi <email address hidden> wrote:
>
> Hi Olivier,
>
> a)
>
> I tried running everything out of the box for p p > e+ ve and python itself
> crashed (so the error was not handled), because I guess one cannot use your
> py-syscalc without lhapdf. I got
>
> terminate called after throwing an instance of 'LHAPDF::ReadError'
> what(): Info file not found for PDF set 'NNPDF23_lo_as_0130_qed'
>
> So the case of lhapdf not being used should be better handled.
>
> Also, in the run_card, the commentary of the entry
>
> NNPDF23_nlo_as_0116_mc = sys_pdf # matching scales
>
> is not correct. And it is not clear what the user is supposed to or can put
> for this field. It seems that putting the LHAID works as well as putting
> directly the name of the pdf. Is this also true when using the old C++
> syscalc?
> The commentary of this entry should clarify this point.
>
> Also, when putting a PDF set in 'sys_pdf' that I don't have, the codes
> doesn't try to download it (unless of course this is also the PDF set that
> I want to use for my central value).
> It would be a nice feature if it did, because as of now the code simply
> crashes when doing the systematics.
> There should at least be a check that the PDF is available.
>
> This comment applies both for LO and NLO.
>
> b)
>
> The help of the command 'systematics' is nice. I would maybe rephrase the
> message (which is tricky for users not used to that)
>
> --together=mur,muf,dyn # which parameter should we vary together to hav...

[Valentin Hirschi (valentin-hirschi)]

> Hi Valentin,
>
> Thanks for this review.
>
> > So the case of lhapdf not being used should be better handled.
>
>
> We actually do not need to run with lhapdf, just to have the lhapdf (version6)
> available.

I know what I meant was that if I don't have lhapdf (or only lhapdf5) on my system, then systematics crashes with the following message:

UnboundLocalError : local variable 'lhapdf_libdir' referenced before assignment

( File "/Users/valentin/Documents/Work/MG5/python_syscalc/madgraph/various/misc.py", line 1637, in import_python_lhapdf
    candidates=[dirname for dirname in os.listdir(lhapdf_libdir+'64') \)

It would be nice if this situation could be detected before the crash and if a message would pop up warning the user that systematics are bypassed because he doesn't have lhapdf and that once he will have installed it, he can still rerun systematics independently on his results already produced by using the interactive madevent interface.

>
> > Also, in the run_card, the commentary of the entry
> > NNPDF23_nlo_as_0116_mc = sys_pdf # matching scales
> > is not correct
>
>
> fixed
>
> > It seems that putting the LHAID works as well as putting
> > directly the name of the pdf. Is this also true when using the old C++
> > syscalc?
> > The commentary of this entry should clarify this point.
>
> This is not a valid input for SysCalc. So I will not mention it here.

Ok, and it would be indeed to complicated to explain in which situation the user can put an lhaid there, I agree. Too bad.

> > Also, when putting a PDF set in 'sys_pdf' that I don't have, the codes
> > doesn't try to download it (unless of course this is also the PDF set that
> > I want to use for my central value).
> > It would be a nice feature if it did, because as of now the code simply
> > crashes when doing the systematics.
>
> Done. I also automatically load the central if the run was not performed with
> lhapdf.
> (knowing that nnpdf is likely to be slightly off due to their bug in lhapdf.)

Excellent, I tested it and it works.

> > The help of the command 'systematics' is nice. I would maybe rephrase the
> > message
>
> Done
>
> > pdfset = lhapdf.mkPDF(int(data)).set()
> > File "lhapdf.pyx", line 469, in lhapdf.mkPDF (lhapdf.cpp:7253)
> > File "lhapdf.pyx", line 454, in lhapdf.mkPDF_lhaid (lhapdf.cpp:7044)
> > RuntimeError: Empty PDF file name given to Info::load
>
> I do not reproduce this. Looks like either some PDF sets are missing or
> (worse) corrupted
> I have face as well some problem with some lhapdf set so I have use another
> method
> to load them and it fixes the problem that i faced.

Ok, now at NLO, I crash on:

AttributeError : 'Event' object has no attribute 'nloweight'
[...]
  File "/Users/valentin/Documents/Work/MG5/python_syscalc/madgraph/various/systematics.py", line 583, in get_nlo_wgt
    nloinfo = event.parse_nlo_weight()
  File "/Users/valentin/Documents/Work/MG5/python_syscalc/madgraph/various/lhe_parser.py", line 1134, in parse_nlo_weight
    return self.nloweight
AttributeError: 'Event' object has no attribute 'nloweight'

This happened both with and without PDF variations in my generation (I assume you require the user t...

[Olivier Mattelaer (olivier-mattelaer)]

> I know what I meant was that if I don't have lhapdf (or only lhapdf5) on my
> system, then systematics crashes with the following message:

ok this should be fixed now.

> Ok, now at NLO, I crash on:
> AttributeError: 'Event' object has no attribute 'nloweight'
>

This is because, you didnot put 'store_rwgt_info' on True in the NLo run_card.
Impossible to compute systematics in that case.
I have put in place a nicer error message for that situation.

> This happened both with and without PDF variations in my generation (I assume
> you require the user to having generated events with both variations turned
> on, it would be great if you could safe-check this).

No you can have both variation turned off.

> By the way, additional general question:
> If there is already a bunch of weights in the event, but only some are missing
> given the inputs provided to 'systematics', will you recompute all of them or
> do you automatically filter them out?
> It's ok if you recompute all of them, but I was wondering if you had optimized
> this bit.

No I'm not smart at all on that point of view, I recompute them and add them to the event file.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/python_syscalc]

298. By Olivier Mattelaer on 2016-08-20

apply change related to round2 of Valentin review

[Eleni Vryonidou (evryonidou)]

Hi Olivier,

I tried getting a fresh branch of python_syscalc but trying to run fails with this error:

ImportError : No module named various.systematics

And in the log file:

 File "/Users/evryonidou/Desktop/python_syscalc/check/bin/internal/common_run_interface.py", line 1336, in do_systematics
    import internal.various.systematics as systematics
ImportError: No module named various.systematics

So I couldn’t go further with any other checks… Do you have any idea what this problem is?

Cheers,

Eleni

> On 12 Aug 2016, at 01:55, Olivier Mattelaer <email address hidden> wrote:
>
> Olivier Mattelaer has proposed merging lp:~maddevelopers/mg5amcnlo/python_syscalc into lp:~maddevelopers/mg5amcnlo/2.5.0.
>
> Requested reviews:
> Eleni Vryonidou (evryonidou)
> Rikkert Frederix (frederix)
> Valentin Hirschi (valentin-hirschi)
>
> For more details, see:
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/python_syscalc/+merge/302735
>
> This version allows to compute systematics error band on both LO and NLO sample.
>
> For LO code, this is run automatically if 'use_syst' is set on True in the run_card.
> For NLO code, this need to be called manually via the command 'systematics RUN_XX OPTIONS'
> The sample need to be generated with the options 'store_rwgt_info=True'
>
> The computation can also be called outside of madgraph via the code
> madgraph/various/systematics.py INPUT.LHE[.GZ] OUTPUT.LHE[.GZ] OPTIONS
>
> Here is the help of the systematics function such that you can see the options/make comment on this help if this is not clear/ comment on the default value.
>
> INFO: syntax: systematics RUN_NAME [OUTPUT] [options]
> INFO: -- Run the systematics run on the run_name run.
> INFO: RUN_NAME can be a path to a lhef file.
> INFO: OUTPUT can be the path to the output lhe file. we overwritte the input file otherwise
> INFO:
> INFO: options: (value written are default)
> INFO:
> INFO: --mur=0.5,1,2 # specify the values for renormalisation scale variation
> INFO: --muf=0.5,1,2 # specify the values for factorisation scale variation
> INFO: --alps=1 # specify the values for MLM emission scale variation
> INFO: --dyn=-1,1,2,3,4 # specify the dynamical schemes to use.
> INFO: # -1 is the one used by the sample.
> INFO: # > 0 correspond to options of dynamical_scale_choice of the run_card.
> INFO: --pdf=errorset # specify the pdfs to use for pdf variation. (see below)
> INFO: --together=mur,muf,dyn # which parameter should we vary together to have all uncorrelated combination of the variations
> INFO: --from_card # use the information from the run_card (LO only).
> INFO:
> INFO: Allowed value for the pdf options:
> INFO: central : Do not perform any pdf variation
> INFO: errorset : runs over the errorset
> INFO: 244800 : runs over the associated set and his errorset
> INFO: 244800@0 : runs over the associated set ONLY
> INFO: CT10 : runs over the associated set and his errorset
> INFO: CT10@0 : runs over the associated set ONLY
> INFO: CT10@X : runs over the Xth set of ...

[lp:~maddevelopers/mg5amcnlo/python_syscalc]

299. By Olivier Mattelaer on 2016-08-22

fix an import problem

[Olivier Mattelaer (olivier-mattelaer)]

Thanks to have spot this.

It should be fix now.

Cheers,

Olivier
> On Aug 22, 2016, at 18:30, Eleni Vryonidou <email address hidden> wrote:
>
> Hi Olivier,
>
> I tried getting a fresh branch of python_syscalc but trying to run fails with this error:
>
> ImportError : No module named various.systematics
>
> And in the log file:
>
> File "/Users/evryonidou/Desktop/python_syscalc/check/bin/internal/common_run_interface.py", line 1336, in do_systematics
> import internal.various.systematics as systematics
> ImportError: No module named various.systematics
>
> So I couldn’t go further with any other checks… Do you have any idea what this problem is?
>
> Cheers,
>
> Eleni
>
>
>> On 12 Aug 2016, at 01:55, Olivier Mattelaer <email address hidden> wrote:
>>
>> Olivier Mattelaer has proposed merging lp:~maddevelopers/mg5amcnlo/python_syscalc into lp:~maddevelopers/mg5amcnlo/2.5.0.
>>
>> Requested reviews:
>> Eleni Vryonidou (evryonidou)
>> Rikkert Frederix (frederix)
>> Valentin Hirschi (valentin-hirschi)
>>
>> For more details, see:
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/python_syscalc/+merge/302735
>>
>> This version allows to compute systematics error band on both LO and NLO sample.
>>
>> For LO code, this is run automatically if 'use_syst' is set on True in the run_card.
>> For NLO code, this need to be called manually via the command 'systematics RUN_XX OPTIONS'
>> The sample need to be generated with the options 'store_rwgt_info=True'
>>
>> The computation can also be called outside of madgraph via the code
>> madgraph/various/systematics.py INPUT.LHE[.GZ] OUTPUT.LHE[.GZ] OPTIONS
>>
>> Here is the help of the systematics function such that you can see the options/make comment on this help if this is not clear/ comment on the default value.
>>
>> INFO: syntax: systematics RUN_NAME [OUTPUT] [options]
>> INFO: -- Run the systematics run on the run_name run.
>> INFO: RUN_NAME can be a path to a lhef file.
>> INFO: OUTPUT can be the path to the output lhe file. we overwritte the input file otherwise
>> INFO:
>> INFO: options: (value written are default)
>> INFO:
>> INFO: --mur=0.5,1,2 # specify the values for renormalisation scale variation
>> INFO: --muf=0.5,1,2 # specify the values for factorisation scale variation
>> INFO: --alps=1 # specify the values for MLM emission scale variation
>> INFO: --dyn=-1,1,2,3,4 # specify the dynamical schemes to use.
>> INFO: # -1 is the one used by the sample.
>> INFO: # > 0 correspond to options of dynamical_scale_choice of the run_card.
>> INFO: --pdf=errorset # specify the pdfs to use for pdf variation. (see below)
>> INFO: --together=mur,muf,dyn # which parameter should we vary together to have all uncorrelated combination of the variations
>> INFO: --from_card # use the information from the run_card (LO only).
>> INFO:
>> INFO: Allowed value for the pdf options:
>> INFO: central : Do not perform any pdf variation
>> INFO: errorset : runs over the errorset
>> INFO: 244800 : ...

[Eleni Vryonidou (evryonidou)]

Hi,

ok that works now.

At NLO if I generate events with the default pdf set (nn23nlo) it crashes with the following error:

PROCNLO_loop_sm_0>systematics run_02
INFO: Running Systematic computation
INFO: Idle: 1, Running: 3, Completed: 0 [ current time: 17h12 ]
WARNING: program /opt/local/Library/Frameworks/Python.framework/Versions/2.7/Resources/Python.app/Contents/MacOS/Python /Users/evryonidou/Desktop/python_syscalc/PROCNLO_loop_sm_0/bin/internal/systematics.py events.lhe.gz ./tmp_0_events.lhe.gz --start_event=0 --stop_event=2500 --result=./log_sys_0.txt --lhapdf_config=lhapdf-config launch ends with non zero status: 1. Stop all computation
WARNING: program /opt/local/Library/Frameworks/Python.framework/Versions/2.7/Resources/Python.app/Contents/MacOS/Python /Users/evryonidou/Desktop/python_syscalc/PROCNLO_loop_sm_0/bin/internal/systematics.py events.lhe.gz ./tmp_1_events.lhe.gz --start_event=2500 --stop_event=5000 --result=./log_sys_1.txt --lhapdf_config=lhapdf-config launch ends with non zero status: 1. Stop all computation
Command "systematics run_02" interrupted with error:
Exception : program /opt/local/Library/Frameworks/Python.framework/Versions/2.7/Resources/Python.app/Contents/MacOS/Python /Users/evryonidou/Desktop/python_syscalc/PROCNLO_loop_sm_0/bin/internal/systematics.py events.lhe.gz ./tmp_0_events.lhe.gz --start_event=0 --stop_event=2500 --result=./log_sys_0.txt --lhapdf_config=lhapdf-config launch ends with non zero status: 1. Stop all computation

If I generate events with lhapdf then it works fine.

Cheers,

Eleni

> On 22 Aug 2016, at 21:00, Olivier Mattelaer <email address hidden> wrote:
>
> Thanks to have spot this.
>
> It should be fix now.
>
> Cheers,
>
> Olivier
>> On Aug 22, 2016, at 18:30, Eleni Vryonidou <email address hidden> wrote:
>>
>> Hi Olivier,
>>
>> I tried getting a fresh branch of python_syscalc but trying to run fails with this error:
>>
>> ImportError : No module named various.systematics
>>
>> And in the log file:
>>
>> File "/Users/evryonidou/Desktop/python_syscalc/check/bin/internal/common_run_interface.py", line 1336, in do_systematics
>> import internal.various.systematics as systematics
>> ImportError: No module named various.systematics
>>
>> So I couldn’t go further with any other checks… Do you have any idea what this problem is?
>>
>> Cheers,
>>
>> Eleni
>>
>>
>>> On 12 Aug 2016, at 01:55, Olivier Mattelaer <email address hidden> wrote:
>>>
>>> Olivier Mattelaer has proposed merging lp:~maddevelopers/mg5amcnlo/python_syscalc into lp:~maddevelopers/mg5amcnlo/2.5.0.
>>>
>>> Requested reviews:
>>> Eleni Vryonidou (evryonidou)
>>> Rikkert Frederix (frederix)
>>> Valentin Hirschi (valentin-hirschi)
>>>
>>> For more details, see:
>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/python_syscalc/+merge/302735
>>>
>>> This version allows to compute systematics error band on both LO and NLO sample.
>>>
>>> For LO code, this is run automatically if 'use_syst' is set on True in the run_card.
>>> For NLO code, this need to be called manually via the command 'systematics RUN_XX OPTIONS'
>>> The sample need to be generated with the options 'store_...

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

This is actually a bug of the lhapdf nn23nlo set.

*** skip this if you do not care to know what the bug is ***

More in details, the version from MG5 is not the same as the one from lhapdf.
Turns out that the one from MG5 is the correct one and that they have to update their set on lhapdf.
Indeed in lhapdf, the central set was generated with a fit for the options ‘always positive’ but the error set was fitted for the options ’can be negative’.
so the set is inconsistent. In MG5aMC, the central set is ‘can be negative’ which is the central value of the associated error set.
The LHAPDF author/NNPDF author are aware of this, but I do not know when and how this error will be solved. (The bug was actually spotted due to this project).
Since this will pass by an update of the pdf, it is likely that people will not perform the update of the set for a long time…

*** end of the bug description ***

I just think of a new work-around which is not perfect but good enough
1) if the central pdf is zero, then use the wgt from the event to get the pdf value.
2) if the pdf is non zero, use it.
in debug mode, this the second will need to a crash since the code will realize that a mismatch happens.
I have updated the code such that in production mode, the test/crash does not happen.

Cheers,

Olivier

PS: I have made a test from re-computing the central pdf from the errorset, and this slow down the code by huge factor (20 times slower) and provide the same result as the above.
PPS: if you want to test, please run in production mode with python -O before the executable.

> On Aug 23, 2016, at 16:24, Eleni Vryonidou <email address hidden> wrote:
>
> Hi,
>
> ok that works now.
>
> At NLO if I generate events with the default pdf set (nn23nlo) it crashes with the following error:
>
> PROCNLO_loop_sm_0>systematics run_02
> INFO: Running Systematic computation
> INFO: Idle: 1, Running: 3, Completed: 0 [ current time: 17h12 ]
> WARNING: program /opt/local/Library/Frameworks/Python.framework/Versions/2.7/Resources/Python.app/Contents/MacOS/Python /Users/evryonidou/Desktop/python_syscalc/PROCNLO_loop_sm_0/bin/internal/systematics.py events.lhe.gz ./tmp_0_events.lhe.gz --start_event=0 --stop_event=2500 --result=./log_sys_0.txt --lhapdf_config=lhapdf-config launch ends with non zero status: 1. Stop all computation
> WARNING: program /opt/local/Library/Frameworks/Python.framework/Versions/2.7/Resources/Python.app/Contents/MacOS/Python /Users/evryonidou/Desktop/python_syscalc/PROCNLO_loop_sm_0/bin/internal/systematics.py events.lhe.gz ./tmp_1_events.lhe.gz --start_event=2500 --stop_event=5000 --result=./log_sys_1.txt --lhapdf_config=lhapdf-config launch ends with non zero status: 1. Stop all computation
> Command "systematics run_02" interrupted with error:
> Exception : program /opt/local/Library/Frameworks/Python.framework/Versions/2.7/Resources/Python.app/Contents/MacOS/Python /Users/evryonidou/Desktop/python_syscalc/PROCNLO_loop_sm_0/bin/internal/systematics.py events.lhe.gz ./tmp_0_events.lhe.gz --start_event=0 --stop_event=2500 --result=./log_sys_0.txt --lhapdf_config=lhapdf-config launch ends with non zero status: 1. Stop all computation
>
> If I...

[lp:~maddevelopers/mg5amcnlo/python_syscalc]

300. By Olivier Mattelaer on 2016-08-23

put in place a nice workaround in presence of the nnpdf23 bug

301. By Olivier Mattelaer on 2016-08-23

add security to check that the card are consistent with the systematic computation

302. By Olivier Mattelaer on 2016-08-24

fixing a bug in the reweight module (incorrect relative path)

303. By Olivier Mattelaer on 2016-08-24

remove a print comment

[Eleni Vryonidou (evryonidou)]

review: Approve

[Rikkert Frederix (frederix)]

Hi Olivier,

Ok. I finally had some time to look into this branch. I haven't really tested anything yet, but looking through the changes made in this branch, I noticed that there might be a mistake related to the possible changes of the dynamical scale setting 'dyn'. In particular, it looks like you simply recompute if from the momenta as given in the event file. This is not correct, because this gives different scales to the (MC) subtraction terms in the S and H events, even if those subtraction terms were computed with exactly the same kinematics.

What we do in normal running of the code is the following:

n-body kinematics:
- Born
- Virtual
- (Integrated) FKS Subtraction terms

(n+1)-body kinematics:
- MC subtraction terms (in S and H events)
- Real-emission

See also the calls to 'set_alphaS' in driver_mintMC where the scales are either set using n-body kinematics ('call set_alphaS(p1_cnt(0,1,0))') or (n+1)-body kinematics ('call set_alphaS(p)'), with the subsequent calls to the matrix elements just after those.

I think this needs to be changed, even though in practice this difference is numerically usually small. But to get it 100% correct these changes are needed.

Cheers,
Rik

review: Approve

[Rikkert Frederix (frederix)]

review: Needs Fixing

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2015-08-11 00:35:53 +0000
+++ Template/LO/Cards/run_card.dat	2016-08-24 22:15:21 +0000
@@ -279,5 +279,5 @@
 %(sys_alpsfact)s = sys_alpsfact  # \alpha_s emission scale factors
 %(sys_matchscale)s = sys_matchscale # variation of merging scale
 # PDF sets and number of members (0 or none for all members).
-%(sys_pdf)s = sys_pdf # matching scales
+%(sys_pdf)s = sys_pdf # separate by && if more than one set.
 # MSTW2008nlo68cl.LHgrid 1  = sys_pdf
 
=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2016-06-07 09:33:55 +0000
+++ Template/LO/SubProcesses/reweight.f	2016-08-24 22:15:21 +0000
@@ -994,21 +994,21 @@
             pt2ijcl(nexternal-2)=q2fact(1)
             if(nexternal.gt.3) pt2ijcl(nexternal-3)=q2fact(1)
          else
-            q2fact(1)=pt2ijcl(nexternal-2)
-            q2fact(2)=q2fact(1)
+            q2fact(1)=scalefact**2*pt2ijcl(nexternal-2)
+            q2fact(2)=scalefact**2*q2fact(1)
          endif
       elseif(jcentral(1).eq.0)then
-            q2fact(1) = pt2ijcl(jfirst(1))
+            q2fact(1) = scalefact**2*pt2ijcl(jfirst(1))
       elseif(jcentral(2).eq.0)then
-            q2fact(2) = pt2ijcl(jfirst(2))
+            q2fact(2) = scalefact**2*pt2ijcl(jfirst(2))
       elseif(ickkw.eq.2.or.(pdfwgt.and.ickkw.gt.0))then
 c     Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
 c     f1(x1,pt2E) is given by DSIG, just need to set scale.
 c     Use the minimum scale found for fact scale in ME
          if(jlast(1).gt.0.and.jfirst(1).le.jlast(1))
-     $        q2fact(1)=min(pt2ijcl(jfirst(1)),q2fact(1))
+     $        q2fact(1)=scalefact**2*min(pt2ijcl(jfirst(1)),q2fact(1))
          if(jlast(2).gt.0.and.jfirst(2).le.jlast(2))
-     $        q2fact(2)=min(pt2ijcl(jfirst(2)),q2fact(2))
+     $        q2fact(2)=scalefact**2*min(pt2ijcl(jfirst(2)),q2fact(2))
       endif
 
 c     Check that factorization scale is >= 2 GeV
 
=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2016-08-17 19:36:40 +0000
+++ UpdateNotes.txt	2016-08-24 22:15:21 +0000
@@ -1,5 +1,6 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
+<<<<<<< TREE
 2.5.0
        OM: Modify the structure of the output format such that all the internal format have the same structure
        OM: Adding the Plugin directory. Three kind of plugin are currently supported
@@ -15,6 +16,20 @@
        OM: add the check that the param_card is compatible with the model restriction.
        OM+VH: 
      
+=======
+
+X.X.X (xx/xx/xx)
+        OM: NLO/LO Re-weighting works in multi-core
+        OM: Introduces a new function for LO/NLO interface "systematics"
+            This function allows to compute systematics uncertainty from the event sample
+            It requires the event sample to have been generated with 
+                - use_syst = T (for LO sample)
+                - store_reweight_info = T (for NLO sample)
+            At LO the code is run automatically if use_syst=T (but if SysCalc is installed)
+        OM: Fix a bug in the reweight_card where relative path was not understood from the local directory 
+            where the program was runned by the user.
+
+>>>>>>> MERGE-SOURCE
 2.4.3 (01/08/16)
         OM: Reduce the amount of log file/output generated for LO run (output can use up to three times less output).
         OM: For the LO combination of events (unweighting) pass to the method previously used for loop-induced.
 
=== modified file 'VERSION'
--- VERSION	2016-08-03 13:42:53 +0000
+++ VERSION	2016-08-24 22:15:21 +0000
@@ -1,5 +1,11 @@
+<<<<<<< TREE
 version = 2.5.0.alpha
 date = 2016-08-01
+=======
+version = new_sys
+date = 2016-08-08
+
+>>>>>>> MERGE-SOURCE
 
 
 
 
=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2016-08-17 16:25:19 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2016-08-24 22:15:21 +0000
@@ -3297,98 +3297,6 @@
         self.update_status('Run complete', level='shower', update_results=True)
 
 
-    ############################################################################
-    def set_run_name(self, name, tag=None, level='parton', reload_card=False):
-        """define the run name, the run_tag, the banner and the results."""
-        
-        # when are we force to change the tag new_run:previous run requiring changes
-        upgrade_tag = {'parton': ['parton','pythia','pgs','delphes','shower'],
-                       'pythia': ['pythia','pgs','delphes'],
-                       'shower': ['shower'],
-                       'pgs': ['pgs'],
-                       'delphes':['delphes'],
-                       'plot':[]}
-        
-        
-
-        if name == self.run_name:        
-            if reload_card:
-                run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
-                self.run_card = banner_mod.RunCardNLO(run_card)
-
-            #check if we need to change the tag
-            if tag:
-                self.run_card['run_tag'] = tag
-                self.run_tag = tag
-                self.results.add_run(self.run_name, self.run_card)
-            else:
-                for tag in upgrade_tag[level]:
-                    if getattr(self.results[self.run_name][-1], tag):
-                        tag = self.get_available_tag()
-                        self.run_card['run_tag'] = tag
-                        self.run_tag = tag
-                        self.results.add_run(self.run_name, self.run_card)                        
-                        break
-            return # Nothing to do anymore
-        
-        # save/clean previous run
-        if self.run_name:
-            self.store_result()
-        # store new name
-        self.run_name = name
-        
-        # Read run_card
-        run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
-        self.run_card = banner_mod.RunCardNLO(run_card)
-
-        new_tag = False
-        # First call for this run -> set the banner
-        self.banner = banner_mod.recover_banner(self.results, level, self.run_name, tag)
-        if tag:
-            self.run_card['run_tag'] = tag
-            new_tag = True
-        elif not self.run_name in self.results and level =='parton':
-            pass # No results yet, so current tag is fine
-        elif not self.run_name in self.results:
-            #This is only for case when you want to trick the interface
-            logger.warning('Trying to run data on unknown run.')
-            self.results.add_run(name, self.run_card)
-            self.results.update('add run %s' % name, 'all', makehtml=True)
-        else:
-            for tag in upgrade_tag[level]:
-                
-                if getattr(self.results[self.run_name][-1], tag):
-                    # LEVEL is already define in the last tag -> need to switch tag
-                    tag = self.get_available_tag()
-                    self.run_card['run_tag'] = tag
-                    new_tag = True
-                    break
-            if not new_tag:
-                # We can add the results to the current run
-                tag = self.results[self.run_name][-1]['tag']
-                self.run_card['run_tag'] = tag # ensure that run_tag is correct                
-             
-                    
-        if name in self.results and not new_tag:
-            self.results.def_current(self.run_name)
-        else:
-            self.results.add_run(self.run_name, self.run_card)
-
-        self.run_tag = self.run_card['run_tag']
-
-        # Return the tag of the previous run having the required data for this
-        # tag/run to working wel.
-        if level == 'parton':
-            return
-        elif level == 'pythia':
-            return self.results[self.run_name][0]['tag']
-        else:
-            for i in range(-1,-len(self.results[self.run_name])-1,-1):
-                tagRun = self.results[self.run_name][i]
-                if tagRun.pythia:
-                    return tagRun['tag']
-
-
     def store_result(self):
         """ tar the pythia results. This is done when we are quite sure that 
         the pythia output will not be use anymore """
@@ -4928,9 +4836,9 @@
             if '--web' in args:
                 i = args.index('--web') 
                 args.pop(i)                                                                                                                                                                     
-                cmd_line =  aMCatNLOCmd(force_run=True)
+                cmd_line =  aMCatNLOCmd(me_dir=os.path.dirname(root_path),force_run=True)
             else:
-                cmd_line =  aMCatNLOCmdShell(force_run=True)
+                cmd_line =  aMCatNLOCmdShell(me_dir=os.path.dirname(root_path),force_run=True)
 
             if not hasattr(cmd_line, 'do_%s' % args[0]):
                 if parser_error:
 
=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2016-08-17 16:25:19 +0000
+++ madgraph/interface/common_run_interface.py	2016-08-24 22:15:21 +0000
@@ -68,6 +68,7 @@
     import internal.files as files
     import internal.save_load_object as save_load_object
     import internal.gen_crossxhtml as gen_crossxhtml
+    import internal.lhe_parser as lhe_parser
     from internal import InvalidCmd, MadGraph5Error
     MADEVENT=True    
 else:
@@ -78,6 +79,7 @@
     import madgraph.various.misc as misc
     import madgraph.iolibs.files as files
     import madgraph.various.cluster as cluster
+    import madgraph.various.lhe_parser as lhe_parser
     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
@@ -450,7 +452,6 @@
         """return the path to the output events
         """
 
-
         if self.mode == 'madevent':
             possible_path = [
                 pjoin(self.me_dir,'Events', run_name, 'unweighted_events.lhe.gz'),
@@ -465,7 +466,10 @@
                 correct_path = path
                 break
         else:
-            raise self.InvalidCmd('No events file corresponding to %s run. ' % run_name)
+            if os.path.exists(run_name):
+                correct_path = run_name
+            else:
+                raise self.InvalidCmd('No events file corresponding to %s run. ' % run_name)
         return correct_path
 
 
@@ -619,6 +623,101 @@
 
         return args
 
+    ############################################################################
+    def set_run_name(self, name, tag=None, level='parton', reload_card=False,
+                     allow_new_tag=True):
+        """define the run name, the run_tag, the banner and the results."""
+        
+        # when are we force to change the tag new_run:previous run requiring changes
+        upgrade_tag = {'parton': ['parton','pythia','pgs','delphes'],
+                       'pythia': ['pythia','pgs','delphes'],
+                       'pgs': ['pgs'],
+                       'delphes':['delphes'],
+                       'plot':[],
+                       'syscalc':[]}
+        
+        
+
+        if name == self.run_name:       
+            if reload_card:
+                run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
+                self.run_card = banner_mod.RunCard(run_card)
+                
+            #check if we need to change the tag
+            if tag:
+                self.run_card['run_tag'] = tag
+                self.run_tag = tag
+                self.results.add_run(self.run_name, self.run_card)
+            else:                
+                for tag in upgrade_tag[level]:
+                    if getattr(self.results[self.run_name][-1], tag):
+                        tag = self.get_available_tag()
+                        self.run_card['run_tag'] = tag
+                        self.run_tag = tag
+                        self.results.add_run(self.run_name, self.run_card)                        
+                        break
+            return # Nothing to do anymore
+        
+        # save/clean previous run
+        if self.run_name:
+            self.store_result()
+        # store new name
+        self.run_name = name
+        
+        new_tag = False
+        # First call for this run -> set the banner
+        self.banner = banner_mod.recover_banner(self.results, level, name)
+        if 'mgruncard' in self.banner:
+            self.run_card = self.banner.charge_card('run_card')
+        else:
+            # Read run_card
+            run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
+            self.run_card = banner_mod.RunCard(run_card)   
+        
+        if tag:
+            self.run_card['run_tag'] = tag
+            new_tag = True
+        elif not self.run_name in self.results and level =='parton':
+            pass # No results yet, so current tag is fine
+        elif not self.run_name in self.results:
+            #This is only for case when you want to trick the interface
+            logger.warning('Trying to run data on unknown run.')
+            self.results.add_run(name, self.run_card)
+            self.results.update('add run %s' % name, 'all', makehtml=False)
+        else:
+            for tag in upgrade_tag[level]:
+                
+                if getattr(self.results[self.run_name][-1], tag):
+                    # LEVEL is already define in the last tag -> need to switch tag
+                    tag = self.get_available_tag()
+                    self.run_card['run_tag'] = tag
+                    new_tag = True
+                    break
+            if not new_tag:
+                # We can add the results to the current run
+                tag = self.results[self.run_name][-1]['tag']
+                self.run_card['run_tag'] = tag # ensure that run_tag is correct                
+                   
+        if allow_new_tag and (name in self.results and not new_tag):
+            self.results.def_current(self.run_name)
+        else:
+            self.results.add_run(self.run_name, self.run_card)
+
+        self.run_tag = self.run_card['run_tag']
+
+        # Return the tag of the previous run having the required data for this
+        # tag/run to working wel.
+        if level == 'parton':
+            return
+        elif level == 'pythia':
+            return self.results[self.run_name][0]['tag']
+        else:
+            for i in range(-1,-len(self.results[self.run_name])-1,-1):
+                tagRun = self.results[self.run_name][i]
+                if tagRun.pythia:
+                    return tagRun['tag']
+            
+
     @misc.multiple_try(nb_try=5, sleep=2)
     def load_results_db(self):
         """load the current results status"""
@@ -1087,7 +1186,270 @@
         """Dummy routine, to be overwritten by daughter classes"""
 
         pass
-    
+
+    ############################################################################
+    def help_systematics(self):
+        """help for systematics command"""
+        logger.info("syntax: systematics RUN_NAME [OUTPUT] [options]",'$MG:color:BLACK')
+        logger.info("-- Run the systematics run on the run_name run.")
+        logger.info("   RUN_NAME can be a path to a lhef file.")
+        logger.info("   OUTPUT can be the path to the output lhe file. we overwritte the input file otherwise") 
+        logger.info("")
+        logger.info("options: (value written are default)", '$MG:color:BLACK')
+        logger.info("")
+        logger.info("   --mur=0.5,1,2    # specify the values for renormalisation scale variation")
+        logger.info("   --muf=0.5,1,2    # specify the values for factorisation scale variation")
+        logger.info("   --alps=1         # specify the values for MLM emission scale variation")
+        logger.info("   --dyn=-1,1,2,3,4 # specify the dynamical schemes to use.")
+        logger.info("                    #   -1 is the one used by the sample.")
+        logger.info("                    #   > 0 correspond to options of dynamical_scale_choice of the run_card.")
+        logger.info("   --pdf=errorset   # specify the pdfs to use for pdf variation. (see below)")
+        logger.info("   --together=mur,muf,dyn # lists the parameter that must be varied simultaneously so as to ")
+        logger.info("                          # compute the weights for all combinations of their variations.")
+        logger.info("   --from_card      # use the information from the run_card (LO only).")
+        logger.info("")
+        logger.info("   Allowed value for the pdf options:", '$MG:color:BLACK')
+        logger.info("       central  : Do not perform any pdf variation"    )
+        logger.info("       errorset : runs over the  errorset")
+        logger.info("       244800   : runs over the associated set and his errorset")
+        logger.info("       244800@0 : runs over the associated set ONLY")
+#        logger.info("       244800@X : runs over the Xth set of the associated error set")
+        logger.info("       CT10     : runs over the associated set and his errorset")
+        logger.info("       CT10@0   : runs over the associated set ONLY")
+        logger.info("       CT10@X   : runs over the Xth set of the associated error set")
+        logger.info("       XX,YY,ZZ : runs over the sets for XX,YY,ZZ (those three follows above syntax)")
+        
+    def complete_systematics(self, text, line, begidx, endidx):
+        """auto completion for the systematics command"""
+ 
+        args = self.split_arg(line[0:begidx], error=False)
+        options = ['--mur=', '--muf=', '--pdf=', '--dyn=','--alps=','--together=','--from_card ']
+        
+        if len(args) == 1 and os.path.sep not in text:
+            #return valid run_name
+            data = misc.glob(pjoin('*','*events.lhe*'), pjoin(self.me_dir, 'Events'))
+            data = [n.rsplit('/',2)[1] for n in data]
+            return  self.list_completion(text, data, line)
+        elif len(args)==1:
+            #logger.warning('1args')
+            return self.path_completion(text,
+                                        os.path.join('.',*[a for a in args \
+                                                    if a.endswith(os.path.sep)]))
+        elif len(args)==2 and os.path.sep in args[1]:
+            #logger.warning('2args %s', args[1])
+            return self.path_completion(text, '.')
+              
+        elif not line.endswith(tuple(options)):
+            return self.list_completion(text, options)
+        
+         
+    ############################################################################
+    def do_systematics(self, line):
+        """ syntax is 'systematics [INPUT [OUTPUT]] OPTIONS'
+            --mur=0.5,1,2
+            --muf=0.5,1,2
+            --alps=1
+            --dyn=-1
+            --together=mur,muf #can be repeated
+            
+            #special options
+            --from_card=
+        """
+    
+        lhapdf = misc.import_python_lhapdf(self.options['lhapdf'])
+        if not lhapdf:
+            logger.info('can not run systematics since can not link python to lhapdf')
+            return
+    
+        self.update_status('Running Systematic computation', level='parton')
+        args = self.split_arg(line)
+        #split arguments and option
+        opts= []
+        args = [a for a in args if not a.startswith('-') or opts.append(a)] 
+
+        #check sanity of options
+        if any(not o.startswith(('--mur=', '--muf=', '--alps=','--dyn=','--together=','--from_card','--pdf='))
+                for o in opts):
+            raise self.InvalidCmd, "command systematics called with invalid option syntax. Please retry."
+        
+        # check that we have define the input
+        if len(args) == 0:
+            if self.run_name:
+                args[0] = self.run_name
+            else:
+                raise self.InvalidCmd, 'no default run. Please specify the run_name'
+                
+        # always pass to a path + get the event size
+        result_file= sys.stdout
+        if not os.path.sep in args[0]:
+            path = [pjoin(self.me_dir, 'Events', args[0], 'unweighted_events.lhe.gz'),
+                    pjoin(self.me_dir, 'Events', args[0], 'unweighted_events.lhe'),
+                    pjoin(self.me_dir, 'Events', args[0], 'events.lhe.gz'),
+                    pjoin(self.me_dir, 'Events', args[0], 'events.lhe')]
+            
+            for p in path:
+                if os.path.exists(p):
+                    nb_event = self.results[args[0]].get_current_info()['nb_event']
+                    
+                    
+                    if self.run_name != args[0]:
+                        tag = self.results[args[0]].tags[0]
+                        self.set_run_name(args[0], tag,'parton', False)
+                    result_file = open(pjoin(self.me_dir,'Events', self.run_name, 'parton_systematics.log'),'w')
+                    args[0] = p
+                    break
+            else:
+                raise self.InvalidCmd, 'Invalid run name. Please retry'
+        elif self.options['nb_core'] != 1:
+            lhe = lhe_parser.EventFile(args)
+            nb_event = len(lhe)
+            lhe.close()
+
+        input = args[0]
+        if len(args)>1:
+            output = pjoin(os.getcwd(),args[1])
+        else:
+            output = input
+    
+        lhaid = [self.run_card.get_lhapdf_id()]
+        if 'store_rwgt_info' in self.run_card and not self.run_card['store_rwgt_info']:
+            raise self.InvalidCmd,  "The events was not generated with store_rwgt_info=True. Can not evaluate systematics error on this event file."
+        elif 'use_syst'  in self.run_card and not self.run_card['use_syst']:
+            raise self.InvalidCmd,  "The events was not generated with use_syst=True. Can not evaluate systematics error on this event file."
+        
+        try:
+            pdfsets_dir = self.get_lhapdf_pdfsetsdir()
+        except Exception, error:
+            logger.debug(str(error))
+            logger.warning('Systematic computation requires lhapdf to run. Bypass Systematics')
+            return
+
+        if '--from_card' in opts:
+            opts.remove('--from_card')
+            opts.append('--from_card=internal')
+            
+            # Check that all pdfset are correctly installed
+            if 'sys_pdf' in self.run_card:
+                sys_pdf = self.run_card['sys_pdf'].split('&&')
+                lhaid += [l.split()[0] for l in sys_pdf]
+
+        else:
+            #check that all p
+            pdf = [a[6:] for a in opts if a.startswith('--pdf=')]
+            lhaid += [t.split('@')[0] for p in pdf for t in p.split(',') 
+                                            if t not in ['errorset', 'central']]
+        
+        # Copy all the relevant PDF sets
+        [self.copy_lhapdf_set([onelha], pdfsets_dir) for onelha in lhaid]
+            
+        
+        if self.options['run_mode'] ==2:
+            nb_submit = min(self.options['nb_core'], nb_event//2500)
+        elif self.options['run_mode'] ==1:
+            nb_submit = min(self.options['cluster_size'], nb_event//25000)
+        else:
+            nb_submit =1 
+
+        if MADEVENT:
+            import internal.systematics as systematics
+        else:
+            import madgraph.various.systematics as systematics
+
+        #one core:
+        if nb_submit in [0,1]:
+            systematics.call_systematics([input, output] + opts, 
+                                         log=lambda x: logger.info(str(x)),
+                                         result=result_file
+                                         )
+            
+        elif self.options['run_mode'] in [1,2]:
+            event_per_job = nb_event // nb_submit
+            nb_job_with_plus_one = nb_event % nb_submit
+            start_event, stop_event = 0,0
+            for i in range(nb_submit):
+                #computing start/stop event
+                event_requested = event_per_job
+                if i < nb_job_with_plus_one:
+                    event_requested += 1
+                start_event = stop_event
+                stop_event = start_event + event_requested
+                    
+                prog = sys.executable
+                input_files = [os.path.basename(input)]
+                output_files = ['./tmp_%s_%s' % (i, os.path.basename(output)),
+                                './log_sys_%s.txt' % (i)]
+                argument = []
+                if not __debug__:
+                    argument.append('-O')
+                argument +=  [pjoin(self.me_dir, 'bin', 'internal', 'systematics.py'),
+                             input_files[0], output_files[0]] + opts +\
+                             ['--start_event=%i' % start_event,
+                              '--stop_event=%i' %stop_event,
+                              '--result=./log_sys_%s.txt' %i,
+                              '--lhapdf_config=%s' % self.options['lhapdf']]
+                required_output = output_files            
+                self.cluster.cluster_submit(prog, argument, 
+                                            input_files=input_files,
+                                            output_files=output_files,
+                                            cwd=os.path.dirname(output),
+                                            required_output=required_output,
+                                            stdout='/dev/null'
+                                            )
+            starttime = time.time()
+            update_status = lambda idle, run, finish: \
+                    self.update_status((idle, run, finish, 'running systematics'), level=None,
+                                       force=False, starttime=starttime)
+
+            self.cluster.wait(os.path.dirname(output), update_status, update_first=update_status)
+            
+            #collect the data
+            all_cross = []
+            for i in range(nb_submit):
+                pos=0
+                for line in open(pjoin(os.path.dirname(output), 'log_sys_%s.txt'%i)):
+                    if line.startswith('#'):
+                        continue
+                    split = line.split()
+                    if len(split) in [0,1]:
+                        continue
+                    key = tuple(float(x) for x in split[:-1])
+                    cross= float(split[-1])
+                    if 'event_norm' in self.run_card and \
+                            self.run_card['event_norm'] in ['average', 'unity']:
+                        cross *= (event_per_job+1 if i <nb_job_with_plus_one else event_per_job)
+                    if len(all_cross) > pos:
+                        all_cross[pos] += cross
+                    else:
+                        all_cross.append(cross)
+                    pos+=1
+                        
+            if 'event_norm' in self.run_card and \
+                                       self.run_card['event_norm'] in ['unity']:
+                all_cross= [cross/nb_event for cross in all_cross]
+                
+            sys_obj = systematics.call_systematics([input, None] + opts, 
+                                         log=lambda x: logger.info(str(x)),
+                                         result=result_file,
+                                         running=False
+                                         )                    
+            sys_obj.print_cross_sections(all_cross, nb_event, result_file)
+            
+            #concatenate the output file
+            subprocess.call(['cat']+\
+                            ['./tmp_%s_%s' % (i, os.path.basename(output)) for i in range(nb_submit)],
+                            stdout=open(output,'w'),
+                            cwd=os.path.dirname(output))
+            for i in range(nb_submit):
+                os.remove('%s/tmp_%s_%s' %(os.path.dirname(output),i,os.path.basename(output)))
+            #    os.remove('%s/log_sys_%s.txt' % (os.path.dirname(output),i))
+                                                  
+
+            
+            
+
+        self.update_status('End of systematics computation', level='parton', makehtml=False)
+        
+        
     ############################################################################
     def do_reweight(self, line):
         """ syntax is "reweight RUN_NAME"
@@ -1096,9 +1458,50 @@
             cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Reweight
         """
         
+        #### Utility function
+        def check_multicore(self):
+            """ determine if the cards are save for multicore use"""
+            card = pjoin(self.me_dir, 'Cards', 'reweight_card.dat')
+
+            multicore = True
+            if self.options['run_mode'] in [0,1]:
+                multicore = False
+
+            lines = [l.strip() for l in open(card) if not l.strip().startswith('#')]
+            while lines and not lines[0].startswith('launch'):
+                line = lines.pop(0)
+                # if not standard output mode forbid multicore mode
+                if line.startswith('change') and line[6:].strip().startswith('output'):
+                    return False
+                if line.startswith('change') and line[6:].strip().startswith('multicore'):
+                    split_line = line.split()
+                    if len(split_line) > 2: 
+                        multicore = bool(split_line[2])
+            # we have reached the first launch in the card ensure that no output change 
+            #are done after that point.
+            lines = [line[6:].strip() for line in lines if line.startswith('change')]
+            for line in lines:
+                if line.startswith(('process','model','output', 'rwgt_dir')):
+                    return False
+                elif line.startswith('multicore'):
+                    split_line = line.split()
+                    if len(split_line) > 1: 
+                        multicore = bool(split_line[1])
+
+            return multicore
+            
+        
+        
         if '-from_cards' in line and not os.path.exists(pjoin(self.me_dir, 'Cards', 'reweight_card.dat')):
             return
-        
+        # option for multicore to avoid that all of them create the same directory
+        if '--multicore=create' in line:
+            multicore='create'
+        elif '--multicore=wait' in line:
+            multicore='wait'
+        else:
+            multicore=False
+            
         # Check that MG5 directory is present .
         if MADEVENT and not self.options['mg5_path']:
             raise self.InvalidCmd, '''The module reweight requires that MG5 is installed on the system.
@@ -1125,6 +1528,10 @@
             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()
+            # ensure that the run_card is present
+            if not hasattr(self, 'run_card'):
+                self.run_card = banner_mod.RunCard(pjoin(self.me_dir, 'Cards', 'run_card.dat'))
+            
             # 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')):
@@ -1134,41 +1541,113 @@
             if not isinstance(self, cmd.CmdShell):
                 command.append('--web')
             command.append('reweight')
-            if self.run_name:
-                command.append(self.run_name)
+            
+            #########   START SINGLE CORE MODE ############
+            if self.options['nb_core']==1 or self.run_card['nevents'] < 101 or not check_multicore(self):
+                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=os.getcwd())
+                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
+                ##########    END SINGLE CORE HANDLING #############
             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")
+                ##########    START MULTI-CORE HANDLING #############
+                if not isinstance(self.cluster, cluster.MultiCore):
+                    mycluster = cluster.MultiCore(nb_core=self.options['nb_core'])
+                else:
+                    mycluster = self.cluster
+                
+                new_args=list(args)
+                self.check_decay_events(new_args) 
+                try:
+                    os.remove(pjoin(self.me_dir,'rw_me','rwgt.pkl'))
+                except Exception, error:
+                    pass
+                # prepare multi-core  job:
+                import madgraph.various.lhe_parser as lhe_parser
+                # args now alway content the path to the valid files
+                if 'nevt_job' in self.run_card and self.run_card['nevt_job'] !=-1:
+                    nevt_job = self.run_card['nevt_job']
+                else:
+                    nevt_job = max(5000, self.run_card['nevents']/50)
+                logger.info("split the event file in bunch of %s events" % nevt_job)
+                nb_file = lhe_parser.EventFile(new_args[0]).split(nevt_job)
+                starttime = time.time()
+                update_status = lambda idle, run, finish: \
+                    self.update_status((idle, run, finish, 'reweight'), level=None,
+                                       force=False, starttime=starttime)
+
+                all_lhe = []
+                devnull= open(os.devnull)
+                for i in range(nb_file):
+                    new_command = list(command) 
+                    new_command.append('%s_%s.lhe' % (new_args[0],i))
+                    all_lhe.append('%s_%s.lhe' % (new_args[0],i))
+                    if '-from_cards' not in command:
+                        new_command.append('-from_cards')
+                    if i==0:
+                        if __debug__:
+                            stdout = None
+                        else:
+                            stdout = open(pjoin(self.me_dir,'Events', self.run_name, 'reweight.log'),'w')
+                        new_command.append('--multicore=create')
+                    else:
+                        stdout = devnull
+                        #stdout = open(pjoin(self.me_dir,'Events', self.run_name, 'reweight%s.log' % i),'w')
+                        new_command.append('--multicore=wait')
+                    mycluster.submit(prog=command[0], argument=new_command[1:], stdout=stdout, cwd=os.getcwd())
+                mycluster.wait(self.me_dir,update_status)
+                devnull.close()
+                
+                lhe = lhe_parser.MultiEventFile(all_lhe, parse=False)
+                nb_event, cross_sections = lhe.write(new_args[0], get_info=True)
+                if any(os.path.exists('%s_%s_debug.log' % (f, self.run_tag)) for f in all_lhe):
+                    for f in all_lhe:
+                        if os.path.exists('%s_%s_debug.log' % (f, self.run_tag)):
+                            raise Exception, "Some of the run failed: Please read %s_%s_debug.log" % (f, self.run_tag) 
+                
+                
+                if 'event_norm' in self.run_card and self.run_card['event_norm'] == 'average':
+                    for key, value in cross_sections.items():
+                        cross_sections[key] = value / (nb_event+1)
+                lhe.remove()
+                for key in cross_sections:
+                    if key == 'orig' or key.isdigit():
+                        continue
+                    logger.info('%s : %s pb' % (key, cross_sections[key]))
                 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              
+            ##########    END MULTI-CORE HANDLING #############
+                              
 
         self.to_store.append('event')
         # forbid this function to create an empty item in results.
@@ -1189,6 +1668,8 @@
         path = pjoin(self.me_dir, 'Cards', 'reweight_card.dat')
         reweight_cmd.raw_input=False
         reweight_cmd.me_dir = self.me_dir
+        reweight_cmd.multicore = multicore #allow the directory creation or not
+        print "We are in mode", multicore
         reweight_cmd.import_command_file(path)
         reweight_cmd.do_quit('')
             
@@ -2319,6 +2800,8 @@
 
         pdfsetname=set()
         for lhaid in lhaid_list:
+            if  isinstance(lhaid, str) and lhaid.isdigit():
+                lhaid = int(lhaid)
             if isinstance(lhaid, (int,float)):
                 try:
                     if lhaid in self.lhapdf_pdfsets:
 
=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2016-08-17 19:36:40 +0000
+++ madgraph/interface/madevent_interface.py	2016-08-24 22:15:21 +0000
@@ -2012,11 +2012,17 @@
             if not bypass_run:
                 self.exec_cmd('refine %s' % nb_event, postcmd=False)
             
-                self.exec_cmd('combine_events', postcmd=False)
+                self.exec_cmd('combine_events', postcmd=False,printcmd=False)
                 self.print_results_in_shell(self.results.current)
 
-            
-                self.run_syscalc('parton')
+                if self.run_card['use_syst']:
+                    scdir = self.options['syscalc_path']
+                    if not scdir or not os.path.exists(scdir):
+                        self.exec_cmd('systematics %s --from_card' % self.run_name, postcmd=False,printcmd=False)    
+                    else:
+                        self.run_syscalc('parton')
+                
+                    
                 self.create_plot('parton')            
                 self.exec_cmd('store_events', postcmd=False)            
                 self.exec_cmd('reweight -from_cards', postcmd=False)            
@@ -3972,106 +3978,6 @@
         self.Gdirs = Gdirs
         return self.Gdirs
                 
-    ############################################################################
-    def set_run_name(self, name, tag=None, level='parton', reload_card=False,
-                     allow_new_tag=True):
-        """define the run name, the run_tag, the banner and the results."""
-        
-        # when are we force to change the tag new_run:previous run requiring changes
-        upgrade_tag = {'parton': ['parton','pythia','pgs','delphes'],
-                       'pythia': ['pythia','pgs','delphes'],
-                       'pgs': ['pgs'],
-                       'delphes':['delphes'],
-                       'plot':[],
-                       'syscalc':[]}
-        
-        
-
-        if name == self.run_name:        
-            if reload_card:
-                run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
-                self.run_card = banner_mod.RunCard(run_card)
-
-            #check if we need to change the tag
-            if tag:
-                self.run_card['run_tag'] = tag
-                self.run_tag = tag
-                self.results.add_run(self.run_name, self.run_card)
-            else:
-                for tag in upgrade_tag[level]:
-                    if getattr(self.results[self.run_name][-1], tag):
-                        tag = self.get_available_tag()
-                        self.run_card['run_tag'] = tag
-                        self.run_tag = tag
-                        self.results.add_run(self.run_name, self.run_card)                        
-                        break
-            return # Nothing to do anymore
-        
-        # save/clean previous run
-        if self.run_name:
-            self.store_result()
-        # store new name
-        self.run_name = name
-        
-        new_tag = False
-        # First call for this run -> set the banner
-        self.banner = banner_mod.recover_banner(self.results, level, name)
-        if 'mgruncard' in self.banner:
-            self.run_card = self.banner.charge_card('run_card')
-        else:
-            # Read run_card
-            run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
-            self.run_card = banner_mod.RunCard(run_card)   
-        
-        if tag:
-            self.run_card['run_tag'] = tag
-            new_tag = True
-        elif not self.run_name in self.results and level =='parton':
-            pass # No results yet, so current tag is fine
-        elif not self.run_name in self.results:
-            #This is only for case when you want to trick the interface
-            logger.warning('Trying to run data on unknown run.')
-            self.results.add_run(name, self.run_card)
-            self.results.update('add run %s' % name, 'all', makehtml=False)
-        else:
-            for tag in upgrade_tag[level]:
-                
-                if getattr(self.results[self.run_name][-1], tag):
-                    # LEVEL is already define in the last tag -> need to switch tag
-                    tag = self.get_available_tag()
-                    self.run_card['run_tag'] = tag
-                    new_tag = True
-                    break
-            if not new_tag:
-                # We can add the results to the current run
-                tag = self.results[self.run_name][-1]['tag']
-                self.run_card['run_tag'] = tag # ensure that run_tag is correct                
-                   
-        if allow_new_tag and (name in self.results and not new_tag):
-            self.results.def_current(self.run_name)
-        else:
-            self.results.add_run(self.run_name, self.run_card)
-
-        self.run_tag = self.run_card['run_tag']
-
-        # Return the tag of the previous run having the required data for this
-        # tag/run to working wel.
-        if level == 'parton':
-            return
-        elif level == 'pythia':
-            return self.results[self.run_name][0]['tag']
-        else:
-            for i in range(-1,-len(self.results[self.run_name])-1,-1):
-                tagRun = self.results[self.run_name][i]
-                if tagRun.pythia:
-                    return tagRun['tag']
-            
-            
-        
-        
-        
-        
-        
 
     ############################################################################
     def find_model_name(self):
@@ -4215,7 +4121,7 @@
         scdir = self.options['syscalc_path']
         if not scdir or not os.path.exists(scdir):
             return
-        logger.info('running syscalc on mode %s' % mode)    
+        logger.info('running SysCalc on mode %s' % mode)    
     
 
         # Check that all pdfset are correctly installed
@@ -4265,7 +4171,9 @@
                 if not (os.path.exists(event_path) or os.path.exists(event_path+".gz")):
                     event_path = pjoin(event_dir, 'unweighted_events.lhe')
                 output = pjoin(event_dir, 'syscalc.lhe')
+                stdout = open(pjoin(event_dir, self.run_name, '%s_systematics.log' % (mode)),'w')
             elif mode == 'Pythia':
+                stdout = open(pjoin(event_dir, self.run_name, '%s_%s_syscalc.log' % (tag,mode)),'w')
                 if 'mgpythiacard' in self.banner:
                     pat = re.compile('''^\s*qcut\s*=\s*([\+\-\d.e]*)''', re.M+re.I)
                     data = pat.search(self.banner['mgpythiacard'])
@@ -4294,7 +4202,7 @@
         try:
             proc = misc.call([os.path.join(scdir, 'sys_calc'),
                                event_path, card, output],
-                            stdout = open(pjoin(event_dir, self.run_name, '%s_%s_syscalc.log' % (tag,mode)),'w'),
+                            stdout = stdout,
                             stderr = subprocess.STDOUT,
                             cwd=event_dir)
             # Wait 5 s to make sure file is finished writing
 
=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2016-08-22 20:43:28 +0000
+++ madgraph/interface/madgraph_interface.py	2016-08-24 22:15:21 +0000
@@ -192,9 +192,9 @@
                             info['date'])
 
         if os.path.exists(pjoin(MG5DIR, '.bzr')):
-            proc = subprocess.Popen(['bzr', 'nick'], stdout=subprocess.PIPE)
+            proc = subprocess.Popen(['bzr', 'nick'], stdout=subprocess.PIPE,cwd=MG5DIR)
             bzrname,_ = proc.communicate()
-            proc = subprocess.Popen(['bzr', 'revno'], stdout=subprocess.PIPE)
+            proc = subprocess.Popen(['bzr', 'revno'], stdout=subprocess.PIPE,cwd=MG5DIR)
             bzrversion,_ = proc.communicate() 
             bzrname, bzrversion = bzrname.strip(), bzrversion.strip() 
             len_name = len(bzrname)
 
=== modified file 'madgraph/interface/reweight_interface.py'
--- madgraph/interface/reweight_interface.py	2016-08-18 22:31:05 +0000
+++ madgraph/interface/reweight_interface.py	2016-08-24 22:15:21 +0000
@@ -82,11 +82,11 @@
         self.model = None
         self.has_standalone_dir = False
         self.mother= mother # calling interface
-            
+        self.multicore=False
         
         self.options = {'curr_dir': os.path.realpath(os.getcwd()),
                         'rwgt_name':None}
-        
+
         self.events_file = None
         self.processes = {}
         self.second_model = None
@@ -201,7 +201,7 @@
         logger.info("options: %s" % option)
 
 
-    def get_LO_definition_from_NLO(self, proc):
+    def get_LO_definition_from_NLO(self, proc, real_only=False):
         """return the LO definitions of the process corresponding to the born/real"""
         
         # split the line definition with the part before and after the NLO tag
@@ -212,7 +212,10 @@
             #NO NLO tag => nothing to do actually return input
             return proc
         elif not order.startswith(('virt','loonly','noborn')):
-            # OK this a standard NLO process
+            # OK this a standard NLO process            
+            if real_only:
+                commandline= '' 
+            
             if '=' in order:
                 # get the type NLO QCD/QED/...
                 order = order.split('=',1)[1]
@@ -358,6 +361,10 @@
             self.rwgt_dir = args[1]
             if not os.path.exists(self.rwgt_dir):
                 os.mkdir(self.rwgt_dir)
+        elif args[0] == 'multicore':
+            pass 
+            # this line is meant to be parsed by common_run_interface and change the way this class is called.
+            #It has no direct impact on this class.
         else:
             logger.critical("unknown option! %s.  Discard line." % args[0])
         
@@ -415,12 +422,25 @@
 
         model_line = self.banner.get('proc_card', 'full_model_line')
 
-        if not self.has_standalone_dir:
+        if not self.has_standalone_dir:                
             if self.rwgt_dir and os.path.exists(pjoin(self.rwgt_dir,'rw_me','rwgt.pkl')):
                 self.load_from_pickle()
-                self.me_dir = self.rwgt_dir
+                if not self.rwgt_dir:
+                    self.me_dir = self.rwgt_dir
+            elif self.multicore == 'wait':
+                while not os.path.exists(pjoin(self.me_dir,'rw_me','rwgt.pkl')):
+                    time.sleep(60)
+                    print 'wait for pickle'                  
+                print "loading from pickle"
+                if not self.rwgt_dir:
+                    self.rwgt_dir = self.me_dir
+                self.load_from_pickle(keep_name=True)
             else:
-                self.create_standalone_directory()        
+                self.create_standalone_directory()  
+                if self.multicore == 'create':
+                    if not self.rwgt_dir:
+                        self.rwgt_dir = self.me_dir
+                    self.save_to_pickle()      
         
         if self.rwgt_dir:
             path_me =self.rwgt_dir
@@ -501,7 +521,6 @@
                 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 \
@@ -779,10 +798,11 @@
                 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)' %\
+        if self.multicore != 'create':
+            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)
@@ -1023,6 +1043,11 @@
     def calculate_matrix_element(self, event, hypp_id, space, scale2=0):
         """routine to return the matrix element"""
         
+        if self.has_nlo:
+            nb_retry, sleep = 10, 60 
+        else:
+            nb_retry, sleep = 5, 20 
+        
         tag, order = event.get_tag_and_order()
         if isinstance(hypp_id, str) and hypp_id.startswith('V'):
             tag = (tag,'V')
@@ -1068,7 +1093,7 @@
                         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)                
+                    misc.multiple_try(nb_retry,sleep)(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)
@@ -1096,10 +1121,9 @@
                 try:
                     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)
+                    misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix%spy.so' % newtag], cwd=Pdir)
                     mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, newtag), globals(), locals(), [],-1)
                  
                 S = mymod.SubProcesses
@@ -1121,7 +1145,7 @@
             Pname = os.path.basename(Pdir)
             os.environ['MENUM'] = '2'
             if not self.rwgt_dir or not os.path.exists(pjoin(Pdir, 'matrix2py.so')):
-                misc.compile(['matrix2py.so'], cwd=pjoin(subdir, Pdir))
+                misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix2py.so'], cwd=pjoin(subdir, Pdir))
             if (Pdir, 1) not in dir_to_f2py_free_mod:
                 metag = 1
                 dir_to_f2py_free_mod[(Pdir,1)] = (metag, nb_f2py_module)
@@ -1134,10 +1158,9 @@
             try:
                 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
                 os.environ['MENUM'] = str(metag)
-                misc.compile(['matrix%spy.so' % metag], cwd=pjoin(subdir, Pdir))
+                misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix%spy.so' % metag], cwd=pjoin(subdir, Pdir))
                 mymod = __import__("%s_second.SubProcesses.%s.matrix%spy" % (base, Pname, metag))
                      
             reload(mymod)
@@ -1211,21 +1234,24 @@
                 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:
-            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] ))  
+                    
+        if not self.multicore == 'create':
+            # No print of results for the multicore mode for the one printed on screen
+            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:
+                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:
+        if self.rwgt_dir and self.multicore == False:
             self.save_to_pickle()
         
         with misc.stdchannel_redirected(sys.stdout, os.devnull):
@@ -1245,7 +1271,7 @@
     @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']
@@ -1350,12 +1376,18 @@
             logger.info('generating the square matrix element for reweighting (second model and/or processes)')
         start = time.time()
         commandline=''
-        for proc in data['processes']:
+        for i,proc in enumerate(data['processes']):
             if '[' not in proc:
                 commandline += "add process %s ;" % proc
             else:
                 has_nlo = True
-                commandline += self.get_LO_definition_from_NLO(proc)
+                if self.banner.get('run_card','ickkw') == 3:
+                    if len(proc) == min([len(p.strip()) for p in data['processes']]):
+                        commandline += self.get_LO_definition_from_NLO(proc)
+                    else:
+                        commandline += self.get_LO_definition_from_NLO(proc, real_only=True)
+                else:
+                    commandline += self.get_LO_definition_from_NLO(proc)
         
         commandline = commandline.replace('add process', 'generate',1)
         logger.info(commandline)
@@ -1434,6 +1466,18 @@
                          pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat'), 
                          commentdefault=False)
             
+            if self.multicore == 'create':
+                print "compile OLP", data['paths'][0]
+                misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][0],'SubProcesses'),
+                             nb_core=self.mother.options['nb_core'])
+        
+        if os.path.exists(pjoin(path_me, data['paths'][1], 'Cards', 'MadLoopParams.dat')):
+            if self.multicore == 'create':
+                print "compile OLP", data['paths'][1]
+                misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'),
+                             nb_core=self.mother.options['nb_core'])
+                
+        
         # 5. create the virtual for NLO reweighting  ---------------------------
         if has_nlo and 'NLO' in self.rwgt_mode:
             # Do not pass here for LO/NLO_tree
@@ -1527,6 +1571,11 @@
                                   self.banner.run_card['lpp2']],
                                  self.banner.run_card.get_lhapdf_id())
                 self.combine_wgt = mymod.get_wgt
+
+            if self.multicore == 'create':
+                print "compile OLP", data['paths'][1]
+                misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'),
+                             nb_core=self.mother.options['nb_core'])
                 
         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
@@ -1630,18 +1679,22 @@
         to_save['rwgt_dir'] = self.rwgt_dir
         to_save['has_nlo'] = self.has_nlo
         to_save['rwgt_mode'] = self.rwgt_mode
+        to_save['rwgt_name'] = self.options['rwgt_name']
 
         name = pjoin(self.rwgt_dir, 'rw_me', 'rwgt.pkl')
         save_load_object.save_to_file(name, to_save)
 
-    def load_from_pickle(self):
+
+    def load_from_pickle(self, keep_name=False):
         import madgraph.iolibs.save_load_object as save_load_object
         
         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()),
-                        'rewgt_name': None}
+                        'rwgt_name': None}
+        if keep_name:
+            self.options['rwgt_name'] = obj['rwgt_name']
         
         old_rwgt = obj['rwgt_dir']
            
@@ -1665,10 +1718,14 @@
         if not self.rwgt_mode:
             self.rwgt_mode = obj['rwgt_mode']
             logger.info("mode set to %s" % self.rwgt_mode)
-        if self.has_nlo:
+        if self.has_nlo and 'NLO' in self.rwgt_mode:
             path = pjoin(obj['rwgt_dir'], 'rw_mevirt','Source')
             sys.path.insert(0, path)
-            mymod = __import__('rwgt2py', globals(), locals())
+            try:
+                mymod = __import__('rwgt2py', globals(), locals())
+            except ImportError:
+                misc.compile(['rwgt2py.so'], cwd=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']],
 
=== modified file 'madgraph/iolibs/export_fks.py'
--- madgraph/iolibs/export_fks.py	2016-08-07 07:16:07 +0000
+++ madgraph/iolibs/export_fks.py	2016-08-24 22:15:21 +0000
@@ -268,7 +268,8 @@
           pjoin('various','FO_analyse_card.py'),
           pjoin('various','histograms.py'),      
           pjoin('various','banner.py'),          
-          pjoin('various','cluster.py'),          
+          pjoin('various','cluster.py'),
+          pjoin('various','systematics.py'),          
           pjoin('various','lhe_parser.py'),
           pjoin('madevent','sum_html.py'),
           pjoin('madevent','gen_crossxhtml.py'),          
 
=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2016-08-15 14:06:30 +0000
+++ madgraph/iolibs/export_v4.py	2016-08-24 22:15:21 +0000
@@ -3227,6 +3227,8 @@
                                 self.dir_path+'/bin/internal/lhe_parser.py')                        
         cp(_file_path+'/various/banner.py', 
                                    self.dir_path+'/bin/internal/banner.py')
+        cp(_file_path+'/various/systematics.py', self.dir_path+'/bin/internal/systematics.py')        
+
         cp(_file_path+'/various/cluster.py', 
                                        self.dir_path+'/bin/internal/cluster.py') 
         cp(_file_path+'/madevent/combine_runs.py', 
 
=== modified file 'madgraph/madevent/gen_crossxhtml.py'
--- madgraph/madevent/gen_crossxhtml.py	2016-04-15 11:14:25 +0000
+++ madgraph/madevent/gen_crossxhtml.py	2016-08-24 22:15:21 +0000
@@ -812,7 +812,7 @@
                 self.parton.append('param_card')
             
             if 'syst' not in self.parton and \
-                                    exists(pjoin(path, "%s_parton_syscalc.log" %self['tag'])):
+                                    exists(pjoin(path, "parton_systematics.log")):
                 self.parton.append('syst')
 
             for kind in ['top','HwU','pdf','ps']:
@@ -1150,8 +1150,8 @@
                     local_dico['err'] = self['error']
                     local_dico['nb_event'] = self['nb_event']
                     if 'syst' in self.parton:
-                        local_dico['syst'] = '<font face=symbol>&#177;</font> <a href="./Events/%(run_name)s/%(tag)s_parton_syscalc.log">systematics</a>' \
-                                             % {'run_name':self['run_name'], 'tag': self['tag']}
+                        local_dico['syst'] = '<font face=symbol>&#177;</font> <a href="./Events/%(run_name)s/parton_systematics.log">systematics</a>' \
+                                             % {'run_name':self['run_name']}
                 elif self['cross_pythia']:
                     if self.parton:
                         local_dico['cross_span'] = nb_line -1
@@ -1173,7 +1173,7 @@
                     local_dico['err'] = self['error']
                     local_dico['nb_event'] = self['nb_event']
                     if 'syst' in self.parton:
-                        local_dico['syst'] = '<font face=symbol>&#177;</font> <a href="./Events/%(run_name)s/%(tag)s_parton_syscalc.log">systematics</a>' \
+                        local_dico['syst'] = '<font face=symbol>&#177;</font> <a href="./Events/%(run_name)s/parton_systematics.log">systematics</a>' \
                                              % {'run_name':self['run_name'], 'tag': self['tag']}
                     
             elif ttype == 'pythia' and self['cross_pythia']:
 
=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2016-08-17 19:36:40 +0000
+++ madgraph/various/banner.py	2016-08-24 22:15:21 +0000
@@ -193,34 +193,6 @@
             return cross, math.sqrt(error)
         
 
-    
-    def modify_init_cross(self, cross):
-        """modify the init information with the associate cross-section"""
-
-        assert isinstance(cross, dict)
-#        assert "all" in cross
-        assert "init" in self
-        
-        all_lines = self["init"].split('\n')
-        new_data = []
-        new_data.append(all_lines[0])
-        for i in range(1, len(all_lines)):
-            line = all_lines[i]
-            split = line.split()
-            if len(split) == 4:
-                xsec, xerr, xmax, pid = split 
-            else:
-                new_data += all_lines[i:]
-                break
-            if int(pid) not in cross:
-                raise Exception
-            pid = int(pid)
-            ratio = cross[pid]/float(xsec)
-            line = "   %+13.7e %+13.7e %+13.7e %i" % \
-                (float(cross[pid]), ratio* float(xerr), ratio*float(xmax), pid)
-            new_data.append(line)
-        self['init'] = '\n'.join(new_data)
-    
     def scale_init_cross(self, ratio):
         """modify the init information with the associate scale"""
 
@@ -244,6 +216,15 @@
             new_data.append(line)
         self['init'] = '\n'.join(new_data)
     
+    def get_pdg_beam(self):
+        """return the pdg of each beam"""
+        
+        assert "init" in self
+        
+        all_lines = self["init"].split('\n')
+        pdg1,pdg2,_ = all_lines[0].split(None, 2)
+        return int(pdg1), int(pdg2)
+    
     def load_basic(self, medir):
         """ Load the proc_card /param_card and run_card """
         
 
=== modified file 'madgraph/various/cluster.py'
--- madgraph/various/cluster.py	2016-06-13 08:27:43 +0000
+++ madgraph/various/cluster.py	2016-08-24 22:15:21 +0000
@@ -618,6 +618,8 @@
                     if isinstance(exe,str):
                         if os.path.exists(exe) and not exe.startswith('/'):
                             exe = './' + exe
+                        if isinstance(opt['stdout'],str):
+                            opt['stdout'] = open(opt['stdout'],'w')
                         if opt['stderr'] == None:
                             opt['stderr'] = subprocess.STDOUT
                         proc = misc.Popen([exe] + arg,  **opt)
@@ -673,7 +675,6 @@
         
         tag = (prog, tuple(argument), cwd, nb_submit)
         if isinstance(prog, str):
-            
     
             opt = {'cwd': cwd, 
                    'stdout':stdout,
 
=== modified file 'madgraph/various/lhe_parser.py'
--- madgraph/various/lhe_parser.py	2016-07-29 15:27:30 +0000
+++ madgraph/various/lhe_parser.py	2016-08-24 22:15:21 +0000
@@ -178,6 +178,8 @@
             if '.gz' in path and isinstance(self, EventFileNoGzip) and\
                 mode == 'r' and os.path.exists(path[:-3]):
                 super(EventFile, self).__init__(path[:-3], mode, *args, **opt)
+            else:
+                raise
                 
         self.banner = ''
         if mode == 'r':
@@ -506,7 +508,26 @@
         else:
             return out
 
+    def split(self, nb_event=0):
+        """split the file in multiple file. Do not change the weight!"""
     
+        nb_file = -1
+        for i, event in enumerate(self):
+            if i % nb_event == 0:
+                if i:
+                    #close previous file
+                    current.write('</LesHouchesEvent>\n')
+                    current.close()
+                # create the new file
+                nb_file +=1
+                current = open('%s_%s.lhe' % (self.name, nb_file),'w')
+                current.write(self.banner)
+            current.write(str(event))
+        if i!=0:
+            current.write('</LesHouchesEvent>\n')
+            current.close()   
+        return nb_file +1
+
     def update_Hwu(self, hwu, fct, name='lhe', keep_wgt=True):
         
         first=True
@@ -620,13 +641,14 @@
        The number of events in each file need to be provide in advance 
        (if not provide the file is first read to find that number"""
     
-    def __new__(cls, start_list=[]):
+    def __new__(cls, start_list=[],parse=True):
         return object.__new__(MultiEventFile)
     
-    def __init__(self, start_list=[]):
+    def __init__(self, start_list=[], parse=True):
         """if trunc_error is define here then this allow
         to only read all the files twice and not three times."""
         self.files = []
+        self.parsefile = parse #if self.files is formatted or just the path
         self.banner = ''
         self.initial_nb_events = []
         self.total_event_in_files = 0
@@ -636,8 +658,11 @@
         self.across = []
         self.scales = []
         if start_list:
-            for p in start_list:
-                self.add(p)
+            if parse:
+                for p in start_list:
+                    self.add(p)
+            else:
+                self.files = start_list
         self._configure = False
         
     def add(self, path, cross, error, across):
@@ -881,6 +906,66 @@
 
         return super(MultiEventFile, self).unweight(outputpath, get_wgt_multi, **opts)
 
+    def write(self, path, random=False, banner=None, get_info=False):
+        """ """
+        
+        if isinstance(path, str):
+            out = EventFile(path, 'w')
+            if self.parsefile and not banner:    
+                banner = self.files[0].banner
+            elif not banner:
+                firstlhe = EventFile(self.files[0])
+                banner = firstlhe.banner                
+        else: 
+            out = path
+        if banner:
+            out.write(banner)
+        nb_event = 0
+        info = collections.defaultdict(float)
+        if random and self.open:
+            for event in self:
+                nb_event +=1
+                out.write(event)
+                if get_info:
+                    event.parse_reweight()
+                    for key, value in event.reweight_data.items():
+                        info[key] += value
+                    info['central'] += event.wgt
+        elif not random:
+            for i,f in enumerate(self.files):
+                #check if we need to parse the file or not
+                if not self.parsefile:
+                    if i==0:
+                        try:
+                            lhe = firstlhe
+                        except:
+                            lhe = EventFile(f)
+                    else:
+                        lhe = EventFile(f)
+                else:
+                    lhe = f
+                for event in lhe:
+                    nb_event +=1
+                    if get_info:
+                        event.parse_reweight()
+                        for key, value in event.reweight_data.items():
+                            info[key] += value
+                        info['central'] += event.wgt
+                    out.write(str(event))
+                lhe.close()
+        out.write("</LesHouchesEvents>\n") 
+        return nb_event, info
+                            
+    def remove(self):
+        """ """
+        if self.parsefile:
+            for f in self.files:
+                os.remove(f.name)
+        else:
+            for f in self.files:
+                os.remove(f)
+            
+        
            
 class Event(list):
     """Class storing a single event information (list of particles + global information)"""
@@ -1025,7 +1110,7 @@
         self.reweight_order = []
         start, stop = self.tag.find('<rwgt>'), self.tag.find('</rwgt>')
         if start != -1 != stop :
-            pattern = re.compile(r'''<\s*wgt id=(?:\'|\")(?P<id>[^\'\"]+)(?:\'|\")\s*>\s*(?P<val>[\ded+-.]*)\s*</wgt>''')
+            pattern = re.compile(r'''<\s*wgt id=(?:\'|\")(?P<id>[^\'\"]+)(?:\'|\")\s*>\s*(?P<val>[\ded+-.]*)\s*</wgt>''',re.I)
             data = pattern.findall(self.tag)
             try:
                 self.reweight_data = dict([(pid, float(value)) for (pid, value) in data
@@ -1046,7 +1131,54 @@
         
             text = self.tag[start+8:stop]
             self.nloweight = NLO_PARTIALWEIGHT(text, self)
-        
+        return self.nloweight
+        
+    def parse_lo_weight(self):
+        """ """
+        if hasattr(self, 'loweight'):
+            return self.loweight
+        
+        start, stop = self.tag.find('<mgrwt>'), self.tag.find('</mgrwt>')
+        
+        if start != -1 != stop :
+            text = self.tag[start+8:stop]
+#<rscale>  3 0.29765919e+03</rscale>
+#<asrwt>0</asrwt>
+#<pdfrwt beam="1">  1       21 0.15134321e+00 0.29765919e+03</pdfrwt>
+#<pdfrwt beam="2">  1       21 0.38683649e-01 0.29765919e+03</pdfrwt>
+#<totfact> 0.17315115e+03</totfact>
+            self.loweight={}
+            for line in text.split('\n'):
+                line = line.replace('<', ' <').replace("'",'"')
+                if 'rscale' in line:
+                    _, nqcd, scale, _ = line.split()
+                    self.loweight['n_qcd'] = int(nqcd)
+                    self.loweight['ren_scale'] = float(scale)
+                elif '<pdfrwt beam="1"' in line:
+                    args = line.split()
+                    self.loweight['n_pdfrw1'] = int(args[2])
+                    npdf = self.loweight['n_pdfrw1']
+                    self.loweight['pdf_pdg_code1'] = [int(i) for i in args[3:3+npdf]]
+                    self.loweight['pdf_x1'] = [float(i) for i in args[3+npdf:3+2*npdf]]
+                    self.loweight['pdf_q1'] = [float(i) for i in args[3+2*npdf:3+3*npdf]]
+                elif '<pdfrwt beam="2"' in line:
+                    args = line.split()
+                    self.loweight['n_pdfrw2'] = int(args[2])
+                    npdf = self.loweight['n_pdfrw2']
+                    self.loweight['pdf_pdg_code2'] = [int(i) for i in args[3:3+npdf]]
+                    self.loweight['pdf_x2'] = [float(i) for i in args[3+npdf:3+2*npdf]]
+                    self.loweight['pdf_q2'] = [float(i) for i in args[3+2*npdf:3+3*npdf]]
+                elif '<asrwt>' in line:
+                    args = line.replace('>','> ').split()
+                    nalps = int(args[1])
+                    self.loweight['asrwt'] = [float(a) for a in args[2:2+nalps]] 
+                    
+                elif 'totfact' in line:
+                    args = line.split() 
+                    self.loweight['tot_fact'] = float(args[1])
+        else:
+            return None
+        return self.loweight
             
     def parse_matching_scale(self):
         """Parse the line containing the starting scale for the shower"""
@@ -1608,10 +1740,40 @@
         for particle in self:
             if particle.status != 1:
                 continue 
-            scale += particle.mass**2 + particle.momentum.pt**2
+            p = FourMomentum(particle)
+            scale += math.sqrt(p.mass_sqr + p.pt**2)
     
         return prefactor * scale
     
+
+    def get_et_scale(self, prefactor=1):
+        
+        scale = 0 
+        for particle in self:
+            if particle.status != 1:
+                continue 
+            p = FourMomentum(particle)
+            pt = p.pt
+            if (pt>0):
+                scale += p.E*pt/math.sqrt(pt**2+p.pz**2)
+    
+        return prefactor * scale    
+    
+    def get_sqrts_scale(self, prefactor=1):
+        
+        scale = 0 
+        init = []
+        for particle in self:
+            if particle.status == -1:
+                init.append(FourMomentum(particle))
+        if len(init) == 1:
+            return init[0].mass
+        elif len(init)==2:
+            return math.sqrt((init[0]+init[1])**2)
+                   
+    
+    
+    
     def get_momenta_str(self, get_order, allow_reversed=True):
         """return the momenta str in the order asked for"""
         
@@ -1785,7 +1947,7 @@
         if power == 1:
             return FourMomentum(self)
         elif power == 2:
-            return self.mass_sqr()
+            return self.mass_sqr
     
     def __repr__(self):
         return 'FourMomentum(%s,%s,%s,%s)' % (self.E, self.px, self.py,self.pz)
@@ -1850,6 +2012,20 @@
         if isinstance(input, str):
             self.parse(input)
         
+    def __str__(self):
+        
+        out = """        pwgt: %(pwgt)s
+        born, real : %(born)s %(real)s
+        pdgs : %(pdgs)s
+        bjks : %(bjks)s
+        scales**2, gs: %(scales2)s %(gs)s
+        born/real related : %(born_related)s %(real_related)s
+        type / nfks : %(type)s  %(nfks)s
+        to merge : %(to_merge_pdg)s in %(merge_new_pdg)s
+        ref_wgt :  %(ref_wgt)s""" % self.__dict__
+        return out
+        
+        
     def parse(self, text):
         """parse the line and create the related object"""
         #0.546601845792D+00 0.000000000000D+00 0.000000000000D+00 0.119210435309D+02 0.000000000000D+00  5 -1 2 -11 12 21 0 0.24546101D-01 0.15706890D-02 0.12586055D+04 0.12586055D+04 0.12586055D+04  1  2  2  2  5  2  2 0.539995789976D+04
@@ -2086,6 +2262,42 @@
         @property
         def aqcd(self):
             return self.event.aqcd
+        
+        def get_ht_scale(self, prefactor=1):
+        
+            scale = 0 
+            for particle in self:
+                p = particle
+                scale += math.sqrt(max(0, p.mass_sqr + p.pt**2))
+        
+            return prefactor * scale
+        
+        def get_et_scale(self, prefactor=1):
+            
+            scale = 0 
+            for particle in self:
+                p = particle
+                pt = p.pt
+                if (pt>0):
+                    scale += p.E*pt/math.sqrt(pt**2+p.pz**2)
+        
+            return prefactor * scale    
+        
+        
+        def get_sqrts_scale(self, event,prefactor=1):
+            
+            scale = 0 
+            nb_init = 0
+            for particle in event:
+                if particle.status == -1:
+                    nb_init+=1
+            if nb_init == 1:
+                return self[0].mass
+            elif nb_init==2:
+                return math.sqrt((self[0]+self[1])**2)
+                   
+    
+        
             
     def __init__(self, input, event):
         
 
=== modified file 'madgraph/various/misc.py'
--- madgraph/various/misc.py	2016-08-07 07:16:07 +0000
+++ madgraph/various/misc.py	2016-08-24 22:15:21 +0000
@@ -1334,6 +1334,7 @@
     if not __debug__:
         return
     
+    use_print = False
     import inspect
     if opt.has_key('cond') and not opt['cond']:
         return
@@ -1346,6 +1347,8 @@
         level = opt['level']
     else:
         level = logging.getLogger('madgraph').level
+        if level == 0:
+            use_print = True
         #print  "madgraph level",level
         #if level == 20:
         #    level = 10 #avoid info level
@@ -1378,12 +1381,18 @@
     else:
         intro = ''
     
-
-    log.log(level, ' '.join([intro]+[str(a) for a in args]) + \
+    if not use_print:
+        log.log(level, ' '.join([intro]+[str(a) for a in args]) + \
                    ' \033[1;30m[%s at line %s]\033[0m' % (os.path.basename(filename), lineno))
+<<<<<<< TREE
 
     if wait:
         raw_input('press_enter to continue')
+=======
+    else:
+        print ' '.join([intro]+[str(a) for a in args]) + \
+                   ' \033[1;30m[%s at line %s]\033[0m' % (os.path.basename(filename), lineno)
+>>>>>>> MERGE-SOURCE
     return 
 
 ################################################################################
@@ -1606,6 +1615,7 @@
 except:
     def apple_notify(subtitle, info_text, userInfo={}):
         return
+<<<<<<< TREE
 ## End apple notify
 
 
@@ -1672,3 +1682,76 @@
     
     
 
+=======
+
+
+python_lhapdf=None
+def import_python_lhapdf(lhapdfconfig):
+    """load the python module of lhapdf return None if it can not be loaded"""
+
+    #save the result to have it faster and avoid the segfault at the second try if lhapdf is not compatible
+    global python_lhapdf
+    if python_lhapdf:
+        if python_lhapdf == -1:
+            return None
+        else:
+            return python_lhapdf
+        
+    use_lhapdf=False
+    try:
+        lhapdf_libdir=subprocess.Popen([lhapdfconfig,'--libdir'],\
+                                           stdout=subprocess.PIPE).stdout.read().strip()
+    except:
+        use_lhapdf=False
+        return False
+    else:
+        try:
+            candidates=[dirname for dirname in os.listdir(lhapdf_libdir) \
+                            if os.path.isdir(os.path.join(lhapdf_libdir,dirname))]
+        except OSError:
+            candidates=[]
+        for candidate in candidates:
+            if os.path.isfile(os.path.join(lhapdf_libdir,candidate,'site-packages','lhapdf.so')):
+                sys.path.insert(0,os.path.join(lhapdf_libdir,candidate,'site-packages'))
+                try:
+                    import lhapdf
+                    use_lhapdf=True
+                    break
+                except ImportError:
+                    sys.path.pop(0)
+                    continue
+    if not use_lhapdf:
+        try:
+            candidates=[dirname for dirname in os.listdir(lhapdf_libdir+'64') \
+                            if os.path.isdir(os.path.join(lhapdf_libdir+'64',dirname))]
+        except OSError:
+            candidates=[]
+        for candidate in candidates:
+            if os.path.isfile(os.path.join(lhapdf_libdir+'64',candidate,'site-packages','lhapdf.so')):
+                sys.path.insert(0,os.path.join(lhapdf_libdir+'64',candidate,'site-packages'))
+                try:
+                    import lhapdf
+                    use_lhapdf=True
+                    break
+                except ImportError:
+                    sys.path.pop(0)
+                    continue
+        if not use_lhapdf:
+            try:
+                import lhapdf
+                use_lhapdf=True
+            except ImportError:
+                print 'fail'
+                logger.warning("Failed to access python version of LHAPDF: "\
+                                   "If the python interface to LHAPDF is available on your system, try "\
+                                   "adding its location to the PYTHONPATH environment variable and the"\
+                                   "LHAPDF library location to LD_LIBRARY_PATH (linux) or DYLD_LIBRARY_PATH (mac os x).")
+        
+    if use_lhapdf:
+        python_lhapdf = lhapdf
+        python_lhapdf.setVerbosity(0)
+    else:
+        python_lhapdf = None
+
+    return python_lhapdf
+>>>>>>> MERGE-SOURCE
 
=== added file 'madgraph/various/systematics.py'
--- madgraph/various/systematics.py	1970-01-01 00:00:00 +0000
+++ madgraph/various/systematics.py	2016-08-24 22:15:21 +0000
@@ -0,0 +1,724 @@
+################################################################################
+#
+# Copyright (c) 2016 The MadGraph5_aMC@NLO Development team and Contributors
+#
+# This file is a part of the MadGraph5_aMC@NLO project, an application which 
+# automatically generates Feynman diagrams and matrix elements for arbitrary
+# high-energy processes in the Standard Model and beyond.
+#
+# It is subject to the MadGraph5_aMC@NLO license which should accompany this 
+# distribution.
+#
+# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
+#
+################################################################################
+from __future__ import division
+if __name__ == "__main__":
+    import sys
+    import os
+    root = os.path.dirname(__file__)
+    if os.path.basename(root) == 'internal':
+        sys.path.append(os.path.dirname(root))
+    else:
+        sys.path.append(os.path.dirname(os.path.dirname(root)))
+        
+import lhe_parser
+import banner
+import banner as banner_mod
+import itertools
+import misc
+import math
+import os 
+import re
+import sys
+import time
+import StringIO
+
+pjoin = os.path.join
+
+class SystematicsError(Exception):
+    pass
+
+class Systematics(object):
+    
+    def __init__(self, input_file, output_file,
+                 start_event=0, stop_event=sys.maxint, write_banner=False,
+                 mur=[0.5,1,2],
+                 muf=[0.5,1,2],
+                 alps=[1],
+                 pdf='errorset', #[(id, subset)]
+                 dyn=[-1,1,2,3,4],
+                 together=[('mur', 'muf', 'dyn')],
+                 lhapdf_config=misc.which('lhapdf-config'),
+                 log=lambda x: sys.stdout.write(str(x)+'\n')
+                 ):
+
+        # INPUT/OUTPUT FILE
+        if isinstance(input_file, str):
+            self.input = lhe_parser.EventFile(input_file)
+        else:
+            self.input = input_file
+        self.output_path = output_file
+        if output_file != None:
+            if isinstance(output_file, str):
+                if output_file == input_file:
+                    directory,name = os.path.split(output_file)
+                    new_name = pjoin(directory, '.tmp_'+name)
+                    self.output =  lhe_parser.EventFile(new_name, 'w')
+                else:
+                    self.output =  lhe_parser.EventFile(output_file, 'w')
+            else:
+                self.output = output_file
+        self.log = log
+        
+        #get some information from the run_card.
+        self.banner = banner_mod.Banner(self.input.banner)  
+        self.force_write_banner = bool(write_banner)
+        self.orig_dyn = self.banner.get('run_card', 'dynamical_scale_choice')
+        self.orig_pdf = self.banner.run_card.get_lhapdf_id()
+    
+        #check for beam
+        beam1, beam2 = self.banner.get_pdg_beam()
+        if abs(beam1) != 2212 and abs(beam2) != 2212:
+            raise SystematicsError, 'can only reweight proton beam'
+        elif abs(beam1) != 2212:
+            self.b1 = 0
+            self.b2 = beam2//2212
+        elif abs(beam2) != 2212:
+            self.b1 = beam1//2212
+            self.b2 = 0
+        else:             
+            self.b1 = beam1//2212
+            self.b2 = beam2//2212
+    
+        if isinstance(self.banner.run_card, banner_mod.RunCardLO):
+            self.is_lo = True
+            if not self.banner.run_card['use_syst']:
+                raise SystematicsError, 'The events was not generated with use_syst=True. Can not evaluate systematics error on this event.'
+        else:
+            self.is_lo = False
+            if not self.banner.run_card['store_rwgt_info']:
+                raise SystematicsError, 'The events was not generated with store_rwgt_info=True. Can not evaluate systematics error on this event.'
+
+        # MUR/MUF/ALPS PARSING
+        if isinstance(mur, str):
+            mur = mur.split(',')
+        self.mur=[float(i) for i in mur]   
+        if isinstance(muf, str):
+            muf = muf.split(',')
+        self.muf=[float(i) for i in muf]
+        
+        if isinstance(alps, str):
+            alps = alps.split(',')
+        self.alps=[float(i) for i in alps]
+
+        # DYNAMICAL SCALE PARSING + together
+        if isinstance(dyn, str):
+            dyn = dyn.split(',')
+        self.dyn=[int(i) for i in dyn]
+ 
+        if isinstance(together, str):
+            self.together = together.split(',')
+        else:
+            self.together = together
+            
+        # START/STOP EVENT                                   
+        self.start_event=int(start_event)
+        self.stop_event=int(stop_event)
+        if start_event != 0:
+            self.log( "#starting from event #%s" % start_event)
+        if stop_event != sys.maxint:
+            self.log( "#stopping at event #%s" % stop_event)
+        
+        # LHAPDF set 
+        if isinstance(lhapdf_config, list):
+            lhapdf_config = lhapdf_config[0]
+        lhapdf = misc.import_python_lhapdf(lhapdf_config)
+        if not lhapdf:
+            return
+        lhapdf.setVerbosity(0)
+        self.pdfsets = {}  
+        if isinstance(pdf, str):
+            pdf = pdf.split(',')
+            
+        if isinstance(pdf,list) and isinstance(pdf[0],(str,int)):
+            self.pdf = []
+            for data in pdf:
+                if data == 'errorset':
+                    data = '%s' % self.orig_pdf
+                if data == 'central':
+                    data = '%s@0' % self.orig_pdf
+                if '@' in data:
+                    #one particular dataset
+                    name, arg = data.rsplit('@',1)
+                    if int(arg) == 0:
+                        if name.isdigit():
+                            self.pdf.append(lhapdf.mkPDF(int(name)))
+                        else:
+                            self.pdf.append(lhapdf.mkPDF(name))
+                    elif name.isdigit():
+                        try:
+                            self.pdf.append(lhapdf.mkPDF(int(name)+int(arg)))
+                        except:
+                            raise Exception, 'invididual error set need to called with name not with lhapdfID'
+                    else:
+                        self.pdf.append(lhapdf.mkPDF(name, int(arg)))
+                else:
+                    if data.isdigit():
+                        pdfset = lhapdf.mkPDF(int(data)).set()
+                    else:
+                        pdfset = lhapdf.mkPDF(data).set()
+                    self.pdfsets[pdfset.lhapdfID] = pdfset 
+                    self.pdf += pdfset.mkPDFs()
+        else:
+            self.pdf = pdf
+            
+        for p in self.pdf:
+            if p.lhapdfID == self.orig_pdf:
+                self.orig_pdf = p
+                break
+        else:
+            self.orig_pdf = lhapdf.mkPDF(self.orig_pdf)
+        self.log( "# events generated with PDF: %s (%s)" %(self.orig_pdf.set().name,self.orig_pdf.lhapdfID ))
+        # create all the function that need to be called
+        self.get_all_fct() # define self.fcts and self.args
+
+    def run(self, stdout=sys.stdout):
+        """ """
+        start_time = time.time()
+        if self.start_event == 0 or self.force_write_banner:
+            lowest_id = self.write_banner(self.output)
+        else:
+            lowest_id = self.get_id()        
+
+        ids = [lowest_id+i for i in range(len(self.args)-1)]
+        all_cross = [0 for i in range(len(self.args))]
+        
+        for nb_event,event in enumerate(self.input):
+            if nb_event < self.start_event:
+                continue
+            elif nb_event >= self.stop_event:
+                if self.force_write_banner:
+                    self.output.write('</LesHouchesEvents>\n')
+                break
+            if self.is_lo:
+                if (nb_event-self.start_event)>=0 and (nb_event-self.start_event) % 2500 ==0:
+                    self.log( '# currently at event %s [ellapsed time: %.2g s]' % (nb_event, time.time()-start_time))
+            else:
+                if (nb_event-self.start_event)>=0 and (nb_event-self.start_event) % 1000 ==0:
+                    self.log( '# currently at event %i [ellapsed time: %.2g s]' % (nb_event, time.time()-start_time))
+                    
+            self.new_event() #re-init the caching of alphas/pdf
+            if self.is_lo:
+                wgts = [self.get_lo_wgt(event, *arg) for arg in self.args]
+            else:
+                wgts = [self.get_nlo_wgt(event, *arg) for arg in self.args]
+            
+            if wgts[0] == 0:
+                print wgts
+                print event
+                raise Exception
+            
+            wgt = [event.wgt*wgts[i]/wgts[0] for i in range(1,len(wgts))]
+            all_cross = [(all_cross[j] + event.wgt*wgts[j]/wgts[0]) for j in range(len(wgts))]
+            
+            rwgt_data = event.parse_reweight()
+            rwgt_data.update(zip(ids, wgt))
+            event.reweight_order += ids
+            # order the 
+            self.output.write(str(event))
+        else:
+            self.output.write('</LesHouchesEvents>\n')
+        self.output.close()
+        self.print_cross_sections(all_cross, min(nb_event,self.stop_event)-self.start_event+1, stdout)
+        
+        if self.output.name != self.output_path:
+            import shutil
+            shutil.move(self.output.name, self.output_path)
+        
+        return all_cross
+        
+    def print_cross_sections(self, all_cross, nb_event, stdout):
+        """print the cross-section."""
+        
+        norm = self.banner.get('run_card', 'event_norm', default='sum')
+        #print "normalisation is ", norm
+        #print "nb_event is ", nb_event
+    
+        max_scale, min_scale = 0,sys.maxint
+        max_alps, min_alps = 0, sys.maxint
+        max_dyn, min_dyn = 0,sys.maxint
+        pdfs = {}
+        dyns = {} # dyn : {'max': , 'min':}
+
+        if norm == 'sum':
+            norm = 1
+        elif norm == 'average':
+            norm = 1./nb_event
+        elif norm == 'unity':
+            norm = 1
+            
+        all_cross = [c*norm for c in all_cross]
+        stdout.write("# mur\t\tmuf\t\talpsfact\tdynamical_scale\tpdf\t\tcross-section\n")
+        for i,arg in enumerate(self.args):
+            
+            to_print = list(arg)
+            to_print[4] = to_print[4].lhapdfID
+            to_print.append(all_cross[i])  
+            to_report = []  
+            stdout.write('%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\n' % tuple(to_print)) 
+            
+            mur, muf, alps, dyn, pdf = arg[:5]
+            if pdf == self.orig_pdf and (dyn==self.orig_dyn or dyn==-1)\
+                                              and (mur!=1 or muf!=1 or alps!=1):
+                max_scale = max(max_scale,all_cross[i])
+                min_scale = min(min_scale,all_cross[i])
+            if pdf == self.orig_pdf and mur==1 and muf==1 and \
+                                    (dyn==self.orig_dyn or dyn==-1) and alps!=1:
+                max_alps = max(max_alps,all_cross[i])
+                min_alps = min(min_alps,all_cross[i]) 
+            
+            if pdf == self.orig_pdf and mur==1 and muf==1 and alps==1:
+                max_dyn = max(max_dyn,all_cross[i])
+                min_dyn = min(min_dyn,all_cross[i])
+                                            
+            if pdf == self.orig_pdf and (alps!=1 or mur!=1 or muf!=1) and \
+                                                (dyn!=self.orig_dyn or dyn!=-1):
+                if dyn not in dyns:
+                    dyns[dyn] = {'max':0, 'min':sys.maxint,'central':0}
+                curr = dyns[dyn]
+                curr['max'] = max(curr['max'],all_cross[i])
+                curr['min'] = min(curr['min'],all_cross[i])
+            if pdf == self.orig_pdf and (alps==1 and mur==1 and muf==1) and \
+                                                (dyn!=self.orig_dyn or dyn!=-1):
+                if dyn not in dyns:
+                    dyns[dyn] = {'max':0, 'min':sys.maxint,'central':all_cross[i]}
+                else:
+                    dyns[dyn]['central'] = all_cross[i]          
+                
+            if alps==1 and mur==1 and muf==1 and (dyn==self.orig_dyn or dyn==-1):
+                pdfset = pdf.set()
+                if pdfset.lhapdfID in self.pdfsets:
+                    if pdfset.lhapdfID not in pdfs :
+                        pdfs[pdfset.lhapdfID] = [0] * pdfset.size
+                    pdfs[pdfset.lhapdfID][pdf.memberID] = all_cross[i]
+                else:
+                    to_report.append('# PDF %s : %s\n' % (pdf.lhapdfID, all_cross[i]))
+  
+        stdout.write('\n') 
+                
+        resume = StringIO.StringIO()
+                
+        resume.write( '#***************************************************************************\n')
+        resume.write( "#\n")
+        resume.write( '# original cross-section: %s\n' % all_cross[0])
+        if max_scale:
+            resume.write( '#     scale variation: +%2.3g%% -%2.3g%%\n' % ((max_scale-all_cross[0])/all_cross[0]*100,(all_cross[0]-min_scale)/all_cross[0]*100))
+        if max_alps:
+            resume.write( '#     emission scale variation: +%2.3g%% -%2.3g%%\n' % ((max_alps-all_cross[0])/all_cross[0]*100,(max_alps-min_scale)/all_cross[0]*100))
+        if max_dyn and (max_dyn!= all_cross[0] or min_dyn != all_cross[0]):
+            resume.write( '#     central scheme variation: +%2.3g%% -%2.3g%%\n' % ((max_dyn-all_cross[0])/all_cross[0]*100,(all_cross[0]-min_dyn)/all_cross[0]*100))
+        if self.orig_pdf.lhapdfID in pdfs:
+            lhapdfid = self.orig_pdf.lhapdfID
+            values = pdfs[lhapdfid]
+            pdfset = self.pdfsets[lhapdfid]
+            pdferr =  pdfset.uncertainty(values)
+            resume.write( '# PDF variation: +%2.3g%% -%2.3g%%\n' % (pdferr.errplus*100/all_cross[0], pdferr.errminus*100/all_cross[0]))       
+        # report error/central not directly linked to the central
+        resume.write( "#\n")        
+        for lhapdfid,values in pdfs.items():
+            if lhapdfid == self.orig_pdf.lhapdfID:
+                continue
+            pdfset = self.pdfsets[lhapdfid]
+            pdferr =  pdfset.uncertainty(values)
+            resume.write( '#PDF %s: %g +%2.3g%% -%2.3g%%\n' % (pdfset.name, pdferr.central,pdferr.errplus*100/all_cross[0], pdferr.errminus*100/all_cross[0]))
+
+        dyn_name = {1: '\sum ET', 2:'\sum\sqrt{m^2+pt^2}', 3:'0.5 \sum\sqrt{m^2+pt^2}',4:'\sqrt{\hat s}' }
+        for key, curr in dyns.items():
+            if key ==-1:
+                continue
+            central, maxvalue, minvalue = curr['central'], curr['max'], curr['min']
+            if central == 0:
+                continue
+            if maxvalue == 0:
+                resume.write("# dynamical scheme # %s : %g # %s\n" %(key, central, dyn_name[key]))
+            else:
+                resume.write("# dynamical scheme # %s : %g +%2.3g%% -%2.3g%% # %s\n" %(key, central, (maxvalue-central)/central*100,(central-minvalue)/central*100, dyn_name[key]))
+      
+        resume.write('\n'.join(to_report))
+
+        resume.write( '#***************************************************************************\n')
+    
+        stdout.write(resume.getvalue())
+        self.log(resume.getvalue())
+    
+    
+    def write_banner(self, fsock):
+        """create the new banner with the information of the weight"""
+
+        cid = self.get_id()
+        lowest_id = cid
+        
+        in_scale = False
+        in_pdf = False
+        in_alps = False
+        
+        text = ''
+        
+        for arg in self.args[1:]:
+            mur, muf, alps, dyn, pdf = arg[:5]
+            if pdf == self.orig_pdf and alps ==1 and (mur!=1 or muf!=1 or dyn!=-1):
+                if not in_scale:
+                    text += "<weightgroup name=\"Central scale variation\" combine=\"envelope\">\n"
+                    in_scale=True
+            elif in_scale:
+                if (pdf == self.orig_pdf and alps ==1):
+                    pass
+                else:
+                    text += "</weightgroup> # scale\n"
+                    in_scale = False
+                
+            if pdf == self.orig_pdf and mur == muf == 1 and dyn==-1 and alps!=1:
+                if not in_alps:
+                    text += "<weightgroup name=\"Emission scale variation\" combine=\"envelope\">\n"
+                    in_alps=True
+            elif in_alps:
+                text += "</weightgroup> # ALPS\n"
+                in_alps=False
+            
+            if pdf != self.orig_pdf and mur == muf == 1 and dyn==-1 and alps ==1:
+                if pdf.lhapdfID in self.pdfsets:
+                    if in_pdf:
+                        text += "</weightgroup> # PDFSET -> PDFSET\n"
+                    pdfset = self.pdfsets[pdf.lhapdfID]
+                    text +="<weightgroup name=\"%s\" combine=\"%s\"> # %s: %s\n" %\
+                            (pdfset.name, pdfset.errorType,pdfset.lhapdfID, pdfset.description)
+                    in_pdf=pdf.lhapdfID 
+                elif pdf.memberID == 1 and (pdf.lhapdfID - pdf.memberID) in self.pdfsets:
+                    if in_pdf:
+                        text += "</weightgroup> # PDFSET -> PDFSET\n"
+                    pdfset = self.pdfsets[pdf.lhapdfID - 1]
+                    text +="<weightgroup name=\"%s\" combine=\"%s\"> # %s: %s\n" %\
+                            (pdfset.name, pdfset.errorType,pdfset.lhapdfID, pdfset.description)
+                    in_pdf=pdfset.lhapdfID 
+                elif in_pdf and pdf.lhapdfID - pdf.memberID != in_pdf:
+                    text += "</weightgroup> # PDFSET -> PDF\n"
+                    in_pdf = False 
+            elif in_pdf:
+                text += "</weightgroup> PDF \n"
+                in_pdf=False                
+                 
+                
+                
+            tag, info = '',''
+            if mur!=1.:
+                tag += 'MUR="%s" ' % mur
+                info += 'MUR=%s ' % mur
+            else:
+                tag += 'MUR="%s" ' % mur
+            if muf!=1.:
+                tag += 'MUF="%s" ' % muf
+                info += 'MUF=%s ' % muf
+            else:
+                tag += 'MUF="%s" ' % muf
+                
+            if alps!=1.:
+                tag += 'ALPSFACT="%s" ' % alps
+                info += 'alpsfact=%s ' % alps
+            if dyn!=-1.:
+                tag += 'DYN_SCALE="%s" ' % dyn
+                info += 'dyn_scale_choice=%s ' % {1:'sum pt', 2:'HT',3:'HT/2',4:'sqrts'}[dyn]
+                                           
+            if pdf != self.orig_pdf:
+                tag += 'LHAPDF="%s" ' % pdf.lhapdfID
+                info += 'LHAPDF=%s MemberID=%s' % (pdf.lhapdfID-pdf.memberID, pdf.memberID)
+            else:
+                tag += 'LHAPDF="%s" ' % pdf.lhapdfID
+                
+            text +='<weight id="%s" %s> %s </weight>\n' % (cid, tag, info)
+            cid+=1
+        
+        if in_scale or in_alps or in_pdf:
+             text += "</weightgroup>\n"
+            
+        if 'initrwgt' in self.banner:
+            self.banner['initrwgt'] += text
+        else:
+            self.banner['initrwgt'] = text
+            
+        
+        self.banner.write(self.output, close_tag=False)
+        
+        return lowest_id
+        
+
+    def get_id(self):
+        
+        if 'initrwgt' in self.banner:
+            pattern = re.compile('<weight id=(?:\'|\")([_\w]+)(?:\'|\")', re.S+re.I+re.M)
+            return  max([int(wid) for wid in  pattern.findall(self.banner['initrwgt'])])+1
+        else:
+            return 1
+        
+        
+
+
+    def get_all_fct(self):
+        
+        all_args = []
+        default = [1.,1.,1.,-1,self.orig_pdf]
+        #all_args.append(default)
+        pos = {'mur':0, 'muf':1, 'alps':2, 'dyn':3, 'pdf':4}
+        done = set()
+        for one_block in self.together:
+            for name in one_block:
+                done.add(name)
+            for together in itertools.product(*[getattr(self,name) for name in one_block]):
+                new_args = list(default)
+                for name,value in zip(one_block, together):
+                    new_args[pos[name]] = value
+                all_args.append(new_args)
+        for name in pos:
+            if name in done:
+                continue
+            for value in getattr(self, name):
+                new_args = list(default)
+                new_args[pos[name]] = value
+                all_args.append(new_args)
+        
+        self.args = [default]+ [arg for arg in all_args if arg!= default]
+
+        self.log( "#Will Compute %s weights per event." % (len(self.args)-1))
+        return
+    
+    def new_event(self):
+        self.alphas = {}
+        self.pdfQ2 = {}
+        
+            
+    def get_pdfQ(self, pdf, pdg, x, scale):
+        
+        if pdg in [-21,-22]:
+            pdg = abs(pdg)
+        elif pdg == 0:
+            return 1
+        
+        f = pdf.xfxQ(pdg, x, scale)/x
+#        if f == 0 and pdf.memberID ==0:
+#            pdfset = pdf.set()
+#            allnumber= [p.xfxQ(pdg, x, scale) for p in pdfset.mkPDFs()]
+#            f = pdfset.uncertainty(allnumber).central /x
+        return f
+
+    def get_pdfQ2(self, pdf, pdg, x, scale):
+
+        if pdg in [-21,-22]:
+            pdg = abs(pdg)
+        elif pdg == 0:
+            return 1
+                
+        if (pdf, pdg,x,scale) in self.pdfQ2:
+            return self.pdfQ2[(pdf, pdg,x,scale)]
+        f = pdf.xfxQ2(pdg, x, scale)/x
+        self.pdfQ2[(pdf, pdg,x,scale)] = f
+        return f        
+        
+        #one method to handle the nnpd2.3 problem -> now just move to central
+        if f == 0 and pdf.memberID ==0:
+            # to avoid problem with nnpdf2.3 in lhapdf6.1.6
+            #print 'central pdf returns 0', pdg, x, scale
+            #print self.pdfsets
+            pdfset = pdf.set()
+            allnumber= [0] + [self.get_pdfQ2(p, pdg, x, scale) for p in pdfset.mkPDFs()[1:]]
+            f = pdfset.uncertainty(allnumber).central
+        self.pdfQ2[(pdf, pdg,x,scale)] = f
+        return f
+                
+    def get_lo_wgt(self,event, Dmur, Dmuf, Dalps, dyn, pdf):
+        """ 
+        pdf is a lhapdf object!"""
+        
+        loinfo = event.parse_lo_weight()
+
+        if dyn == -1:
+            mur = loinfo['ren_scale']
+            muf1 = loinfo['pdf_q1'][-1]
+            muf2 = loinfo['pdf_q2'][-1]
+        else:
+            if dyn == 1: 
+                mur = event.get_et_scale(1.)
+            elif dyn == 2:
+                mur = event.get_ht_scale(1.)
+            elif dyn == 3:
+                mur = event.get_ht_scale(0.5)
+            elif dyn == 4:
+                mur = event.get_sqrts_scale(1.)
+            muf1 = mur
+            muf2 = mur
+            loinfo['pdf_q1'][-1] = mur
+            loinfo['pdf_q2'][-1] = mur
+            
+        
+        # MUR part
+        wgt = pdf.alphasQ(Dmur*mur)**loinfo['n_qcd']
+        # MUF/PDF part
+        wgt *= self.get_pdfQ(pdf, self.b1*loinfo['pdf_pdg_code1'][-1], loinfo['pdf_x1'][-1], Dmuf*muf1) 
+        wgt *= self.get_pdfQ(pdf, self.b2*loinfo['pdf_pdg_code2'][-1], loinfo['pdf_x2'][-1], Dmuf*muf2) 
+        
+        for scale in loinfo['asrwt']:
+            wgt *= pdf.alphasQ(Dalps*scale)
+        
+        # ALS part
+        for i in range(loinfo['n_pdfrw1']-1):
+            scale = min(Dalps*loinfo['pdf_q1'][i], Dmuf*muf1)
+            wgt *= self.get_pdfQ(pdf, self.b1*loinfo['pdf_pdg_code1'][i], loinfo['pdf_x1'][i], scale)
+            wgt /= self.get_pdfQ(pdf, self.b1*loinfo['pdf_pdg_code1'][i], loinfo['pdf_x1'][i+1], scale)
+
+        for i in range(loinfo['n_pdfrw2']-1):
+            scale = min(Dalps*loinfo['pdf_q2'][i], Dmuf*muf2)
+            wgt *= self.get_pdfQ(pdf, self.b2*loinfo['pdf_pdg_code2'][i], loinfo['pdf_x2'][i], scale)
+            wgt /= self.get_pdfQ(pdf, self.b2*loinfo['pdf_pdg_code2'][i], loinfo['pdf_x2'][i+1], scale)            
+        
+        return wgt
+
+    def get_nlo_wgt(self,event, Dmur, Dmuf, Dalps, dyn, pdf):
+        """return the new weight for NLO event --with weight information-- """
+        
+        wgt = 0 
+        nloinfo = event.parse_nlo_weight()
+        for cevent in nloinfo.cevents:
+            if dyn == 1: 
+                mur2 = cevent.get_et_scale(1.)**2
+            elif dyn == 2:
+                mur2 = cevent.get_ht_scale(1.)**2
+            elif dyn == 3:
+                mur2 = cevent.get_ht_scale(0.5)**2
+            elif dyn == 4:
+                mur2 = cevent.get_sqrts_scale(event,1)**2
+            else:
+                mur2 = 0
+            muf2 = mur2
+            
+            for onewgt in cevent.wgts:
+                if dyn == -1:
+                    mur2 = onewgt.scales2[1]
+                    muf2 = onewgt.scales2[2]
+                Q2 = onewgt.scales2[0] # Ellis-Sexton scale
+                
+                wgtpdf = self.get_pdfQ2(pdf, self.b1*onewgt.pdgs[0], onewgt.bjks[0],
+                                      Dmuf**2 * muf2)
+                wgtpdf *= self.get_pdfQ2(pdf, self.b2*onewgt.pdgs[1], onewgt.bjks[1],
+                                      Dmuf**2 * muf2)
+                
+                tmp = onewgt.pwgt[0]
+                tmp += onewgt.pwgt[1] * math.log(Dmur**2 * mur2/ Q2)
+                tmp += onewgt.pwgt[2] * math.log(Dmuf**2 * muf2/ Q2)
+                tmp *= math.sqrt(4*math.pi*pdf.alphasQ2(Dmur**2*mur2))**onewgt.qcdpower
+                
+                if wgtpdf == 0: #happens for nn23pdf due to wrong set in lhapdf
+                    key = (self.b1*onewgt.pdgs[0], self.b2*onewgt.pdgs[1], onewgt.bjks[0],onewgt.bjks[1], muf2)
+                    if dyn== -1 and Dmuf==1 and Dmur==1 and pdf==self.orig_pdf:
+                        wgtpdf = onewgt.ref_wgt / tmp
+                        self.pdfQ2[key] = wgtpdf
+                    elif key in self.pdfQ2:
+                        wgtpdf = self.pdfQ2[key]
+                    else:
+                        # real zero!
+                        wgtpdf = 0
+
+                tmp *= wgtpdf                
+                wgt += tmp
+                
+                if __debug__ and dyn== -1 and Dmur==1 and Dmuf==1 and pdf==self.orig_pdf:
+                    if not misc.equal(tmp, onewgt.ref_wgt, sig_fig=4):
+                        misc.sprint(tmp, onewgt.ref_wgt, (tmp-onewgt.ref_wgt)/tmp)
+                        misc.sprint(onewgt)
+                        misc.sprint(cevent)
+                        raise Exception, 'not enough agreement between stored value and computed one'
+                
+                
+        return wgt
+                            
+
+def call_systematics(args, result=sys.stdout, running=True,
+                     log=lambda x:sys.stdout.write(str(x)+'\n')):
+    """calling systematics from a list of arguments"""            
+
+    input, output = args[0:2]
+    opts = {}
+    for arg in args[2:]:
+        if '=' in arg:
+            key,values= arg.split('=')
+            key = key.replace('-','')
+            values = values.split(',')
+            if key == 'together':
+                if key in opts:
+                    opts[key].append(tuple(values))
+                else:
+                    opts[key]=[tuple(values)]
+            elif key == 'result':
+                result = open(values[0],'w')
+            elif key in ['start_event', 'stop_event']:
+                opts[key] = int(values[0])
+            elif key == 'write_banner':
+                opts[key] = banner_mod.ConfigFile.format_variable(values[0], bool, 'write_banner')
+            else:
+                if key in opts:
+                    opts[key] += values
+                else:
+                    opts[key] = values
+        else:
+            raise SystematicsError, "unknow argument", arg
+
+    #load run_card and extract parameter if needed.
+    if 'from_card' in opts:
+        if opts['from_card'] != ['internal']:
+            card = banner.RunCard(opts['from_card'][0])
+        else:
+            lhe = lhe_parser.EventFile(input)
+            card = banner.RunCard(banner.Banner(lhe.banner)['mgruncard'])
+            lhe.close()
+            
+        if isinstance(card, banner.RunCardLO):
+            # LO case
+            opts['mur'] = [float(x) for x in card['sys_scalefact'].split()]
+            opts['muf'] = opts['mur']
+            if card['sys_alpsfact'] != 'None':
+                opts['alps'] = [float(x) for x in card['sys_alpsfact'].split()]
+            else:
+                opts['alps'] = [1.0]
+            opts['together'] = [('mur','muf','alps','dyn')]
+            pdfs =  card['sys_pdf'].split('&&')
+            opts['dyn'] = [-1,1,2,3,4]
+            opts['pdf'] = []
+            for pdf in pdfs:
+                split = pdf.split()
+                if len(split)==1:
+                    opts['pdf'].append('%s' %pdf)
+                else:
+                    pdf,nb = split
+                    for i in range(int(nb)):
+                        opts['pdf'].append('%s@%s' % (pdf, i))
+            if not opts['pdf']:
+                opts['pdf'] = 'central'
+        else:
+            #NLO case
+            raise Exception
+        del opts['from_card']
+    
+
+    obj = Systematics(input, output, log=log,**opts)
+    if running:
+        obj.run(result)  
+    return obj
+
+if __name__ == "__main__":
+    call_systematics(sys.argv[1:])
+   
+   
+        
+        
+        
+    
+    
+        
 
=== modified file 'tests/acceptance_tests/test_cmd_madevent.py'
--- tests/acceptance_tests/test_cmd_madevent.py	2016-08-17 16:25:19 +0000
+++ tests/acceptance_tests/test_cmd_madevent.py	2016-08-24 22:15:21 +0000
@@ -117,9 +117,10 @@
         interface.exec_cmd('set madanalysis_path %s --no_save' % pjoin(MG5DIR, 'MadAnalysis'))
         interface.onecmd('output madevent %s -f' % self.run_dir)            
         
-        if not os.path.exists(pjoin(interface.options['syscalc_path'],'sys_calc')):
-            print "install SysCalc"
-            interface.onecmd('install SysCalc')
+        if os.path.exists(pjoin(interface.options['syscalc_path'],'sys_calc')):
+            shutil.rmtree(interface.options['syscalc_path'])
+            #print "install SysCalc"
+            #interface.onecmd('install SysCalc')
         
         
         self.cmd_line = MECmd.MadEventCmdShell(me_dir=self.run_dir)
@@ -475,10 +476,10 @@
             # check that the html has the information
             self.assertTrue('syst' in data[0].parton)
             # check that the code was runned correctly
-            fsock = open('%s/Events/%s/%s_parton_syscalc.log' % \
-                  (self.run_dir, data[0]['run_name'], data[0]['tag']),'r')
+            fsock = open('%s/Events/%s/parton_systematics.log' % \
+                  (self.run_dir, data[0]['run_name']),'r')
             text = fsock.read()
-            self.assertTrue(text.count('cross-section') >= 3)
+            self.assertTrue(text.count('dynamical scheme') >= 3)
         
                 
     def check_pythia_output(self, run_name='run_01', syst=False):
@@ -488,9 +489,9 @@
         self.assertTrue('lhe' in data[0].pythia)
         self.assertTrue('log' in data[0].pythia)
 
-        if syst:
-            # check that the html has the information
-            self.assertTrue('rwt' in data[0].pythia)
+#        if syst:
+#            # check that the html has the information
+#            self.assertTrue('rwt' in data[0].pythia)
 
     def check_matched_plot(self, run_name='run_01', mintime=None, tag='fermi'):
         """ """
 
=== modified file 'tests/input_files/run_card_matching.dat'
--- tests/input_files/run_card_matching.dat	2013-11-22 13:16:07 +0000
+++ tests/input_files/run_card_matching.dat	2016-08-24 22:15:21 +0000
@@ -234,9 +234,9 @@
 #  will be use by SysCalc (if install)
 #**************************************                                     
 #
-0.5 = sys_scalefact  # Central scale factors
-0.5 = sys_alpsfact  # \alpha_s emission scale factors
+0.5 1 2= sys_scalefact  # Central scale factors
+1= sys_alpsfact  # \alpha_s emission scale factors
 30 = sys_matchscale # variation of merging scale
 # PDF sets and number of members (0 or none for all members).
-CT10nlo.LHgrid 1 = sys_pdf # matching scales
+CT10nlo 1 = sys_pdf # matching scales
 # MSTW2008nlo68cl.LHgrid 1  = sys_pdf