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

Proposed by Olivier Mattelaer
Status: Merged
Merged at revision: 291
Proposed branch: lp:~maddevelopers/mg5amcnlo/python_syscalc
Merge into: lp:~maddevelopers/mg5amcnlo/2.5.0
Diff against target: 2693 lines (+1721/-339) (has conflicts)
19 files modified
Template/LO/Cards/run_card.dat (+1/-1)
Template/LO/SubProcesses/reweight.f (+6/-6)
UpdateNotes.txt (+15/-0)
VERSION (+6/-0)
madgraph/interface/amcatnlo_run_interface.py (+2/-94)
madgraph/interface/common_run_interface.py (+520/-37)
madgraph/interface/madevent_interface.py (+13/-105)
madgraph/interface/madgraph_interface.py (+2/-2)
madgraph/interface/reweight_interface.py (+95/-38)
madgraph/iolibs/export_fks.py (+2/-1)
madgraph/iolibs/export_v4.py (+2/-0)
madgraph/madevent/gen_crossxhtml.py (+4/-4)
madgraph/various/banner.py (+9/-28)
madgraph/various/cluster.py (+2/-1)
madgraph/various/lhe_parser.py (+220/-8)
madgraph/various/misc.py (+85/-2)
madgraph/various/systematics.py (+724/-0)
tests/acceptance_tests/test_cmd_madevent.py (+10/-9)
tests/input_files/run_card_matching.dat (+3/-3)
Text conflict in UpdateNotes.txt
Text conflict in VERSION
Text conflict in madgraph/various/misc.py
To merge this branch: bzr merge lp:~maddevelopers/mg5amcnlo/python_syscalc
Reviewer Review Type Date Requested Status
Rikkert Frederix Needs Fixing
Eleni Vryonidou Approve
Valentin Hirschi Pending
Review via email: mp+302735@code.launchpad.net

Description of the change

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

To post a comment you must log in.
294. By Olivier Mattelaer

making the code runnning correctly on cluster mode

295. By Olivier Mattelaer

forget to push the last minut changes in the pdf syntax

296. By Olivier Mattelaer

merge with latest 2.5.0 version

Revision history for this message
Valentin Hirschi (valentin-hirschi) wrote :
Download full text (105.4 KiB)

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

297. By Olivier Mattelaer

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

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) wrote :
Download full text (109.7 KiB)

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

Revision history for this message
Valentin Hirschi (valentin-hirschi) wrote :
Download full text (116.9 KiB)

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

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) wrote :

> 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

298. By Olivier Mattelaer

apply change related to round2 of Valentin review

Revision history for this message
Eleni Vryonidou (evryonidou) wrote :
Download full text (293.5 KiB)

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

299. By Olivier Mattelaer

fix an import problem

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) wrote :
Download full text (301.7 KiB)

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

Revision history for this message
Eleni Vryonidou (evryonidou) wrote :
Download full text (303.0 KiB)

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

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) wrote :
Download full text (304.9 KiB)

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

300. By Olivier Mattelaer

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

301. By Olivier Mattelaer

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

302. By Olivier Mattelaer

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

303. By Olivier Mattelaer

remove a print comment

Revision history for this message
Eleni Vryonidou (evryonidou) :
review: Approve
Revision history for this message
Rikkert Frederix (frederix) wrote :

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
Revision history for this message
Rikkert Frederix (frederix) :
review: Needs Fixing

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'Template/LO/Cards/run_card.dat'
2--- Template/LO/Cards/run_card.dat 2015-08-11 00:35:53 +0000
3+++ Template/LO/Cards/run_card.dat 2016-08-24 22:15:21 +0000
4@@ -279,5 +279,5 @@
5 %(sys_alpsfact)s = sys_alpsfact # \alpha_s emission scale factors
6 %(sys_matchscale)s = sys_matchscale # variation of merging scale
7 # PDF sets and number of members (0 or none for all members).
8-%(sys_pdf)s = sys_pdf # matching scales
9+%(sys_pdf)s = sys_pdf # separate by && if more than one set.
10 # MSTW2008nlo68cl.LHgrid 1 = sys_pdf
11
12=== modified file 'Template/LO/SubProcesses/reweight.f'
13--- Template/LO/SubProcesses/reweight.f 2016-06-07 09:33:55 +0000
14+++ Template/LO/SubProcesses/reweight.f 2016-08-24 22:15:21 +0000
15@@ -994,21 +994,21 @@
16 pt2ijcl(nexternal-2)=q2fact(1)
17 if(nexternal.gt.3) pt2ijcl(nexternal-3)=q2fact(1)
18 else
19- q2fact(1)=pt2ijcl(nexternal-2)
20- q2fact(2)=q2fact(1)
21+ q2fact(1)=scalefact**2*pt2ijcl(nexternal-2)
22+ q2fact(2)=scalefact**2*q2fact(1)
23 endif
24 elseif(jcentral(1).eq.0)then
25- q2fact(1) = pt2ijcl(jfirst(1))
26+ q2fact(1) = scalefact**2*pt2ijcl(jfirst(1))
27 elseif(jcentral(2).eq.0)then
28- q2fact(2) = pt2ijcl(jfirst(2))
29+ q2fact(2) = scalefact**2*pt2ijcl(jfirst(2))
30 elseif(ickkw.eq.2.or.(pdfwgt.and.ickkw.gt.0))then
31 c Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
32 c f1(x1,pt2E) is given by DSIG, just need to set scale.
33 c Use the minimum scale found for fact scale in ME
34 if(jlast(1).gt.0.and.jfirst(1).le.jlast(1))
35- $ q2fact(1)=min(pt2ijcl(jfirst(1)),q2fact(1))
36+ $ q2fact(1)=scalefact**2*min(pt2ijcl(jfirst(1)),q2fact(1))
37 if(jlast(2).gt.0.and.jfirst(2).le.jlast(2))
38- $ q2fact(2)=min(pt2ijcl(jfirst(2)),q2fact(2))
39+ $ q2fact(2)=scalefact**2*min(pt2ijcl(jfirst(2)),q2fact(2))
40 endif
41
42 c Check that factorization scale is >= 2 GeV
43
44=== modified file 'UpdateNotes.txt'
45--- UpdateNotes.txt 2016-08-17 19:36:40 +0000
46+++ UpdateNotes.txt 2016-08-24 22:15:21 +0000
47@@ -1,5 +1,6 @@
48 Update notes for MadGraph5_aMC@NLO (in reverse time order)
49
50+<<<<<<< TREE
51 2.5.0
52 OM: Modify the structure of the output format such that all the internal format have the same structure
53 OM: Adding the Plugin directory. Three kind of plugin are currently supported
54@@ -15,6 +16,20 @@
55 OM: add the check that the param_card is compatible with the model restriction.
56 OM+VH:
57
58+=======
59+
60+X.X.X (xx/xx/xx)
61+ OM: NLO/LO Re-weighting works in multi-core
62+ OM: Introduces a new function for LO/NLO interface "systematics"
63+ This function allows to compute systematics uncertainty from the event sample
64+ It requires the event sample to have been generated with
65+ - use_syst = T (for LO sample)
66+ - store_reweight_info = T (for NLO sample)
67+ At LO the code is run automatically if use_syst=T (but if SysCalc is installed)
68+ OM: Fix a bug in the reweight_card where relative path was not understood from the local directory
69+ where the program was runned by the user.
70+
71+>>>>>>> MERGE-SOURCE
72 2.4.3 (01/08/16)
73 OM: Reduce the amount of log file/output generated for LO run (output can use up to three times less output).
74 OM: For the LO combination of events (unweighting) pass to the method previously used for loop-induced.
75
76=== modified file 'VERSION'
77--- VERSION 2016-08-03 13:42:53 +0000
78+++ VERSION 2016-08-24 22:15:21 +0000
79@@ -1,5 +1,11 @@
80+<<<<<<< TREE
81 version = 2.5.0.alpha
82 date = 2016-08-01
83+=======
84+version = new_sys
85+date = 2016-08-08
86+
87+>>>>>>> MERGE-SOURCE
88
89
90
91
92=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
93--- madgraph/interface/amcatnlo_run_interface.py 2016-08-17 16:25:19 +0000
94+++ madgraph/interface/amcatnlo_run_interface.py 2016-08-24 22:15:21 +0000
95@@ -3297,98 +3297,6 @@
96 self.update_status('Run complete', level='shower', update_results=True)
97
98
99- ############################################################################
100- def set_run_name(self, name, tag=None, level='parton', reload_card=False):
101- """define the run name, the run_tag, the banner and the results."""
102-
103- # when are we force to change the tag new_run:previous run requiring changes
104- upgrade_tag = {'parton': ['parton','pythia','pgs','delphes','shower'],
105- 'pythia': ['pythia','pgs','delphes'],
106- 'shower': ['shower'],
107- 'pgs': ['pgs'],
108- 'delphes':['delphes'],
109- 'plot':[]}
110-
111-
112-
113- if name == self.run_name:
114- if reload_card:
115- run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
116- self.run_card = banner_mod.RunCardNLO(run_card)
117-
118- #check if we need to change the tag
119- if tag:
120- self.run_card['run_tag'] = tag
121- self.run_tag = tag
122- self.results.add_run(self.run_name, self.run_card)
123- else:
124- for tag in upgrade_tag[level]:
125- if getattr(self.results[self.run_name][-1], tag):
126- tag = self.get_available_tag()
127- self.run_card['run_tag'] = tag
128- self.run_tag = tag
129- self.results.add_run(self.run_name, self.run_card)
130- break
131- return # Nothing to do anymore
132-
133- # save/clean previous run
134- if self.run_name:
135- self.store_result()
136- # store new name
137- self.run_name = name
138-
139- # Read run_card
140- run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
141- self.run_card = banner_mod.RunCardNLO(run_card)
142-
143- new_tag = False
144- # First call for this run -> set the banner
145- self.banner = banner_mod.recover_banner(self.results, level, self.run_name, tag)
146- if tag:
147- self.run_card['run_tag'] = tag
148- new_tag = True
149- elif not self.run_name in self.results and level =='parton':
150- pass # No results yet, so current tag is fine
151- elif not self.run_name in self.results:
152- #This is only for case when you want to trick the interface
153- logger.warning('Trying to run data on unknown run.')
154- self.results.add_run(name, self.run_card)
155- self.results.update('add run %s' % name, 'all', makehtml=True)
156- else:
157- for tag in upgrade_tag[level]:
158-
159- if getattr(self.results[self.run_name][-1], tag):
160- # LEVEL is already define in the last tag -> need to switch tag
161- tag = self.get_available_tag()
162- self.run_card['run_tag'] = tag
163- new_tag = True
164- break
165- if not new_tag:
166- # We can add the results to the current run
167- tag = self.results[self.run_name][-1]['tag']
168- self.run_card['run_tag'] = tag # ensure that run_tag is correct
169-
170-
171- if name in self.results and not new_tag:
172- self.results.def_current(self.run_name)
173- else:
174- self.results.add_run(self.run_name, self.run_card)
175-
176- self.run_tag = self.run_card['run_tag']
177-
178- # Return the tag of the previous run having the required data for this
179- # tag/run to working wel.
180- if level == 'parton':
181- return
182- elif level == 'pythia':
183- return self.results[self.run_name][0]['tag']
184- else:
185- for i in range(-1,-len(self.results[self.run_name])-1,-1):
186- tagRun = self.results[self.run_name][i]
187- if tagRun.pythia:
188- return tagRun['tag']
189-
190-
191 def store_result(self):
192 """ tar the pythia results. This is done when we are quite sure that
193 the pythia output will not be use anymore """
194@@ -4928,9 +4836,9 @@
195 if '--web' in args:
196 i = args.index('--web')
197 args.pop(i)
198- cmd_line = aMCatNLOCmd(force_run=True)
199+ cmd_line = aMCatNLOCmd(me_dir=os.path.dirname(root_path),force_run=True)
200 else:
201- cmd_line = aMCatNLOCmdShell(force_run=True)
202+ cmd_line = aMCatNLOCmdShell(me_dir=os.path.dirname(root_path),force_run=True)
203
204 if not hasattr(cmd_line, 'do_%s' % args[0]):
205 if parser_error:
206
207=== modified file 'madgraph/interface/common_run_interface.py'
208--- madgraph/interface/common_run_interface.py 2016-08-17 16:25:19 +0000
209+++ madgraph/interface/common_run_interface.py 2016-08-24 22:15:21 +0000
210@@ -68,6 +68,7 @@
211 import internal.files as files
212 import internal.save_load_object as save_load_object
213 import internal.gen_crossxhtml as gen_crossxhtml
214+ import internal.lhe_parser as lhe_parser
215 from internal import InvalidCmd, MadGraph5Error
216 MADEVENT=True
217 else:
218@@ -78,6 +79,7 @@
219 import madgraph.various.misc as misc
220 import madgraph.iolibs.files as files
221 import madgraph.various.cluster as cluster
222+ import madgraph.various.lhe_parser as lhe_parser
223 import madgraph.iolibs.save_load_object as save_load_object
224 import madgraph.madevent.gen_crossxhtml as gen_crossxhtml
225 import models.check_param_card as check_param_card
226@@ -450,7 +452,6 @@
227 """return the path to the output events
228 """
229
230-
231 if self.mode == 'madevent':
232 possible_path = [
233 pjoin(self.me_dir,'Events', run_name, 'unweighted_events.lhe.gz'),
234@@ -465,7 +466,10 @@
235 correct_path = path
236 break
237 else:
238- raise self.InvalidCmd('No events file corresponding to %s run. ' % run_name)
239+ if os.path.exists(run_name):
240+ correct_path = run_name
241+ else:
242+ raise self.InvalidCmd('No events file corresponding to %s run. ' % run_name)
243 return correct_path
244
245
246@@ -619,6 +623,101 @@
247
248 return args
249
250+ ############################################################################
251+ def set_run_name(self, name, tag=None, level='parton', reload_card=False,
252+ allow_new_tag=True):
253+ """define the run name, the run_tag, the banner and the results."""
254+
255+ # when are we force to change the tag new_run:previous run requiring changes
256+ upgrade_tag = {'parton': ['parton','pythia','pgs','delphes'],
257+ 'pythia': ['pythia','pgs','delphes'],
258+ 'pgs': ['pgs'],
259+ 'delphes':['delphes'],
260+ 'plot':[],
261+ 'syscalc':[]}
262+
263+
264+
265+ if name == self.run_name:
266+ if reload_card:
267+ run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
268+ self.run_card = banner_mod.RunCard(run_card)
269+
270+ #check if we need to change the tag
271+ if tag:
272+ self.run_card['run_tag'] = tag
273+ self.run_tag = tag
274+ self.results.add_run(self.run_name, self.run_card)
275+ else:
276+ for tag in upgrade_tag[level]:
277+ if getattr(self.results[self.run_name][-1], tag):
278+ tag = self.get_available_tag()
279+ self.run_card['run_tag'] = tag
280+ self.run_tag = tag
281+ self.results.add_run(self.run_name, self.run_card)
282+ break
283+ return # Nothing to do anymore
284+
285+ # save/clean previous run
286+ if self.run_name:
287+ self.store_result()
288+ # store new name
289+ self.run_name = name
290+
291+ new_tag = False
292+ # First call for this run -> set the banner
293+ self.banner = banner_mod.recover_banner(self.results, level, name)
294+ if 'mgruncard' in self.banner:
295+ self.run_card = self.banner.charge_card('run_card')
296+ else:
297+ # Read run_card
298+ run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
299+ self.run_card = banner_mod.RunCard(run_card)
300+
301+ if tag:
302+ self.run_card['run_tag'] = tag
303+ new_tag = True
304+ elif not self.run_name in self.results and level =='parton':
305+ pass # No results yet, so current tag is fine
306+ elif not self.run_name in self.results:
307+ #This is only for case when you want to trick the interface
308+ logger.warning('Trying to run data on unknown run.')
309+ self.results.add_run(name, self.run_card)
310+ self.results.update('add run %s' % name, 'all', makehtml=False)
311+ else:
312+ for tag in upgrade_tag[level]:
313+
314+ if getattr(self.results[self.run_name][-1], tag):
315+ # LEVEL is already define in the last tag -> need to switch tag
316+ tag = self.get_available_tag()
317+ self.run_card['run_tag'] = tag
318+ new_tag = True
319+ break
320+ if not new_tag:
321+ # We can add the results to the current run
322+ tag = self.results[self.run_name][-1]['tag']
323+ self.run_card['run_tag'] = tag # ensure that run_tag is correct
324+
325+ if allow_new_tag and (name in self.results and not new_tag):
326+ self.results.def_current(self.run_name)
327+ else:
328+ self.results.add_run(self.run_name, self.run_card)
329+
330+ self.run_tag = self.run_card['run_tag']
331+
332+ # Return the tag of the previous run having the required data for this
333+ # tag/run to working wel.
334+ if level == 'parton':
335+ return
336+ elif level == 'pythia':
337+ return self.results[self.run_name][0]['tag']
338+ else:
339+ for i in range(-1,-len(self.results[self.run_name])-1,-1):
340+ tagRun = self.results[self.run_name][i]
341+ if tagRun.pythia:
342+ return tagRun['tag']
343+
344+
345 @misc.multiple_try(nb_try=5, sleep=2)
346 def load_results_db(self):
347 """load the current results status"""
348@@ -1087,7 +1186,270 @@
349 """Dummy routine, to be overwritten by daughter classes"""
350
351 pass
352-
353+
354+ ############################################################################
355+ def help_systematics(self):
356+ """help for systematics command"""
357+ logger.info("syntax: systematics RUN_NAME [OUTPUT] [options]",'$MG:color:BLACK')
358+ logger.info("-- Run the systematics run on the run_name run.")
359+ logger.info(" RUN_NAME can be a path to a lhef file.")
360+ logger.info(" OUTPUT can be the path to the output lhe file. we overwritte the input file otherwise")
361+ logger.info("")
362+ logger.info("options: (value written are default)", '$MG:color:BLACK')
363+ logger.info("")
364+ logger.info(" --mur=0.5,1,2 # specify the values for renormalisation scale variation")
365+ logger.info(" --muf=0.5,1,2 # specify the values for factorisation scale variation")
366+ logger.info(" --alps=1 # specify the values for MLM emission scale variation")
367+ logger.info(" --dyn=-1,1,2,3,4 # specify the dynamical schemes to use.")
368+ logger.info(" # -1 is the one used by the sample.")
369+ logger.info(" # > 0 correspond to options of dynamical_scale_choice of the run_card.")
370+ logger.info(" --pdf=errorset # specify the pdfs to use for pdf variation. (see below)")
371+ logger.info(" --together=mur,muf,dyn # lists the parameter that must be varied simultaneously so as to ")
372+ logger.info(" # compute the weights for all combinations of their variations.")
373+ logger.info(" --from_card # use the information from the run_card (LO only).")
374+ logger.info("")
375+ logger.info(" Allowed value for the pdf options:", '$MG:color:BLACK')
376+ logger.info(" central : Do not perform any pdf variation" )
377+ logger.info(" errorset : runs over the errorset")
378+ logger.info(" 244800 : runs over the associated set and his errorset")
379+ logger.info(" 244800@0 : runs over the associated set ONLY")
380+# logger.info(" 244800@X : runs over the Xth set of the associated error set")
381+ logger.info(" CT10 : runs over the associated set and his errorset")
382+ logger.info(" CT10@0 : runs over the associated set ONLY")
383+ logger.info(" CT10@X : runs over the Xth set of the associated error set")
384+ logger.info(" XX,YY,ZZ : runs over the sets for XX,YY,ZZ (those three follows above syntax)")
385+
386+ def complete_systematics(self, text, line, begidx, endidx):
387+ """auto completion for the systematics command"""
388+
389+ args = self.split_arg(line[0:begidx], error=False)
390+ options = ['--mur=', '--muf=', '--pdf=', '--dyn=','--alps=','--together=','--from_card ']
391+
392+ if len(args) == 1 and os.path.sep not in text:
393+ #return valid run_name
394+ data = misc.glob(pjoin('*','*events.lhe*'), pjoin(self.me_dir, 'Events'))
395+ data = [n.rsplit('/',2)[1] for n in data]
396+ return self.list_completion(text, data, line)
397+ elif len(args)==1:
398+ #logger.warning('1args')
399+ return self.path_completion(text,
400+ os.path.join('.',*[a for a in args \
401+ if a.endswith(os.path.sep)]))
402+ elif len(args)==2 and os.path.sep in args[1]:
403+ #logger.warning('2args %s', args[1])
404+ return self.path_completion(text, '.')
405+
406+ elif not line.endswith(tuple(options)):
407+ return self.list_completion(text, options)
408+
409+
410+ ############################################################################
411+ def do_systematics(self, line):
412+ """ syntax is 'systematics [INPUT [OUTPUT]] OPTIONS'
413+ --mur=0.5,1,2
414+ --muf=0.5,1,2
415+ --alps=1
416+ --dyn=-1
417+ --together=mur,muf #can be repeated
418+
419+ #special options
420+ --from_card=
421+ """
422+
423+ lhapdf = misc.import_python_lhapdf(self.options['lhapdf'])
424+ if not lhapdf:
425+ logger.info('can not run systematics since can not link python to lhapdf')
426+ return
427+
428+ self.update_status('Running Systematic computation', level='parton')
429+ args = self.split_arg(line)
430+ #split arguments and option
431+ opts= []
432+ args = [a for a in args if not a.startswith('-') or opts.append(a)]
433+
434+ #check sanity of options
435+ if any(not o.startswith(('--mur=', '--muf=', '--alps=','--dyn=','--together=','--from_card','--pdf='))
436+ for o in opts):
437+ raise self.InvalidCmd, "command systematics called with invalid option syntax. Please retry."
438+
439+ # check that we have define the input
440+ if len(args) == 0:
441+ if self.run_name:
442+ args[0] = self.run_name
443+ else:
444+ raise self.InvalidCmd, 'no default run. Please specify the run_name'
445+
446+ # always pass to a path + get the event size
447+ result_file= sys.stdout
448+ if not os.path.sep in args[0]:
449+ path = [pjoin(self.me_dir, 'Events', args[0], 'unweighted_events.lhe.gz'),
450+ pjoin(self.me_dir, 'Events', args[0], 'unweighted_events.lhe'),
451+ pjoin(self.me_dir, 'Events', args[0], 'events.lhe.gz'),
452+ pjoin(self.me_dir, 'Events', args[0], 'events.lhe')]
453+
454+ for p in path:
455+ if os.path.exists(p):
456+ nb_event = self.results[args[0]].get_current_info()['nb_event']
457+
458+
459+ if self.run_name != args[0]:
460+ tag = self.results[args[0]].tags[0]
461+ self.set_run_name(args[0], tag,'parton', False)
462+ result_file = open(pjoin(self.me_dir,'Events', self.run_name, 'parton_systematics.log'),'w')
463+ args[0] = p
464+ break
465+ else:
466+ raise self.InvalidCmd, 'Invalid run name. Please retry'
467+ elif self.options['nb_core'] != 1:
468+ lhe = lhe_parser.EventFile(args)
469+ nb_event = len(lhe)
470+ lhe.close()
471+
472+ input = args[0]
473+ if len(args)>1:
474+ output = pjoin(os.getcwd(),args[1])
475+ else:
476+ output = input
477+
478+ lhaid = [self.run_card.get_lhapdf_id()]
479+ if 'store_rwgt_info' in self.run_card and not self.run_card['store_rwgt_info']:
480+ raise self.InvalidCmd, "The events was not generated with store_rwgt_info=True. Can not evaluate systematics error on this event file."
481+ elif 'use_syst' in self.run_card and not self.run_card['use_syst']:
482+ raise self.InvalidCmd, "The events was not generated with use_syst=True. Can not evaluate systematics error on this event file."
483+
484+ try:
485+ pdfsets_dir = self.get_lhapdf_pdfsetsdir()
486+ except Exception, error:
487+ logger.debug(str(error))
488+ logger.warning('Systematic computation requires lhapdf to run. Bypass Systematics')
489+ return
490+
491+ if '--from_card' in opts:
492+ opts.remove('--from_card')
493+ opts.append('--from_card=internal')
494+
495+ # Check that all pdfset are correctly installed
496+ if 'sys_pdf' in self.run_card:
497+ sys_pdf = self.run_card['sys_pdf'].split('&&')
498+ lhaid += [l.split()[0] for l in sys_pdf]
499+
500+ else:
501+ #check that all p
502+ pdf = [a[6:] for a in opts if a.startswith('--pdf=')]
503+ lhaid += [t.split('@')[0] for p in pdf for t in p.split(',')
504+ if t not in ['errorset', 'central']]
505+
506+ # Copy all the relevant PDF sets
507+ [self.copy_lhapdf_set([onelha], pdfsets_dir) for onelha in lhaid]
508+
509+
510+ if self.options['run_mode'] ==2:
511+ nb_submit = min(self.options['nb_core'], nb_event//2500)
512+ elif self.options['run_mode'] ==1:
513+ nb_submit = min(self.options['cluster_size'], nb_event//25000)
514+ else:
515+ nb_submit =1
516+
517+ if MADEVENT:
518+ import internal.systematics as systematics
519+ else:
520+ import madgraph.various.systematics as systematics
521+
522+ #one core:
523+ if nb_submit in [0,1]:
524+ systematics.call_systematics([input, output] + opts,
525+ log=lambda x: logger.info(str(x)),
526+ result=result_file
527+ )
528+
529+ elif self.options['run_mode'] in [1,2]:
530+ event_per_job = nb_event // nb_submit
531+ nb_job_with_plus_one = nb_event % nb_submit
532+ start_event, stop_event = 0,0
533+ for i in range(nb_submit):
534+ #computing start/stop event
535+ event_requested = event_per_job
536+ if i < nb_job_with_plus_one:
537+ event_requested += 1
538+ start_event = stop_event
539+ stop_event = start_event + event_requested
540+
541+ prog = sys.executable
542+ input_files = [os.path.basename(input)]
543+ output_files = ['./tmp_%s_%s' % (i, os.path.basename(output)),
544+ './log_sys_%s.txt' % (i)]
545+ argument = []
546+ if not __debug__:
547+ argument.append('-O')
548+ argument += [pjoin(self.me_dir, 'bin', 'internal', 'systematics.py'),
549+ input_files[0], output_files[0]] + opts +\
550+ ['--start_event=%i' % start_event,
551+ '--stop_event=%i' %stop_event,
552+ '--result=./log_sys_%s.txt' %i,
553+ '--lhapdf_config=%s' % self.options['lhapdf']]
554+ required_output = output_files
555+ self.cluster.cluster_submit(prog, argument,
556+ input_files=input_files,
557+ output_files=output_files,
558+ cwd=os.path.dirname(output),
559+ required_output=required_output,
560+ stdout='/dev/null'
561+ )
562+ starttime = time.time()
563+ update_status = lambda idle, run, finish: \
564+ self.update_status((idle, run, finish, 'running systematics'), level=None,
565+ force=False, starttime=starttime)
566+
567+ self.cluster.wait(os.path.dirname(output), update_status, update_first=update_status)
568+
569+ #collect the data
570+ all_cross = []
571+ for i in range(nb_submit):
572+ pos=0
573+ for line in open(pjoin(os.path.dirname(output), 'log_sys_%s.txt'%i)):
574+ if line.startswith('#'):
575+ continue
576+ split = line.split()
577+ if len(split) in [0,1]:
578+ continue
579+ key = tuple(float(x) for x in split[:-1])
580+ cross= float(split[-1])
581+ if 'event_norm' in self.run_card and \
582+ self.run_card['event_norm'] in ['average', 'unity']:
583+ cross *= (event_per_job+1 if i <nb_job_with_plus_one else event_per_job)
584+ if len(all_cross) > pos:
585+ all_cross[pos] += cross
586+ else:
587+ all_cross.append(cross)
588+ pos+=1
589+
590+ if 'event_norm' in self.run_card and \
591+ self.run_card['event_norm'] in ['unity']:
592+ all_cross= [cross/nb_event for cross in all_cross]
593+
594+ sys_obj = systematics.call_systematics([input, None] + opts,
595+ log=lambda x: logger.info(str(x)),
596+ result=result_file,
597+ running=False
598+ )
599+ sys_obj.print_cross_sections(all_cross, nb_event, result_file)
600+
601+ #concatenate the output file
602+ subprocess.call(['cat']+\
603+ ['./tmp_%s_%s' % (i, os.path.basename(output)) for i in range(nb_submit)],
604+ stdout=open(output,'w'),
605+ cwd=os.path.dirname(output))
606+ for i in range(nb_submit):
607+ os.remove('%s/tmp_%s_%s' %(os.path.dirname(output),i,os.path.basename(output)))
608+ # os.remove('%s/log_sys_%s.txt' % (os.path.dirname(output),i))
609+
610+
611+
612+
613+
614+ self.update_status('End of systematics computation', level='parton', makehtml=False)
615+
616+
617 ############################################################################
618 def do_reweight(self, line):
619 """ syntax is "reweight RUN_NAME"
620@@ -1096,9 +1458,50 @@
621 cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Reweight
622 """
623
624+ #### Utility function
625+ def check_multicore(self):
626+ """ determine if the cards are save for multicore use"""
627+ card = pjoin(self.me_dir, 'Cards', 'reweight_card.dat')
628+
629+ multicore = True
630+ if self.options['run_mode'] in [0,1]:
631+ multicore = False
632+
633+ lines = [l.strip() for l in open(card) if not l.strip().startswith('#')]
634+ while lines and not lines[0].startswith('launch'):
635+ line = lines.pop(0)
636+ # if not standard output mode forbid multicore mode
637+ if line.startswith('change') and line[6:].strip().startswith('output'):
638+ return False
639+ if line.startswith('change') and line[6:].strip().startswith('multicore'):
640+ split_line = line.split()
641+ if len(split_line) > 2:
642+ multicore = bool(split_line[2])
643+ # we have reached the first launch in the card ensure that no output change
644+ #are done after that point.
645+ lines = [line[6:].strip() for line in lines if line.startswith('change')]
646+ for line in lines:
647+ if line.startswith(('process','model','output', 'rwgt_dir')):
648+ return False
649+ elif line.startswith('multicore'):
650+ split_line = line.split()
651+ if len(split_line) > 1:
652+ multicore = bool(split_line[1])
653+
654+ return multicore
655+
656+
657+
658 if '-from_cards' in line and not os.path.exists(pjoin(self.me_dir, 'Cards', 'reweight_card.dat')):
659 return
660-
661+ # option for multicore to avoid that all of them create the same directory
662+ if '--multicore=create' in line:
663+ multicore='create'
664+ elif '--multicore=wait' in line:
665+ multicore='wait'
666+ else:
667+ multicore=False
668+
669 # Check that MG5 directory is present .
670 if MADEVENT and not self.options['mg5_path']:
671 raise self.InvalidCmd, '''The module reweight requires that MG5 is installed on the system.
672@@ -1125,6 +1528,10 @@
673 if self.run_name and self.results.current and self.results.current['cross'] == 0:
674 self.results.delete_run(self.run_name, self.run_tag)
675 self.results.save()
676+ # ensure that the run_card is present
677+ if not hasattr(self, 'run_card'):
678+ self.run_card = banner_mod.RunCard(pjoin(self.me_dir, 'Cards', 'run_card.dat'))
679+
680 # we want to run this in a separate shell to avoid hard f2py crash
681 command = [sys.executable]
682 if os.path.exists(pjoin(self.me_dir, 'bin', 'madevent')):
683@@ -1134,41 +1541,113 @@
684 if not isinstance(self, cmd.CmdShell):
685 command.append('--web')
686 command.append('reweight')
687- if self.run_name:
688- command.append(self.run_name)
689+
690+ ######### START SINGLE CORE MODE ############
691+ if self.options['nb_core']==1 or self.run_card['nevents'] < 101 or not check_multicore(self):
692+ if self.run_name:
693+ command.append(self.run_name)
694+ else:
695+ command += args
696+ if '-from_cards' not in command:
697+ command.append('-from_cards')
698+ p = misc.Popen(command, stdout = subprocess.PIPE, stderr = subprocess.STDOUT, cwd=os.getcwd())
699+ while p.poll() is None:
700+ line = p.stdout.readline()
701+ if any(t in line for t in ['INFO:', 'WARNING:', 'CRITICAL:', 'ERROR:', 'root:','KEEP:']) and \
702+ not '***********' in line:
703+ print line[:-1].replace('INFO', 'REWEIGHT').replace('KEEP:','')
704+ elif __debug__ and line:
705+ logger.debug(line[:-1])
706+ if p.returncode !=0:
707+ logger.error("Reweighting failed")
708+ return
709+ self.results = self.load_results_db()
710+ # forbid this function to create an empty item in results.
711+ try:
712+ if self.results[self.run_name][-2]['cross']==0:
713+ self.results.delete_run(self.run_name,self.results[self.run_name][-2]['tag'])
714+ except:
715+ pass
716+ try:
717+ if self.results.current['cross'] == 0 and self.run_name:
718+ self.results.delete_run(self.run_name, self.run_tag)
719+ except:
720+ pass
721+ # re-define current run
722+ try:
723+ self.results.def_current(self.run_name, self.run_tag)
724+ except Exception:
725+ pass
726+ return
727+ ########## END SINGLE CORE HANDLING #############
728 else:
729- command += args
730- if '-from_cards' not in command:
731- command.append('-from_cards')
732- p = misc.Popen(command, stdout = subprocess.PIPE, stderr = subprocess.STDOUT, cwd=self.me_dir)
733- while p.poll() is None:
734- line = p.stdout.readline()
735- if any(t in line for t in ['INFO:', 'WARNING:', 'CRITICAL:', 'ERROR:', 'root:','KEEP:']) and \
736- not '***********' in line:
737- print line[:-1].replace('INFO', 'REWEIGHT').replace('KEEP:','')
738- elif __debug__ and line:
739- logger.debug(line[:-1])
740- if p.returncode !=0:
741- logger.error("Reweighting failed")
742+ ########## START MULTI-CORE HANDLING #############
743+ if not isinstance(self.cluster, cluster.MultiCore):
744+ mycluster = cluster.MultiCore(nb_core=self.options['nb_core'])
745+ else:
746+ mycluster = self.cluster
747+
748+ new_args=list(args)
749+ self.check_decay_events(new_args)
750+ try:
751+ os.remove(pjoin(self.me_dir,'rw_me','rwgt.pkl'))
752+ except Exception, error:
753+ pass
754+ # prepare multi-core job:
755+ import madgraph.various.lhe_parser as lhe_parser
756+ # args now alway content the path to the valid files
757+ if 'nevt_job' in self.run_card and self.run_card['nevt_job'] !=-1:
758+ nevt_job = self.run_card['nevt_job']
759+ else:
760+ nevt_job = max(5000, self.run_card['nevents']/50)
761+ logger.info("split the event file in bunch of %s events" % nevt_job)
762+ nb_file = lhe_parser.EventFile(new_args[0]).split(nevt_job)
763+ starttime = time.time()
764+ update_status = lambda idle, run, finish: \
765+ self.update_status((idle, run, finish, 'reweight'), level=None,
766+ force=False, starttime=starttime)
767+
768+ all_lhe = []
769+ devnull= open(os.devnull)
770+ for i in range(nb_file):
771+ new_command = list(command)
772+ new_command.append('%s_%s.lhe' % (new_args[0],i))
773+ all_lhe.append('%s_%s.lhe' % (new_args[0],i))
774+ if '-from_cards' not in command:
775+ new_command.append('-from_cards')
776+ if i==0:
777+ if __debug__:
778+ stdout = None
779+ else:
780+ stdout = open(pjoin(self.me_dir,'Events', self.run_name, 'reweight.log'),'w')
781+ new_command.append('--multicore=create')
782+ else:
783+ stdout = devnull
784+ #stdout = open(pjoin(self.me_dir,'Events', self.run_name, 'reweight%s.log' % i),'w')
785+ new_command.append('--multicore=wait')
786+ mycluster.submit(prog=command[0], argument=new_command[1:], stdout=stdout, cwd=os.getcwd())
787+ mycluster.wait(self.me_dir,update_status)
788+ devnull.close()
789+
790+ lhe = lhe_parser.MultiEventFile(all_lhe, parse=False)
791+ nb_event, cross_sections = lhe.write(new_args[0], get_info=True)
792+ if any(os.path.exists('%s_%s_debug.log' % (f, self.run_tag)) for f in all_lhe):
793+ for f in all_lhe:
794+ if os.path.exists('%s_%s_debug.log' % (f, self.run_tag)):
795+ raise Exception, "Some of the run failed: Please read %s_%s_debug.log" % (f, self.run_tag)
796+
797+
798+ if 'event_norm' in self.run_card and self.run_card['event_norm'] == 'average':
799+ for key, value in cross_sections.items():
800+ cross_sections[key] = value / (nb_event+1)
801+ lhe.remove()
802+ for key in cross_sections:
803+ if key == 'orig' or key.isdigit():
804+ continue
805+ logger.info('%s : %s pb' % (key, cross_sections[key]))
806 return
807- self.results = self.load_results_db()
808- # forbid this function to create an empty item in results.
809- try:
810- if self.results[self.run_name][-2]['cross']==0:
811- self.results.delete_run(self.run_name,self.results[self.run_name][-2]['tag'])
812- except:
813- pass
814- try:
815- if self.results.current['cross'] == 0 and self.run_name:
816- self.results.delete_run(self.run_name, self.run_tag)
817- except:
818- pass
819- # re-define current run
820- try:
821- self.results.def_current(self.run_name, self.run_tag)
822- except Exception:
823- pass
824- return
825+ ########## END MULTI-CORE HANDLING #############
826+
827
828 self.to_store.append('event')
829 # forbid this function to create an empty item in results.
830@@ -1189,6 +1668,8 @@
831 path = pjoin(self.me_dir, 'Cards', 'reweight_card.dat')
832 reweight_cmd.raw_input=False
833 reweight_cmd.me_dir = self.me_dir
834+ reweight_cmd.multicore = multicore #allow the directory creation or not
835+ print "We are in mode", multicore
836 reweight_cmd.import_command_file(path)
837 reweight_cmd.do_quit('')
838
839@@ -2319,6 +2800,8 @@
840
841 pdfsetname=set()
842 for lhaid in lhaid_list:
843+ if isinstance(lhaid, str) and lhaid.isdigit():
844+ lhaid = int(lhaid)
845 if isinstance(lhaid, (int,float)):
846 try:
847 if lhaid in self.lhapdf_pdfsets:
848
849=== modified file 'madgraph/interface/madevent_interface.py'
850--- madgraph/interface/madevent_interface.py 2016-08-17 19:36:40 +0000
851+++ madgraph/interface/madevent_interface.py 2016-08-24 22:15:21 +0000
852@@ -2012,11 +2012,17 @@
853 if not bypass_run:
854 self.exec_cmd('refine %s' % nb_event, postcmd=False)
855
856- self.exec_cmd('combine_events', postcmd=False)
857+ self.exec_cmd('combine_events', postcmd=False,printcmd=False)
858 self.print_results_in_shell(self.results.current)
859
860-
861- self.run_syscalc('parton')
862+ if self.run_card['use_syst']:
863+ scdir = self.options['syscalc_path']
864+ if not scdir or not os.path.exists(scdir):
865+ self.exec_cmd('systematics %s --from_card' % self.run_name, postcmd=False,printcmd=False)
866+ else:
867+ self.run_syscalc('parton')
868+
869+
870 self.create_plot('parton')
871 self.exec_cmd('store_events', postcmd=False)
872 self.exec_cmd('reweight -from_cards', postcmd=False)
873@@ -3972,106 +3978,6 @@
874 self.Gdirs = Gdirs
875 return self.Gdirs
876
877- ############################################################################
878- def set_run_name(self, name, tag=None, level='parton', reload_card=False,
879- allow_new_tag=True):
880- """define the run name, the run_tag, the banner and the results."""
881-
882- # when are we force to change the tag new_run:previous run requiring changes
883- upgrade_tag = {'parton': ['parton','pythia','pgs','delphes'],
884- 'pythia': ['pythia','pgs','delphes'],
885- 'pgs': ['pgs'],
886- 'delphes':['delphes'],
887- 'plot':[],
888- 'syscalc':[]}
889-
890-
891-
892- if name == self.run_name:
893- if reload_card:
894- run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
895- self.run_card = banner_mod.RunCard(run_card)
896-
897- #check if we need to change the tag
898- if tag:
899- self.run_card['run_tag'] = tag
900- self.run_tag = tag
901- self.results.add_run(self.run_name, self.run_card)
902- else:
903- for tag in upgrade_tag[level]:
904- if getattr(self.results[self.run_name][-1], tag):
905- tag = self.get_available_tag()
906- self.run_card['run_tag'] = tag
907- self.run_tag = tag
908- self.results.add_run(self.run_name, self.run_card)
909- break
910- return # Nothing to do anymore
911-
912- # save/clean previous run
913- if self.run_name:
914- self.store_result()
915- # store new name
916- self.run_name = name
917-
918- new_tag = False
919- # First call for this run -> set the banner
920- self.banner = banner_mod.recover_banner(self.results, level, name)
921- if 'mgruncard' in self.banner:
922- self.run_card = self.banner.charge_card('run_card')
923- else:
924- # Read run_card
925- run_card = pjoin(self.me_dir, 'Cards','run_card.dat')
926- self.run_card = banner_mod.RunCard(run_card)
927-
928- if tag:
929- self.run_card['run_tag'] = tag
930- new_tag = True
931- elif not self.run_name in self.results and level =='parton':
932- pass # No results yet, so current tag is fine
933- elif not self.run_name in self.results:
934- #This is only for case when you want to trick the interface
935- logger.warning('Trying to run data on unknown run.')
936- self.results.add_run(name, self.run_card)
937- self.results.update('add run %s' % name, 'all', makehtml=False)
938- else:
939- for tag in upgrade_tag[level]:
940-
941- if getattr(self.results[self.run_name][-1], tag):
942- # LEVEL is already define in the last tag -> need to switch tag
943- tag = self.get_available_tag()
944- self.run_card['run_tag'] = tag
945- new_tag = True
946- break
947- if not new_tag:
948- # We can add the results to the current run
949- tag = self.results[self.run_name][-1]['tag']
950- self.run_card['run_tag'] = tag # ensure that run_tag is correct
951-
952- if allow_new_tag and (name in self.results and not new_tag):
953- self.results.def_current(self.run_name)
954- else:
955- self.results.add_run(self.run_name, self.run_card)
956-
957- self.run_tag = self.run_card['run_tag']
958-
959- # Return the tag of the previous run having the required data for this
960- # tag/run to working wel.
961- if level == 'parton':
962- return
963- elif level == 'pythia':
964- return self.results[self.run_name][0]['tag']
965- else:
966- for i in range(-1,-len(self.results[self.run_name])-1,-1):
967- tagRun = self.results[self.run_name][i]
968- if tagRun.pythia:
969- return tagRun['tag']
970-
971-
972-
973-
974-
975-
976-
977
978 ############################################################################
979 def find_model_name(self):
980@@ -4215,7 +4121,7 @@
981 scdir = self.options['syscalc_path']
982 if not scdir or not os.path.exists(scdir):
983 return
984- logger.info('running syscalc on mode %s' % mode)
985+ logger.info('running SysCalc on mode %s' % mode)
986
987
988 # Check that all pdfset are correctly installed
989@@ -4265,7 +4171,9 @@
990 if not (os.path.exists(event_path) or os.path.exists(event_path+".gz")):
991 event_path = pjoin(event_dir, 'unweighted_events.lhe')
992 output = pjoin(event_dir, 'syscalc.lhe')
993+ stdout = open(pjoin(event_dir, self.run_name, '%s_systematics.log' % (mode)),'w')
994 elif mode == 'Pythia':
995+ stdout = open(pjoin(event_dir, self.run_name, '%s_%s_syscalc.log' % (tag,mode)),'w')
996 if 'mgpythiacard' in self.banner:
997 pat = re.compile('''^\s*qcut\s*=\s*([\+\-\d.e]*)''', re.M+re.I)
998 data = pat.search(self.banner['mgpythiacard'])
999@@ -4294,7 +4202,7 @@
1000 try:
1001 proc = misc.call([os.path.join(scdir, 'sys_calc'),
1002 event_path, card, output],
1003- stdout = open(pjoin(event_dir, self.run_name, '%s_%s_syscalc.log' % (tag,mode)),'w'),
1004+ stdout = stdout,
1005 stderr = subprocess.STDOUT,
1006 cwd=event_dir)
1007 # Wait 5 s to make sure file is finished writing
1008
1009=== modified file 'madgraph/interface/madgraph_interface.py'
1010--- madgraph/interface/madgraph_interface.py 2016-08-22 20:43:28 +0000
1011+++ madgraph/interface/madgraph_interface.py 2016-08-24 22:15:21 +0000
1012@@ -192,9 +192,9 @@
1013 info['date'])
1014
1015 if os.path.exists(pjoin(MG5DIR, '.bzr')):
1016- proc = subprocess.Popen(['bzr', 'nick'], stdout=subprocess.PIPE)
1017+ proc = subprocess.Popen(['bzr', 'nick'], stdout=subprocess.PIPE,cwd=MG5DIR)
1018 bzrname,_ = proc.communicate()
1019- proc = subprocess.Popen(['bzr', 'revno'], stdout=subprocess.PIPE)
1020+ proc = subprocess.Popen(['bzr', 'revno'], stdout=subprocess.PIPE,cwd=MG5DIR)
1021 bzrversion,_ = proc.communicate()
1022 bzrname, bzrversion = bzrname.strip(), bzrversion.strip()
1023 len_name = len(bzrname)
1024
1025=== modified file 'madgraph/interface/reweight_interface.py'
1026--- madgraph/interface/reweight_interface.py 2016-08-18 22:31:05 +0000
1027+++ madgraph/interface/reweight_interface.py 2016-08-24 22:15:21 +0000
1028@@ -82,11 +82,11 @@
1029 self.model = None
1030 self.has_standalone_dir = False
1031 self.mother= mother # calling interface
1032-
1033+ self.multicore=False
1034
1035 self.options = {'curr_dir': os.path.realpath(os.getcwd()),
1036 'rwgt_name':None}
1037-
1038+
1039 self.events_file = None
1040 self.processes = {}
1041 self.second_model = None
1042@@ -201,7 +201,7 @@
1043 logger.info("options: %s" % option)
1044
1045
1046- def get_LO_definition_from_NLO(self, proc):
1047+ def get_LO_definition_from_NLO(self, proc, real_only=False):
1048 """return the LO definitions of the process corresponding to the born/real"""
1049
1050 # split the line definition with the part before and after the NLO tag
1051@@ -212,7 +212,10 @@
1052 #NO NLO tag => nothing to do actually return input
1053 return proc
1054 elif not order.startswith(('virt','loonly','noborn')):
1055- # OK this a standard NLO process
1056+ # OK this a standard NLO process
1057+ if real_only:
1058+ commandline= ''
1059+
1060 if '=' in order:
1061 # get the type NLO QCD/QED/...
1062 order = order.split('=',1)[1]
1063@@ -358,6 +361,10 @@
1064 self.rwgt_dir = args[1]
1065 if not os.path.exists(self.rwgt_dir):
1066 os.mkdir(self.rwgt_dir)
1067+ elif args[0] == 'multicore':
1068+ pass
1069+ # this line is meant to be parsed by common_run_interface and change the way this class is called.
1070+ #It has no direct impact on this class.
1071 else:
1072 logger.critical("unknown option! %s. Discard line." % args[0])
1073
1074@@ -415,12 +422,25 @@
1075
1076 model_line = self.banner.get('proc_card', 'full_model_line')
1077
1078- if not self.has_standalone_dir:
1079+ if not self.has_standalone_dir:
1080 if self.rwgt_dir and os.path.exists(pjoin(self.rwgt_dir,'rw_me','rwgt.pkl')):
1081 self.load_from_pickle()
1082- self.me_dir = self.rwgt_dir
1083+ if not self.rwgt_dir:
1084+ self.me_dir = self.rwgt_dir
1085+ elif self.multicore == 'wait':
1086+ while not os.path.exists(pjoin(self.me_dir,'rw_me','rwgt.pkl')):
1087+ time.sleep(60)
1088+ print 'wait for pickle'
1089+ print "loading from pickle"
1090+ if not self.rwgt_dir:
1091+ self.rwgt_dir = self.me_dir
1092+ self.load_from_pickle(keep_name=True)
1093 else:
1094- self.create_standalone_directory()
1095+ self.create_standalone_directory()
1096+ if self.multicore == 'create':
1097+ if not self.rwgt_dir:
1098+ self.rwgt_dir = self.me_dir
1099+ self.save_to_pickle()
1100
1101 if self.rwgt_dir:
1102 path_me =self.rwgt_dir
1103@@ -501,7 +521,6 @@
1104 rewgtid = maxid
1105 if self.options['rwgt_name']:
1106 #ensure that the entry is not already define if so overwrites it
1107- misc.sprint(mg_rwgt_info)
1108 for (i, nlotype, diff) in mg_rwgt_info[:]:
1109 for flag in type_rwgt:
1110 if 'rwgt_%s' % i == '%s%s' %(self.options['rwgt_name'],flag) or \
1111@@ -779,10 +798,11 @@
1112 results.current.parton.append('lhe')
1113 logger.info('Event %s is now created with new central weight' % output2.name)
1114
1115- for name in cross:
1116- if name == 'orig':
1117- continue
1118- logger.info('new cross-section is %s: %g pb (indicative error: %g pb)' %\
1119+ if self.multicore != 'create':
1120+ for name in cross:
1121+ if name == 'orig':
1122+ continue
1123+ logger.info('new cross-section is %s: %g pb (indicative error: %g pb)' %\
1124 ('(%s)' %name if name else '',cross[name], error[name]))
1125
1126 self.terminate_fortran_executables(new_card_only=True)
1127@@ -1023,6 +1043,11 @@
1128 def calculate_matrix_element(self, event, hypp_id, space, scale2=0):
1129 """routine to return the matrix element"""
1130
1131+ if self.has_nlo:
1132+ nb_retry, sleep = 10, 60
1133+ else:
1134+ nb_retry, sleep = 5, 20
1135+
1136 tag, order = event.get_tag_and_order()
1137 if isinstance(hypp_id, str) and hypp_id.startswith('V'):
1138 tag = (tag,'V')
1139@@ -1068,7 +1093,7 @@
1140 dir_to_f2py_free_mod[(Pdir,0)] = (metag, nb_f2py_module)
1141 os.environ['MENUM'] = '2'
1142 if not self.rwgt_dir or not os.path.exists(pjoin(Pdir, 'matrix2py.so')):
1143- misc.compile(['matrix2py.so'], cwd=Pdir)
1144+ misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix2py.so'], cwd=Pdir)
1145 self.rename_f2py_lib(Pdir, 2*metag)
1146 try:
1147 mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, 2*metag), globals(), locals(), [],-1)
1148@@ -1096,10 +1121,9 @@
1149 try:
1150 mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, newtag), globals(), locals(), [],-1)
1151 except Exception, error:
1152- os.remove(pjoin(Pdir, 'matrix%spy.so' % newtag))
1153 newtag = "L%s" % newtag
1154 os.environ['MENUM'] = newtag
1155- misc.compile(['matrix%spy.so' % newtag], cwd=Pdir)
1156+ misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix%spy.so' % newtag], cwd=Pdir)
1157 mymod = __import__('%s.SubProcesses.%s.matrix%spy' % (base, Pname, newtag), globals(), locals(), [],-1)
1158
1159 S = mymod.SubProcesses
1160@@ -1121,7 +1145,7 @@
1161 Pname = os.path.basename(Pdir)
1162 os.environ['MENUM'] = '2'
1163 if not self.rwgt_dir or not os.path.exists(pjoin(Pdir, 'matrix2py.so')):
1164- misc.compile(['matrix2py.so'], cwd=pjoin(subdir, Pdir))
1165+ misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix2py.so'], cwd=pjoin(subdir, Pdir))
1166 if (Pdir, 1) not in dir_to_f2py_free_mod:
1167 metag = 1
1168 dir_to_f2py_free_mod[(Pdir,1)] = (metag, nb_f2py_module)
1169@@ -1134,10 +1158,9 @@
1170 try:
1171 mymod = __import__("%s_second.SubProcesses.%s.matrix%spy" % (base, Pname, metag))
1172 except ImportError:
1173- os.remove(pjoin(Pdir, 'matrix%spy.so' % metag ))
1174 metag = "L%s" % metag
1175 os.environ['MENUM'] = str(metag)
1176- misc.compile(['matrix%spy.so' % metag], cwd=pjoin(subdir, Pdir))
1177+ misc.multiple_try(nb_retry,sleep)(misc.compile)(['matrix%spy.so' % metag], cwd=pjoin(subdir, Pdir))
1178 mymod = __import__("%s_second.SubProcesses.%s.matrix%spy" % (base, Pname, metag))
1179
1180 reload(mymod)
1181@@ -1211,21 +1234,24 @@
1182 split = line.split()
1183 if len(split) == 4:
1184 cross, error = float(split[0]), float(split[1])
1185- if 'orig' not in self.all_cross_section:
1186- logger.info('Original cross-section: %s +- %s pb' % (cross, error))
1187- else:
1188- logger.info('Original cross-section: %s +- %s pb (cross-section from sum of weights: %s)' % (cross, error, self.all_cross_section['orig'][0]))
1189- logger.info('Computed cross-section:')
1190- keys = self.all_cross_section.keys()
1191- keys.sort()
1192- for key in keys:
1193- if key == 'orig':
1194- continue
1195- logger.info('%s : %s +- %s pb' % (key[0] if not key[1] else '%s%s' % key,
1196- self.all_cross_section[key][0],self.all_cross_section[key][1] ))
1197+
1198+ if not self.multicore == 'create':
1199+ # No print of results for the multicore mode for the one printed on screen
1200+ if 'orig' not in self.all_cross_section:
1201+ logger.info('Original cross-section: %s +- %s pb' % (cross, error))
1202+ else:
1203+ logger.info('Original cross-section: %s +- %s pb (cross-section from sum of weights: %s)' % (cross, error, self.all_cross_section['orig'][0]))
1204+ logger.info('Computed cross-section:')
1205+ keys = self.all_cross_section.keys()
1206+ keys.sort()
1207+ for key in keys:
1208+ if key == 'orig':
1209+ continue
1210+ logger.info('%s : %s +- %s pb' % (key[0] if not key[1] else '%s%s' % key,
1211+ self.all_cross_section[key][0],self.all_cross_section[key][1] ))
1212 self.terminate_fortran_executables()
1213
1214- if self.rwgt_dir:
1215+ if self.rwgt_dir and self.multicore == False:
1216 self.save_to_pickle()
1217
1218 with misc.stdchannel_redirected(sys.stdout, os.devnull):
1219@@ -1245,7 +1271,7 @@
1220 @misc.mute_logger()
1221 def create_standalone_directory(self, second=False):
1222 """generate the various directory for the weight evaluation"""
1223-
1224+
1225 data={}
1226 if not second:
1227 data['paths'] = ['rw_me', 'rw_mevirt']
1228@@ -1350,12 +1376,18 @@
1229 logger.info('generating the square matrix element for reweighting (second model and/or processes)')
1230 start = time.time()
1231 commandline=''
1232- for proc in data['processes']:
1233+ for i,proc in enumerate(data['processes']):
1234 if '[' not in proc:
1235 commandline += "add process %s ;" % proc
1236 else:
1237 has_nlo = True
1238- commandline += self.get_LO_definition_from_NLO(proc)
1239+ if self.banner.get('run_card','ickkw') == 3:
1240+ if len(proc) == min([len(p.strip()) for p in data['processes']]):
1241+ commandline += self.get_LO_definition_from_NLO(proc)
1242+ else:
1243+ commandline += self.get_LO_definition_from_NLO(proc, real_only=True)
1244+ else:
1245+ commandline += self.get_LO_definition_from_NLO(proc)
1246
1247 commandline = commandline.replace('add process', 'generate',1)
1248 logger.info(commandline)
1249@@ -1434,6 +1466,18 @@
1250 pjoin(path_me, data['paths'][0], 'Cards', 'MadLoopParams.dat'),
1251 commentdefault=False)
1252
1253+ if self.multicore == 'create':
1254+ print "compile OLP", data['paths'][0]
1255+ misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][0],'SubProcesses'),
1256+ nb_core=self.mother.options['nb_core'])
1257+
1258+ if os.path.exists(pjoin(path_me, data['paths'][1], 'Cards', 'MadLoopParams.dat')):
1259+ if self.multicore == 'create':
1260+ print "compile OLP", data['paths'][1]
1261+ misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'),
1262+ nb_core=self.mother.options['nb_core'])
1263+
1264+
1265 # 5. create the virtual for NLO reweighting ---------------------------
1266 if has_nlo and 'NLO' in self.rwgt_mode:
1267 # Do not pass here for LO/NLO_tree
1268@@ -1527,6 +1571,11 @@
1269 self.banner.run_card['lpp2']],
1270 self.banner.run_card.get_lhapdf_id())
1271 self.combine_wgt = mymod.get_wgt
1272+
1273+ if self.multicore == 'create':
1274+ print "compile OLP", data['paths'][1]
1275+ misc.compile(['OLP_static'], cwd=pjoin(path_me, data['paths'][1],'SubProcesses'),
1276+ nb_core=self.mother.options['nb_core'])
1277
1278 elif has_nlo and not second and self.rwgt_mode == ['NLO_tree']:
1279 # We do not have any virtual reweighting to do but we still have to
1280@@ -1630,18 +1679,22 @@
1281 to_save['rwgt_dir'] = self.rwgt_dir
1282 to_save['has_nlo'] = self.has_nlo
1283 to_save['rwgt_mode'] = self.rwgt_mode
1284+ to_save['rwgt_name'] = self.options['rwgt_name']
1285
1286 name = pjoin(self.rwgt_dir, 'rw_me', 'rwgt.pkl')
1287 save_load_object.save_to_file(name, to_save)
1288
1289- def load_from_pickle(self):
1290+
1291+ def load_from_pickle(self, keep_name=False):
1292 import madgraph.iolibs.save_load_object as save_load_object
1293
1294 obj = save_load_object.load_from_file( pjoin(self.rwgt_dir, 'rw_me', 'rwgt.pkl'))
1295
1296 self.has_standalone_dir = True
1297 self.options = {'curr_dir': os.path.realpath(os.getcwd()),
1298- 'rewgt_name': None}
1299+ 'rwgt_name': None}
1300+ if keep_name:
1301+ self.options['rwgt_name'] = obj['rwgt_name']
1302
1303 old_rwgt = obj['rwgt_dir']
1304
1305@@ -1665,10 +1718,14 @@
1306 if not self.rwgt_mode:
1307 self.rwgt_mode = obj['rwgt_mode']
1308 logger.info("mode set to %s" % self.rwgt_mode)
1309- if self.has_nlo:
1310+ if self.has_nlo and 'NLO' in self.rwgt_mode:
1311 path = pjoin(obj['rwgt_dir'], 'rw_mevirt','Source')
1312 sys.path.insert(0, path)
1313- mymod = __import__('rwgt2py', globals(), locals())
1314+ try:
1315+ mymod = __import__('rwgt2py', globals(), locals())
1316+ except ImportError:
1317+ misc.compile(['rwgt2py.so'], cwd=path)
1318+ mymod = __import__('rwgt2py', globals(), locals())
1319 with misc.stdchannel_redirected(sys.stdout, os.devnull):
1320 mymod.initialise([self.banner.run_card['lpp1'],
1321 self.banner.run_card['lpp2']],
1322
1323=== modified file 'madgraph/iolibs/export_fks.py'
1324--- madgraph/iolibs/export_fks.py 2016-08-07 07:16:07 +0000
1325+++ madgraph/iolibs/export_fks.py 2016-08-24 22:15:21 +0000
1326@@ -268,7 +268,8 @@
1327 pjoin('various','FO_analyse_card.py'),
1328 pjoin('various','histograms.py'),
1329 pjoin('various','banner.py'),
1330- pjoin('various','cluster.py'),
1331+ pjoin('various','cluster.py'),
1332+ pjoin('various','systematics.py'),
1333 pjoin('various','lhe_parser.py'),
1334 pjoin('madevent','sum_html.py'),
1335 pjoin('madevent','gen_crossxhtml.py'),
1336
1337=== modified file 'madgraph/iolibs/export_v4.py'
1338--- madgraph/iolibs/export_v4.py 2016-08-15 14:06:30 +0000
1339+++ madgraph/iolibs/export_v4.py 2016-08-24 22:15:21 +0000
1340@@ -3227,6 +3227,8 @@
1341 self.dir_path+'/bin/internal/lhe_parser.py')
1342 cp(_file_path+'/various/banner.py',
1343 self.dir_path+'/bin/internal/banner.py')
1344+ cp(_file_path+'/various/systematics.py', self.dir_path+'/bin/internal/systematics.py')
1345+
1346 cp(_file_path+'/various/cluster.py',
1347 self.dir_path+'/bin/internal/cluster.py')
1348 cp(_file_path+'/madevent/combine_runs.py',
1349
1350=== modified file 'madgraph/madevent/gen_crossxhtml.py'
1351--- madgraph/madevent/gen_crossxhtml.py 2016-04-15 11:14:25 +0000
1352+++ madgraph/madevent/gen_crossxhtml.py 2016-08-24 22:15:21 +0000
1353@@ -812,7 +812,7 @@
1354 self.parton.append('param_card')
1355
1356 if 'syst' not in self.parton and \
1357- exists(pjoin(path, "%s_parton_syscalc.log" %self['tag'])):
1358+ exists(pjoin(path, "parton_systematics.log")):
1359 self.parton.append('syst')
1360
1361 for kind in ['top','HwU','pdf','ps']:
1362@@ -1150,8 +1150,8 @@
1363 local_dico['err'] = self['error']
1364 local_dico['nb_event'] = self['nb_event']
1365 if 'syst' in self.parton:
1366- local_dico['syst'] = '<font face=symbol>&#177;</font> <a href="./Events/%(run_name)s/%(tag)s_parton_syscalc.log">systematics</a>' \
1367- % {'run_name':self['run_name'], 'tag': self['tag']}
1368+ local_dico['syst'] = '<font face=symbol>&#177;</font> <a href="./Events/%(run_name)s/parton_systematics.log">systematics</a>' \
1369+ % {'run_name':self['run_name']}
1370 elif self['cross_pythia']:
1371 if self.parton:
1372 local_dico['cross_span'] = nb_line -1
1373@@ -1173,7 +1173,7 @@
1374 local_dico['err'] = self['error']
1375 local_dico['nb_event'] = self['nb_event']
1376 if 'syst' in self.parton:
1377- local_dico['syst'] = '<font face=symbol>&#177;</font> <a href="./Events/%(run_name)s/%(tag)s_parton_syscalc.log">systematics</a>' \
1378+ local_dico['syst'] = '<font face=symbol>&#177;</font> <a href="./Events/%(run_name)s/parton_systematics.log">systematics</a>' \
1379 % {'run_name':self['run_name'], 'tag': self['tag']}
1380
1381 elif ttype == 'pythia' and self['cross_pythia']:
1382
1383=== modified file 'madgraph/various/banner.py'
1384--- madgraph/various/banner.py 2016-08-17 19:36:40 +0000
1385+++ madgraph/various/banner.py 2016-08-24 22:15:21 +0000
1386@@ -193,34 +193,6 @@
1387 return cross, math.sqrt(error)
1388
1389
1390-
1391- def modify_init_cross(self, cross):
1392- """modify the init information with the associate cross-section"""
1393-
1394- assert isinstance(cross, dict)
1395-# assert "all" in cross
1396- assert "init" in self
1397-
1398- all_lines = self["init"].split('\n')
1399- new_data = []
1400- new_data.append(all_lines[0])
1401- for i in range(1, len(all_lines)):
1402- line = all_lines[i]
1403- split = line.split()
1404- if len(split) == 4:
1405- xsec, xerr, xmax, pid = split
1406- else:
1407- new_data += all_lines[i:]
1408- break
1409- if int(pid) not in cross:
1410- raise Exception
1411- pid = int(pid)
1412- ratio = cross[pid]/float(xsec)
1413- line = " %+13.7e %+13.7e %+13.7e %i" % \
1414- (float(cross[pid]), ratio* float(xerr), ratio*float(xmax), pid)
1415- new_data.append(line)
1416- self['init'] = '\n'.join(new_data)
1417-
1418 def scale_init_cross(self, ratio):
1419 """modify the init information with the associate scale"""
1420
1421@@ -244,6 +216,15 @@
1422 new_data.append(line)
1423 self['init'] = '\n'.join(new_data)
1424
1425+ def get_pdg_beam(self):
1426+ """return the pdg of each beam"""
1427+
1428+ assert "init" in self
1429+
1430+ all_lines = self["init"].split('\n')
1431+ pdg1,pdg2,_ = all_lines[0].split(None, 2)
1432+ return int(pdg1), int(pdg2)
1433+
1434 def load_basic(self, medir):
1435 """ Load the proc_card /param_card and run_card """
1436
1437
1438=== modified file 'madgraph/various/cluster.py'
1439--- madgraph/various/cluster.py 2016-06-13 08:27:43 +0000
1440+++ madgraph/various/cluster.py 2016-08-24 22:15:21 +0000
1441@@ -618,6 +618,8 @@
1442 if isinstance(exe,str):
1443 if os.path.exists(exe) and not exe.startswith('/'):
1444 exe = './' + exe
1445+ if isinstance(opt['stdout'],str):
1446+ opt['stdout'] = open(opt['stdout'],'w')
1447 if opt['stderr'] == None:
1448 opt['stderr'] = subprocess.STDOUT
1449 proc = misc.Popen([exe] + arg, **opt)
1450@@ -673,7 +675,6 @@
1451
1452 tag = (prog, tuple(argument), cwd, nb_submit)
1453 if isinstance(prog, str):
1454-
1455
1456 opt = {'cwd': cwd,
1457 'stdout':stdout,
1458
1459=== modified file 'madgraph/various/lhe_parser.py'
1460--- madgraph/various/lhe_parser.py 2016-07-29 15:27:30 +0000
1461+++ madgraph/various/lhe_parser.py 2016-08-24 22:15:21 +0000
1462@@ -178,6 +178,8 @@
1463 if '.gz' in path and isinstance(self, EventFileNoGzip) and\
1464 mode == 'r' and os.path.exists(path[:-3]):
1465 super(EventFile, self).__init__(path[:-3], mode, *args, **opt)
1466+ else:
1467+ raise
1468
1469 self.banner = ''
1470 if mode == 'r':
1471@@ -506,7 +508,26 @@
1472 else:
1473 return out
1474
1475+ def split(self, nb_event=0):
1476+ """split the file in multiple file. Do not change the weight!"""
1477
1478+ nb_file = -1
1479+ for i, event in enumerate(self):
1480+ if i % nb_event == 0:
1481+ if i:
1482+ #close previous file
1483+ current.write('</LesHouchesEvent>\n')
1484+ current.close()
1485+ # create the new file
1486+ nb_file +=1
1487+ current = open('%s_%s.lhe' % (self.name, nb_file),'w')
1488+ current.write(self.banner)
1489+ current.write(str(event))
1490+ if i!=0:
1491+ current.write('</LesHouchesEvent>\n')
1492+ current.close()
1493+ return nb_file +1
1494+
1495 def update_Hwu(self, hwu, fct, name='lhe', keep_wgt=True):
1496
1497 first=True
1498@@ -620,13 +641,14 @@
1499 The number of events in each file need to be provide in advance
1500 (if not provide the file is first read to find that number"""
1501
1502- def __new__(cls, start_list=[]):
1503+ def __new__(cls, start_list=[],parse=True):
1504 return object.__new__(MultiEventFile)
1505
1506- def __init__(self, start_list=[]):
1507+ def __init__(self, start_list=[], parse=True):
1508 """if trunc_error is define here then this allow
1509 to only read all the files twice and not three times."""
1510 self.files = []
1511+ self.parsefile = parse #if self.files is formatted or just the path
1512 self.banner = ''
1513 self.initial_nb_events = []
1514 self.total_event_in_files = 0
1515@@ -636,8 +658,11 @@
1516 self.across = []
1517 self.scales = []
1518 if start_list:
1519- for p in start_list:
1520- self.add(p)
1521+ if parse:
1522+ for p in start_list:
1523+ self.add(p)
1524+ else:
1525+ self.files = start_list
1526 self._configure = False
1527
1528 def add(self, path, cross, error, across):
1529@@ -881,6 +906,66 @@
1530
1531 return super(MultiEventFile, self).unweight(outputpath, get_wgt_multi, **opts)
1532
1533+ def write(self, path, random=False, banner=None, get_info=False):
1534+ """ """
1535+
1536+ if isinstance(path, str):
1537+ out = EventFile(path, 'w')
1538+ if self.parsefile and not banner:
1539+ banner = self.files[0].banner
1540+ elif not banner:
1541+ firstlhe = EventFile(self.files[0])
1542+ banner = firstlhe.banner
1543+ else:
1544+ out = path
1545+ if banner:
1546+ out.write(banner)
1547+ nb_event = 0
1548+ info = collections.defaultdict(float)
1549+ if random and self.open:
1550+ for event in self:
1551+ nb_event +=1
1552+ out.write(event)
1553+ if get_info:
1554+ event.parse_reweight()
1555+ for key, value in event.reweight_data.items():
1556+ info[key] += value
1557+ info['central'] += event.wgt
1558+ elif not random:
1559+ for i,f in enumerate(self.files):
1560+ #check if we need to parse the file or not
1561+ if not self.parsefile:
1562+ if i==0:
1563+ try:
1564+ lhe = firstlhe
1565+ except:
1566+ lhe = EventFile(f)
1567+ else:
1568+ lhe = EventFile(f)
1569+ else:
1570+ lhe = f
1571+ for event in lhe:
1572+ nb_event +=1
1573+ if get_info:
1574+ event.parse_reweight()
1575+ for key, value in event.reweight_data.items():
1576+ info[key] += value
1577+ info['central'] += event.wgt
1578+ out.write(str(event))
1579+ lhe.close()
1580+ out.write("</LesHouchesEvents>\n")
1581+ return nb_event, info
1582+
1583+ def remove(self):
1584+ """ """
1585+ if self.parsefile:
1586+ for f in self.files:
1587+ os.remove(f.name)
1588+ else:
1589+ for f in self.files:
1590+ os.remove(f)
1591+
1592+
1593
1594 class Event(list):
1595 """Class storing a single event information (list of particles + global information)"""
1596@@ -1025,7 +1110,7 @@
1597 self.reweight_order = []
1598 start, stop = self.tag.find('<rwgt>'), self.tag.find('</rwgt>')
1599 if start != -1 != stop :
1600- pattern = re.compile(r'''<\s*wgt id=(?:\'|\")(?P<id>[^\'\"]+)(?:\'|\")\s*>\s*(?P<val>[\ded+-.]*)\s*</wgt>''')
1601+ pattern = re.compile(r'''<\s*wgt id=(?:\'|\")(?P<id>[^\'\"]+)(?:\'|\")\s*>\s*(?P<val>[\ded+-.]*)\s*</wgt>''',re.I)
1602 data = pattern.findall(self.tag)
1603 try:
1604 self.reweight_data = dict([(pid, float(value)) for (pid, value) in data
1605@@ -1046,7 +1131,54 @@
1606
1607 text = self.tag[start+8:stop]
1608 self.nloweight = NLO_PARTIALWEIGHT(text, self)
1609-
1610+ return self.nloweight
1611+
1612+ def parse_lo_weight(self):
1613+ """ """
1614+ if hasattr(self, 'loweight'):
1615+ return self.loweight
1616+
1617+ start, stop = self.tag.find('<mgrwt>'), self.tag.find('</mgrwt>')
1618+
1619+ if start != -1 != stop :
1620+ text = self.tag[start+8:stop]
1621+#<rscale> 3 0.29765919e+03</rscale>
1622+#<asrwt>0</asrwt>
1623+#<pdfrwt beam="1"> 1 21 0.15134321e+00 0.29765919e+03</pdfrwt>
1624+#<pdfrwt beam="2"> 1 21 0.38683649e-01 0.29765919e+03</pdfrwt>
1625+#<totfact> 0.17315115e+03</totfact>
1626+ self.loweight={}
1627+ for line in text.split('\n'):
1628+ line = line.replace('<', ' <').replace("'",'"')
1629+ if 'rscale' in line:
1630+ _, nqcd, scale, _ = line.split()
1631+ self.loweight['n_qcd'] = int(nqcd)
1632+ self.loweight['ren_scale'] = float(scale)
1633+ elif '<pdfrwt beam="1"' in line:
1634+ args = line.split()
1635+ self.loweight['n_pdfrw1'] = int(args[2])
1636+ npdf = self.loweight['n_pdfrw1']
1637+ self.loweight['pdf_pdg_code1'] = [int(i) for i in args[3:3+npdf]]
1638+ self.loweight['pdf_x1'] = [float(i) for i in args[3+npdf:3+2*npdf]]
1639+ self.loweight['pdf_q1'] = [float(i) for i in args[3+2*npdf:3+3*npdf]]
1640+ elif '<pdfrwt beam="2"' in line:
1641+ args = line.split()
1642+ self.loweight['n_pdfrw2'] = int(args[2])
1643+ npdf = self.loweight['n_pdfrw2']
1644+ self.loweight['pdf_pdg_code2'] = [int(i) for i in args[3:3+npdf]]
1645+ self.loweight['pdf_x2'] = [float(i) for i in args[3+npdf:3+2*npdf]]
1646+ self.loweight['pdf_q2'] = [float(i) for i in args[3+2*npdf:3+3*npdf]]
1647+ elif '<asrwt>' in line:
1648+ args = line.replace('>','> ').split()
1649+ nalps = int(args[1])
1650+ self.loweight['asrwt'] = [float(a) for a in args[2:2+nalps]]
1651+
1652+ elif 'totfact' in line:
1653+ args = line.split()
1654+ self.loweight['tot_fact'] = float(args[1])
1655+ else:
1656+ return None
1657+ return self.loweight
1658
1659 def parse_matching_scale(self):
1660 """Parse the line containing the starting scale for the shower"""
1661@@ -1608,10 +1740,40 @@
1662 for particle in self:
1663 if particle.status != 1:
1664 continue
1665- scale += particle.mass**2 + particle.momentum.pt**2
1666+ p = FourMomentum(particle)
1667+ scale += math.sqrt(p.mass_sqr + p.pt**2)
1668
1669 return prefactor * scale
1670
1671+
1672+ def get_et_scale(self, prefactor=1):
1673+
1674+ scale = 0
1675+ for particle in self:
1676+ if particle.status != 1:
1677+ continue
1678+ p = FourMomentum(particle)
1679+ pt = p.pt
1680+ if (pt>0):
1681+ scale += p.E*pt/math.sqrt(pt**2+p.pz**2)
1682+
1683+ return prefactor * scale
1684+
1685+ def get_sqrts_scale(self, prefactor=1):
1686+
1687+ scale = 0
1688+ init = []
1689+ for particle in self:
1690+ if particle.status == -1:
1691+ init.append(FourMomentum(particle))
1692+ if len(init) == 1:
1693+ return init[0].mass
1694+ elif len(init)==2:
1695+ return math.sqrt((init[0]+init[1])**2)
1696+
1697+
1698+
1699+
1700 def get_momenta_str(self, get_order, allow_reversed=True):
1701 """return the momenta str in the order asked for"""
1702
1703@@ -1785,7 +1947,7 @@
1704 if power == 1:
1705 return FourMomentum(self)
1706 elif power == 2:
1707- return self.mass_sqr()
1708+ return self.mass_sqr
1709
1710 def __repr__(self):
1711 return 'FourMomentum(%s,%s,%s,%s)' % (self.E, self.px, self.py,self.pz)
1712@@ -1850,6 +2012,20 @@
1713 if isinstance(input, str):
1714 self.parse(input)
1715
1716+ def __str__(self):
1717+
1718+ out = """ pwgt: %(pwgt)s
1719+ born, real : %(born)s %(real)s
1720+ pdgs : %(pdgs)s
1721+ bjks : %(bjks)s
1722+ scales**2, gs: %(scales2)s %(gs)s
1723+ born/real related : %(born_related)s %(real_related)s
1724+ type / nfks : %(type)s %(nfks)s
1725+ to merge : %(to_merge_pdg)s in %(merge_new_pdg)s
1726+ ref_wgt : %(ref_wgt)s""" % self.__dict__
1727+ return out
1728+
1729+
1730 def parse(self, text):
1731 """parse the line and create the related object"""
1732 #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
1733@@ -2086,6 +2262,42 @@
1734 @property
1735 def aqcd(self):
1736 return self.event.aqcd
1737+
1738+ def get_ht_scale(self, prefactor=1):
1739+
1740+ scale = 0
1741+ for particle in self:
1742+ p = particle
1743+ scale += math.sqrt(max(0, p.mass_sqr + p.pt**2))
1744+
1745+ return prefactor * scale
1746+
1747+ def get_et_scale(self, prefactor=1):
1748+
1749+ scale = 0
1750+ for particle in self:
1751+ p = particle
1752+ pt = p.pt
1753+ if (pt>0):
1754+ scale += p.E*pt/math.sqrt(pt**2+p.pz**2)
1755+
1756+ return prefactor * scale
1757+
1758+
1759+ def get_sqrts_scale(self, event,prefactor=1):
1760+
1761+ scale = 0
1762+ nb_init = 0
1763+ for particle in event:
1764+ if particle.status == -1:
1765+ nb_init+=1
1766+ if nb_init == 1:
1767+ return self[0].mass
1768+ elif nb_init==2:
1769+ return math.sqrt((self[0]+self[1])**2)
1770+
1771+
1772+
1773
1774 def __init__(self, input, event):
1775
1776
1777=== modified file 'madgraph/various/misc.py'
1778--- madgraph/various/misc.py 2016-08-07 07:16:07 +0000
1779+++ madgraph/various/misc.py 2016-08-24 22:15:21 +0000
1780@@ -1334,6 +1334,7 @@
1781 if not __debug__:
1782 return
1783
1784+ use_print = False
1785 import inspect
1786 if opt.has_key('cond') and not opt['cond']:
1787 return
1788@@ -1346,6 +1347,8 @@
1789 level = opt['level']
1790 else:
1791 level = logging.getLogger('madgraph').level
1792+ if level == 0:
1793+ use_print = True
1794 #print "madgraph level",level
1795 #if level == 20:
1796 # level = 10 #avoid info level
1797@@ -1378,12 +1381,18 @@
1798 else:
1799 intro = ''
1800
1801-
1802- log.log(level, ' '.join([intro]+[str(a) for a in args]) + \
1803+ if not use_print:
1804+ log.log(level, ' '.join([intro]+[str(a) for a in args]) + \
1805 ' \033[1;30m[%s at line %s]\033[0m' % (os.path.basename(filename), lineno))
1806+<<<<<<< TREE
1807
1808 if wait:
1809 raw_input('press_enter to continue')
1810+=======
1811+ else:
1812+ print ' '.join([intro]+[str(a) for a in args]) + \
1813+ ' \033[1;30m[%s at line %s]\033[0m' % (os.path.basename(filename), lineno)
1814+>>>>>>> MERGE-SOURCE
1815 return
1816
1817 ################################################################################
1818@@ -1606,6 +1615,7 @@
1819 except:
1820 def apple_notify(subtitle, info_text, userInfo={}):
1821 return
1822+<<<<<<< TREE
1823 ## End apple notify
1824
1825
1826@@ -1672,3 +1682,76 @@
1827
1828
1829
1830+=======
1831+
1832+
1833+python_lhapdf=None
1834+def import_python_lhapdf(lhapdfconfig):
1835+ """load the python module of lhapdf return None if it can not be loaded"""
1836+
1837+ #save the result to have it faster and avoid the segfault at the second try if lhapdf is not compatible
1838+ global python_lhapdf
1839+ if python_lhapdf:
1840+ if python_lhapdf == -1:
1841+ return None
1842+ else:
1843+ return python_lhapdf
1844+
1845+ use_lhapdf=False
1846+ try:
1847+ lhapdf_libdir=subprocess.Popen([lhapdfconfig,'--libdir'],\
1848+ stdout=subprocess.PIPE).stdout.read().strip()
1849+ except:
1850+ use_lhapdf=False
1851+ return False
1852+ else:
1853+ try:
1854+ candidates=[dirname for dirname in os.listdir(lhapdf_libdir) \
1855+ if os.path.isdir(os.path.join(lhapdf_libdir,dirname))]
1856+ except OSError:
1857+ candidates=[]
1858+ for candidate in candidates:
1859+ if os.path.isfile(os.path.join(lhapdf_libdir,candidate,'site-packages','lhapdf.so')):
1860+ sys.path.insert(0,os.path.join(lhapdf_libdir,candidate,'site-packages'))
1861+ try:
1862+ import lhapdf
1863+ use_lhapdf=True
1864+ break
1865+ except ImportError:
1866+ sys.path.pop(0)
1867+ continue
1868+ if not use_lhapdf:
1869+ try:
1870+ candidates=[dirname for dirname in os.listdir(lhapdf_libdir+'64') \
1871+ if os.path.isdir(os.path.join(lhapdf_libdir+'64',dirname))]
1872+ except OSError:
1873+ candidates=[]
1874+ for candidate in candidates:
1875+ if os.path.isfile(os.path.join(lhapdf_libdir+'64',candidate,'site-packages','lhapdf.so')):
1876+ sys.path.insert(0,os.path.join(lhapdf_libdir+'64',candidate,'site-packages'))
1877+ try:
1878+ import lhapdf
1879+ use_lhapdf=True
1880+ break
1881+ except ImportError:
1882+ sys.path.pop(0)
1883+ continue
1884+ if not use_lhapdf:
1885+ try:
1886+ import lhapdf
1887+ use_lhapdf=True
1888+ except ImportError:
1889+ print 'fail'
1890+ logger.warning("Failed to access python version of LHAPDF: "\
1891+ "If the python interface to LHAPDF is available on your system, try "\
1892+ "adding its location to the PYTHONPATH environment variable and the"\
1893+ "LHAPDF library location to LD_LIBRARY_PATH (linux) or DYLD_LIBRARY_PATH (mac os x).")
1894+
1895+ if use_lhapdf:
1896+ python_lhapdf = lhapdf
1897+ python_lhapdf.setVerbosity(0)
1898+ else:
1899+ python_lhapdf = None
1900+
1901+ return python_lhapdf
1902+>>>>>>> MERGE-SOURCE
1903
1904=== added file 'madgraph/various/systematics.py'
1905--- madgraph/various/systematics.py 1970-01-01 00:00:00 +0000
1906+++ madgraph/various/systematics.py 2016-08-24 22:15:21 +0000
1907@@ -0,0 +1,724 @@
1908+################################################################################
1909+#
1910+# Copyright (c) 2016 The MadGraph5_aMC@NLO Development team and Contributors
1911+#
1912+# This file is a part of the MadGraph5_aMC@NLO project, an application which
1913+# automatically generates Feynman diagrams and matrix elements for arbitrary
1914+# high-energy processes in the Standard Model and beyond.
1915+#
1916+# It is subject to the MadGraph5_aMC@NLO license which should accompany this
1917+# distribution.
1918+#
1919+# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
1920+#
1921+################################################################################
1922+from __future__ import division
1923+if __name__ == "__main__":
1924+ import sys
1925+ import os
1926+ root = os.path.dirname(__file__)
1927+ if os.path.basename(root) == 'internal':
1928+ sys.path.append(os.path.dirname(root))
1929+ else:
1930+ sys.path.append(os.path.dirname(os.path.dirname(root)))
1931+
1932+import lhe_parser
1933+import banner
1934+import banner as banner_mod
1935+import itertools
1936+import misc
1937+import math
1938+import os
1939+import re
1940+import sys
1941+import time
1942+import StringIO
1943+
1944+pjoin = os.path.join
1945+
1946+class SystematicsError(Exception):
1947+ pass
1948+
1949+class Systematics(object):
1950+
1951+ def __init__(self, input_file, output_file,
1952+ start_event=0, stop_event=sys.maxint, write_banner=False,
1953+ mur=[0.5,1,2],
1954+ muf=[0.5,1,2],
1955+ alps=[1],
1956+ pdf='errorset', #[(id, subset)]
1957+ dyn=[-1,1,2,3,4],
1958+ together=[('mur', 'muf', 'dyn')],
1959+ lhapdf_config=misc.which('lhapdf-config'),
1960+ log=lambda x: sys.stdout.write(str(x)+'\n')
1961+ ):
1962+
1963+ # INPUT/OUTPUT FILE
1964+ if isinstance(input_file, str):
1965+ self.input = lhe_parser.EventFile(input_file)
1966+ else:
1967+ self.input = input_file
1968+ self.output_path = output_file
1969+ if output_file != None:
1970+ if isinstance(output_file, str):
1971+ if output_file == input_file:
1972+ directory,name = os.path.split(output_file)
1973+ new_name = pjoin(directory, '.tmp_'+name)
1974+ self.output = lhe_parser.EventFile(new_name, 'w')
1975+ else:
1976+ self.output = lhe_parser.EventFile(output_file, 'w')
1977+ else:
1978+ self.output = output_file
1979+ self.log = log
1980+
1981+ #get some information from the run_card.
1982+ self.banner = banner_mod.Banner(self.input.banner)
1983+ self.force_write_banner = bool(write_banner)
1984+ self.orig_dyn = self.banner.get('run_card', 'dynamical_scale_choice')
1985+ self.orig_pdf = self.banner.run_card.get_lhapdf_id()
1986+
1987+ #check for beam
1988+ beam1, beam2 = self.banner.get_pdg_beam()
1989+ if abs(beam1) != 2212 and abs(beam2) != 2212:
1990+ raise SystematicsError, 'can only reweight proton beam'
1991+ elif abs(beam1) != 2212:
1992+ self.b1 = 0
1993+ self.b2 = beam2//2212
1994+ elif abs(beam2) != 2212:
1995+ self.b1 = beam1//2212
1996+ self.b2 = 0
1997+ else:
1998+ self.b1 = beam1//2212
1999+ self.b2 = beam2//2212
2000+
2001+ if isinstance(self.banner.run_card, banner_mod.RunCardLO):
2002+ self.is_lo = True
2003+ if not self.banner.run_card['use_syst']:
2004+ raise SystematicsError, 'The events was not generated with use_syst=True. Can not evaluate systematics error on this event.'
2005+ else:
2006+ self.is_lo = False
2007+ if not self.banner.run_card['store_rwgt_info']:
2008+ raise SystematicsError, 'The events was not generated with store_rwgt_info=True. Can not evaluate systematics error on this event.'
2009+
2010+ # MUR/MUF/ALPS PARSING
2011+ if isinstance(mur, str):
2012+ mur = mur.split(',')
2013+ self.mur=[float(i) for i in mur]
2014+ if isinstance(muf, str):
2015+ muf = muf.split(',')
2016+ self.muf=[float(i) for i in muf]
2017+
2018+ if isinstance(alps, str):
2019+ alps = alps.split(',')
2020+ self.alps=[float(i) for i in alps]
2021+
2022+ # DYNAMICAL SCALE PARSING + together
2023+ if isinstance(dyn, str):
2024+ dyn = dyn.split(',')
2025+ self.dyn=[int(i) for i in dyn]
2026+
2027+ if isinstance(together, str):
2028+ self.together = together.split(',')
2029+ else:
2030+ self.together = together
2031+
2032+ # START/STOP EVENT
2033+ self.start_event=int(start_event)
2034+ self.stop_event=int(stop_event)
2035+ if start_event != 0:
2036+ self.log( "#starting from event #%s" % start_event)
2037+ if stop_event != sys.maxint:
2038+ self.log( "#stopping at event #%s" % stop_event)
2039+
2040+ # LHAPDF set
2041+ if isinstance(lhapdf_config, list):
2042+ lhapdf_config = lhapdf_config[0]
2043+ lhapdf = misc.import_python_lhapdf(lhapdf_config)
2044+ if not lhapdf:
2045+ return
2046+ lhapdf.setVerbosity(0)
2047+ self.pdfsets = {}
2048+ if isinstance(pdf, str):
2049+ pdf = pdf.split(',')
2050+
2051+ if isinstance(pdf,list) and isinstance(pdf[0],(str,int)):
2052+ self.pdf = []
2053+ for data in pdf:
2054+ if data == 'errorset':
2055+ data = '%s' % self.orig_pdf
2056+ if data == 'central':
2057+ data = '%s@0' % self.orig_pdf
2058+ if '@' in data:
2059+ #one particular dataset
2060+ name, arg = data.rsplit('@',1)
2061+ if int(arg) == 0:
2062+ if name.isdigit():
2063+ self.pdf.append(lhapdf.mkPDF(int(name)))
2064+ else:
2065+ self.pdf.append(lhapdf.mkPDF(name))
2066+ elif name.isdigit():
2067+ try:
2068+ self.pdf.append(lhapdf.mkPDF(int(name)+int(arg)))
2069+ except:
2070+ raise Exception, 'invididual error set need to called with name not with lhapdfID'
2071+ else:
2072+ self.pdf.append(lhapdf.mkPDF(name, int(arg)))
2073+ else:
2074+ if data.isdigit():
2075+ pdfset = lhapdf.mkPDF(int(data)).set()
2076+ else:
2077+ pdfset = lhapdf.mkPDF(data).set()
2078+ self.pdfsets[pdfset.lhapdfID] = pdfset
2079+ self.pdf += pdfset.mkPDFs()
2080+ else:
2081+ self.pdf = pdf
2082+
2083+ for p in self.pdf:
2084+ if p.lhapdfID == self.orig_pdf:
2085+ self.orig_pdf = p
2086+ break
2087+ else:
2088+ self.orig_pdf = lhapdf.mkPDF(self.orig_pdf)
2089+ self.log( "# events generated with PDF: %s (%s)" %(self.orig_pdf.set().name,self.orig_pdf.lhapdfID ))
2090+ # create all the function that need to be called
2091+ self.get_all_fct() # define self.fcts and self.args
2092+
2093+ def run(self, stdout=sys.stdout):
2094+ """ """
2095+ start_time = time.time()
2096+ if self.start_event == 0 or self.force_write_banner:
2097+ lowest_id = self.write_banner(self.output)
2098+ else:
2099+ lowest_id = self.get_id()
2100+
2101+ ids = [lowest_id+i for i in range(len(self.args)-1)]
2102+ all_cross = [0 for i in range(len(self.args))]
2103+
2104+ for nb_event,event in enumerate(self.input):
2105+ if nb_event < self.start_event:
2106+ continue
2107+ elif nb_event >= self.stop_event:
2108+ if self.force_write_banner:
2109+ self.output.write('</LesHouchesEvents>\n')
2110+ break
2111+ if self.is_lo:
2112+ if (nb_event-self.start_event)>=0 and (nb_event-self.start_event) % 2500 ==0:
2113+ self.log( '# currently at event %s [ellapsed time: %.2g s]' % (nb_event, time.time()-start_time))
2114+ else:
2115+ if (nb_event-self.start_event)>=0 and (nb_event-self.start_event) % 1000 ==0:
2116+ self.log( '# currently at event %i [ellapsed time: %.2g s]' % (nb_event, time.time()-start_time))
2117+
2118+ self.new_event() #re-init the caching of alphas/pdf
2119+ if self.is_lo:
2120+ wgts = [self.get_lo_wgt(event, *arg) for arg in self.args]
2121+ else:
2122+ wgts = [self.get_nlo_wgt(event, *arg) for arg in self.args]
2123+
2124+ if wgts[0] == 0:
2125+ print wgts
2126+ print event
2127+ raise Exception
2128+
2129+ wgt = [event.wgt*wgts[i]/wgts[0] for i in range(1,len(wgts))]
2130+ all_cross = [(all_cross[j] + event.wgt*wgts[j]/wgts[0]) for j in range(len(wgts))]
2131+
2132+ rwgt_data = event.parse_reweight()
2133+ rwgt_data.update(zip(ids, wgt))
2134+ event.reweight_order += ids
2135+ # order the
2136+ self.output.write(str(event))
2137+ else:
2138+ self.output.write('</LesHouchesEvents>\n')
2139+ self.output.close()
2140+ self.print_cross_sections(all_cross, min(nb_event,self.stop_event)-self.start_event+1, stdout)
2141+
2142+ if self.output.name != self.output_path:
2143+ import shutil
2144+ shutil.move(self.output.name, self.output_path)
2145+
2146+ return all_cross
2147+
2148+ def print_cross_sections(self, all_cross, nb_event, stdout):
2149+ """print the cross-section."""
2150+
2151+ norm = self.banner.get('run_card', 'event_norm', default='sum')
2152+ #print "normalisation is ", norm
2153+ #print "nb_event is ", nb_event
2154+
2155+ max_scale, min_scale = 0,sys.maxint
2156+ max_alps, min_alps = 0, sys.maxint
2157+ max_dyn, min_dyn = 0,sys.maxint
2158+ pdfs = {}
2159+ dyns = {} # dyn : {'max': , 'min':}
2160+
2161+ if norm == 'sum':
2162+ norm = 1
2163+ elif norm == 'average':
2164+ norm = 1./nb_event
2165+ elif norm == 'unity':
2166+ norm = 1
2167+
2168+ all_cross = [c*norm for c in all_cross]
2169+ stdout.write("# mur\t\tmuf\t\talpsfact\tdynamical_scale\tpdf\t\tcross-section\n")
2170+ for i,arg in enumerate(self.args):
2171+
2172+ to_print = list(arg)
2173+ to_print[4] = to_print[4].lhapdfID
2174+ to_print.append(all_cross[i])
2175+ to_report = []
2176+ stdout.write('%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\n' % tuple(to_print))
2177+
2178+ mur, muf, alps, dyn, pdf = arg[:5]
2179+ if pdf == self.orig_pdf and (dyn==self.orig_dyn or dyn==-1)\
2180+ and (mur!=1 or muf!=1 or alps!=1):
2181+ max_scale = max(max_scale,all_cross[i])
2182+ min_scale = min(min_scale,all_cross[i])
2183+ if pdf == self.orig_pdf and mur==1 and muf==1 and \
2184+ (dyn==self.orig_dyn or dyn==-1) and alps!=1:
2185+ max_alps = max(max_alps,all_cross[i])
2186+ min_alps = min(min_alps,all_cross[i])
2187+
2188+ if pdf == self.orig_pdf and mur==1 and muf==1 and alps==1:
2189+ max_dyn = max(max_dyn,all_cross[i])
2190+ min_dyn = min(min_dyn,all_cross[i])
2191+
2192+ if pdf == self.orig_pdf and (alps!=1 or mur!=1 or muf!=1) and \
2193+ (dyn!=self.orig_dyn or dyn!=-1):
2194+ if dyn not in dyns:
2195+ dyns[dyn] = {'max':0, 'min':sys.maxint,'central':0}
2196+ curr = dyns[dyn]
2197+ curr['max'] = max(curr['max'],all_cross[i])
2198+ curr['min'] = min(curr['min'],all_cross[i])
2199+ if pdf == self.orig_pdf and (alps==1 and mur==1 and muf==1) and \
2200+ (dyn!=self.orig_dyn or dyn!=-1):
2201+ if dyn not in dyns:
2202+ dyns[dyn] = {'max':0, 'min':sys.maxint,'central':all_cross[i]}
2203+ else:
2204+ dyns[dyn]['central'] = all_cross[i]
2205+
2206+ if alps==1 and mur==1 and muf==1 and (dyn==self.orig_dyn or dyn==-1):
2207+ pdfset = pdf.set()
2208+ if pdfset.lhapdfID in self.pdfsets:
2209+ if pdfset.lhapdfID not in pdfs :
2210+ pdfs[pdfset.lhapdfID] = [0] * pdfset.size
2211+ pdfs[pdfset.lhapdfID][pdf.memberID] = all_cross[i]
2212+ else:
2213+ to_report.append('# PDF %s : %s\n' % (pdf.lhapdfID, all_cross[i]))
2214+
2215+ stdout.write('\n')
2216+
2217+ resume = StringIO.StringIO()
2218+
2219+ resume.write( '#***************************************************************************\n')
2220+ resume.write( "#\n")
2221+ resume.write( '# original cross-section: %s\n' % all_cross[0])
2222+ if max_scale:
2223+ 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))
2224+ if max_alps:
2225+ 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))
2226+ if max_dyn and (max_dyn!= all_cross[0] or min_dyn != all_cross[0]):
2227+ 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))
2228+ if self.orig_pdf.lhapdfID in pdfs:
2229+ lhapdfid = self.orig_pdf.lhapdfID
2230+ values = pdfs[lhapdfid]
2231+ pdfset = self.pdfsets[lhapdfid]
2232+ pdferr = pdfset.uncertainty(values)
2233+ resume.write( '# PDF variation: +%2.3g%% -%2.3g%%\n' % (pdferr.errplus*100/all_cross[0], pdferr.errminus*100/all_cross[0]))
2234+ # report error/central not directly linked to the central
2235+ resume.write( "#\n")
2236+ for lhapdfid,values in pdfs.items():
2237+ if lhapdfid == self.orig_pdf.lhapdfID:
2238+ continue
2239+ pdfset = self.pdfsets[lhapdfid]
2240+ pdferr = pdfset.uncertainty(values)
2241+ 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]))
2242+
2243+ 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}' }
2244+ for key, curr in dyns.items():
2245+ if key ==-1:
2246+ continue
2247+ central, maxvalue, minvalue = curr['central'], curr['max'], curr['min']
2248+ if central == 0:
2249+ continue
2250+ if maxvalue == 0:
2251+ resume.write("# dynamical scheme # %s : %g # %s\n" %(key, central, dyn_name[key]))
2252+ else:
2253+ 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]))
2254+
2255+ resume.write('\n'.join(to_report))
2256+
2257+ resume.write( '#***************************************************************************\n')
2258+
2259+ stdout.write(resume.getvalue())
2260+ self.log(resume.getvalue())
2261+
2262+
2263+ def write_banner(self, fsock):
2264+ """create the new banner with the information of the weight"""
2265+
2266+ cid = self.get_id()
2267+ lowest_id = cid
2268+
2269+ in_scale = False
2270+ in_pdf = False
2271+ in_alps = False
2272+
2273+ text = ''
2274+
2275+ for arg in self.args[1:]:
2276+ mur, muf, alps, dyn, pdf = arg[:5]
2277+ if pdf == self.orig_pdf and alps ==1 and (mur!=1 or muf!=1 or dyn!=-1):
2278+ if not in_scale:
2279+ text += "<weightgroup name=\"Central scale variation\" combine=\"envelope\">\n"
2280+ in_scale=True
2281+ elif in_scale:
2282+ if (pdf == self.orig_pdf and alps ==1):
2283+ pass
2284+ else:
2285+ text += "</weightgroup> # scale\n"
2286+ in_scale = False
2287+
2288+ if pdf == self.orig_pdf and mur == muf == 1 and dyn==-1 and alps!=1:
2289+ if not in_alps:
2290+ text += "<weightgroup name=\"Emission scale variation\" combine=\"envelope\">\n"
2291+ in_alps=True
2292+ elif in_alps:
2293+ text += "</weightgroup> # ALPS\n"
2294+ in_alps=False
2295+
2296+ if pdf != self.orig_pdf and mur == muf == 1 and dyn==-1 and alps ==1:
2297+ if pdf.lhapdfID in self.pdfsets:
2298+ if in_pdf:
2299+ text += "</weightgroup> # PDFSET -> PDFSET\n"
2300+ pdfset = self.pdfsets[pdf.lhapdfID]
2301+ text +="<weightgroup name=\"%s\" combine=\"%s\"> # %s: %s\n" %\
2302+ (pdfset.name, pdfset.errorType,pdfset.lhapdfID, pdfset.description)
2303+ in_pdf=pdf.lhapdfID
2304+ elif pdf.memberID == 1 and (pdf.lhapdfID - pdf.memberID) in self.pdfsets:
2305+ if in_pdf:
2306+ text += "</weightgroup> # PDFSET -> PDFSET\n"
2307+ pdfset = self.pdfsets[pdf.lhapdfID - 1]
2308+ text +="<weightgroup name=\"%s\" combine=\"%s\"> # %s: %s\n" %\
2309+ (pdfset.name, pdfset.errorType,pdfset.lhapdfID, pdfset.description)
2310+ in_pdf=pdfset.lhapdfID
2311+ elif in_pdf and pdf.lhapdfID - pdf.memberID != in_pdf:
2312+ text += "</weightgroup> # PDFSET -> PDF\n"
2313+ in_pdf = False
2314+ elif in_pdf:
2315+ text += "</weightgroup> PDF \n"
2316+ in_pdf=False
2317+
2318+
2319+
2320+ tag, info = '',''
2321+ if mur!=1.:
2322+ tag += 'MUR="%s" ' % mur
2323+ info += 'MUR=%s ' % mur
2324+ else:
2325+ tag += 'MUR="%s" ' % mur
2326+ if muf!=1.:
2327+ tag += 'MUF="%s" ' % muf
2328+ info += 'MUF=%s ' % muf
2329+ else:
2330+ tag += 'MUF="%s" ' % muf
2331+
2332+ if alps!=1.:
2333+ tag += 'ALPSFACT="%s" ' % alps
2334+ info += 'alpsfact=%s ' % alps
2335+ if dyn!=-1.:
2336+ tag += 'DYN_SCALE="%s" ' % dyn
2337+ info += 'dyn_scale_choice=%s ' % {1:'sum pt', 2:'HT',3:'HT/2',4:'sqrts'}[dyn]
2338+
2339+ if pdf != self.orig_pdf:
2340+ tag += 'LHAPDF="%s" ' % pdf.lhapdfID
2341+ info += 'LHAPDF=%s MemberID=%s' % (pdf.lhapdfID-pdf.memberID, pdf.memberID)
2342+ else:
2343+ tag += 'LHAPDF="%s" ' % pdf.lhapdfID
2344+
2345+ text +='<weight id="%s" %s> %s </weight>\n' % (cid, tag, info)
2346+ cid+=1
2347+
2348+ if in_scale or in_alps or in_pdf:
2349+ text += "</weightgroup>\n"
2350+
2351+ if 'initrwgt' in self.banner:
2352+ self.banner['initrwgt'] += text
2353+ else:
2354+ self.banner['initrwgt'] = text
2355+
2356+
2357+ self.banner.write(self.output, close_tag=False)
2358+
2359+ return lowest_id
2360+
2361+
2362+ def get_id(self):
2363+
2364+ if 'initrwgt' in self.banner:
2365+ pattern = re.compile('<weight id=(?:\'|\")([_\w]+)(?:\'|\")', re.S+re.I+re.M)
2366+ return max([int(wid) for wid in pattern.findall(self.banner['initrwgt'])])+1
2367+ else:
2368+ return 1
2369+
2370+
2371+
2372+
2373+ def get_all_fct(self):
2374+
2375+ all_args = []
2376+ default = [1.,1.,1.,-1,self.orig_pdf]
2377+ #all_args.append(default)
2378+ pos = {'mur':0, 'muf':1, 'alps':2, 'dyn':3, 'pdf':4}
2379+ done = set()
2380+ for one_block in self.together:
2381+ for name in one_block:
2382+ done.add(name)
2383+ for together in itertools.product(*[getattr(self,name) for name in one_block]):
2384+ new_args = list(default)
2385+ for name,value in zip(one_block, together):
2386+ new_args[pos[name]] = value
2387+ all_args.append(new_args)
2388+ for name in pos:
2389+ if name in done:
2390+ continue
2391+ for value in getattr(self, name):
2392+ new_args = list(default)
2393+ new_args[pos[name]] = value
2394+ all_args.append(new_args)
2395+
2396+ self.args = [default]+ [arg for arg in all_args if arg!= default]
2397+
2398+ self.log( "#Will Compute %s weights per event." % (len(self.args)-1))
2399+ return
2400+
2401+ def new_event(self):
2402+ self.alphas = {}
2403+ self.pdfQ2 = {}
2404+
2405+
2406+ def get_pdfQ(self, pdf, pdg, x, scale):
2407+
2408+ if pdg in [-21,-22]:
2409+ pdg = abs(pdg)
2410+ elif pdg == 0:
2411+ return 1
2412+
2413+ f = pdf.xfxQ(pdg, x, scale)/x
2414+# if f == 0 and pdf.memberID ==0:
2415+# pdfset = pdf.set()
2416+# allnumber= [p.xfxQ(pdg, x, scale) for p in pdfset.mkPDFs()]
2417+# f = pdfset.uncertainty(allnumber).central /x
2418+ return f
2419+
2420+ def get_pdfQ2(self, pdf, pdg, x, scale):
2421+
2422+ if pdg in [-21,-22]:
2423+ pdg = abs(pdg)
2424+ elif pdg == 0:
2425+ return 1
2426+
2427+ if (pdf, pdg,x,scale) in self.pdfQ2:
2428+ return self.pdfQ2[(pdf, pdg,x,scale)]
2429+ f = pdf.xfxQ2(pdg, x, scale)/x
2430+ self.pdfQ2[(pdf, pdg,x,scale)] = f
2431+ return f
2432+
2433+ #one method to handle the nnpd2.3 problem -> now just move to central
2434+ if f == 0 and pdf.memberID ==0:
2435+ # to avoid problem with nnpdf2.3 in lhapdf6.1.6
2436+ #print 'central pdf returns 0', pdg, x, scale
2437+ #print self.pdfsets
2438+ pdfset = pdf.set()
2439+ allnumber= [0] + [self.get_pdfQ2(p, pdg, x, scale) for p in pdfset.mkPDFs()[1:]]
2440+ f = pdfset.uncertainty(allnumber).central
2441+ self.pdfQ2[(pdf, pdg,x,scale)] = f
2442+ return f
2443+
2444+ def get_lo_wgt(self,event, Dmur, Dmuf, Dalps, dyn, pdf):
2445+ """
2446+ pdf is a lhapdf object!"""
2447+
2448+ loinfo = event.parse_lo_weight()
2449+
2450+ if dyn == -1:
2451+ mur = loinfo['ren_scale']
2452+ muf1 = loinfo['pdf_q1'][-1]
2453+ muf2 = loinfo['pdf_q2'][-1]
2454+ else:
2455+ if dyn == 1:
2456+ mur = event.get_et_scale(1.)
2457+ elif dyn == 2:
2458+ mur = event.get_ht_scale(1.)
2459+ elif dyn == 3:
2460+ mur = event.get_ht_scale(0.5)
2461+ elif dyn == 4:
2462+ mur = event.get_sqrts_scale(1.)
2463+ muf1 = mur
2464+ muf2 = mur
2465+ loinfo['pdf_q1'][-1] = mur
2466+ loinfo['pdf_q2'][-1] = mur
2467+
2468+
2469+ # MUR part
2470+ wgt = pdf.alphasQ(Dmur*mur)**loinfo['n_qcd']
2471+ # MUF/PDF part
2472+ wgt *= self.get_pdfQ(pdf, self.b1*loinfo['pdf_pdg_code1'][-1], loinfo['pdf_x1'][-1], Dmuf*muf1)
2473+ wgt *= self.get_pdfQ(pdf, self.b2*loinfo['pdf_pdg_code2'][-1], loinfo['pdf_x2'][-1], Dmuf*muf2)
2474+
2475+ for scale in loinfo['asrwt']:
2476+ wgt *= pdf.alphasQ(Dalps*scale)
2477+
2478+ # ALS part
2479+ for i in range(loinfo['n_pdfrw1']-1):
2480+ scale = min(Dalps*loinfo['pdf_q1'][i], Dmuf*muf1)
2481+ wgt *= self.get_pdfQ(pdf, self.b1*loinfo['pdf_pdg_code1'][i], loinfo['pdf_x1'][i], scale)
2482+ wgt /= self.get_pdfQ(pdf, self.b1*loinfo['pdf_pdg_code1'][i], loinfo['pdf_x1'][i+1], scale)
2483+
2484+ for i in range(loinfo['n_pdfrw2']-1):
2485+ scale = min(Dalps*loinfo['pdf_q2'][i], Dmuf*muf2)
2486+ wgt *= self.get_pdfQ(pdf, self.b2*loinfo['pdf_pdg_code2'][i], loinfo['pdf_x2'][i], scale)
2487+ wgt /= self.get_pdfQ(pdf, self.b2*loinfo['pdf_pdg_code2'][i], loinfo['pdf_x2'][i+1], scale)
2488+
2489+ return wgt
2490+
2491+ def get_nlo_wgt(self,event, Dmur, Dmuf, Dalps, dyn, pdf):
2492+ """return the new weight for NLO event --with weight information-- """
2493+
2494+ wgt = 0
2495+ nloinfo = event.parse_nlo_weight()
2496+ for cevent in nloinfo.cevents:
2497+ if dyn == 1:
2498+ mur2 = cevent.get_et_scale(1.)**2
2499+ elif dyn == 2:
2500+ mur2 = cevent.get_ht_scale(1.)**2
2501+ elif dyn == 3:
2502+ mur2 = cevent.get_ht_scale(0.5)**2
2503+ elif dyn == 4:
2504+ mur2 = cevent.get_sqrts_scale(event,1)**2
2505+ else:
2506+ mur2 = 0
2507+ muf2 = mur2
2508+
2509+ for onewgt in cevent.wgts:
2510+ if dyn == -1:
2511+ mur2 = onewgt.scales2[1]
2512+ muf2 = onewgt.scales2[2]
2513+ Q2 = onewgt.scales2[0] # Ellis-Sexton scale
2514+
2515+ wgtpdf = self.get_pdfQ2(pdf, self.b1*onewgt.pdgs[0], onewgt.bjks[0],
2516+ Dmuf**2 * muf2)
2517+ wgtpdf *= self.get_pdfQ2(pdf, self.b2*onewgt.pdgs[1], onewgt.bjks[1],
2518+ Dmuf**2 * muf2)
2519+
2520+ tmp = onewgt.pwgt[0]
2521+ tmp += onewgt.pwgt[1] * math.log(Dmur**2 * mur2/ Q2)
2522+ tmp += onewgt.pwgt[2] * math.log(Dmuf**2 * muf2/ Q2)
2523+ tmp *= math.sqrt(4*math.pi*pdf.alphasQ2(Dmur**2*mur2))**onewgt.qcdpower
2524+
2525+ if wgtpdf == 0: #happens for nn23pdf due to wrong set in lhapdf
2526+ key = (self.b1*onewgt.pdgs[0], self.b2*onewgt.pdgs[1], onewgt.bjks[0],onewgt.bjks[1], muf2)
2527+ if dyn== -1 and Dmuf==1 and Dmur==1 and pdf==self.orig_pdf:
2528+ wgtpdf = onewgt.ref_wgt / tmp
2529+ self.pdfQ2[key] = wgtpdf
2530+ elif key in self.pdfQ2:
2531+ wgtpdf = self.pdfQ2[key]
2532+ else:
2533+ # real zero!
2534+ wgtpdf = 0
2535+
2536+ tmp *= wgtpdf
2537+ wgt += tmp
2538+
2539+ if __debug__ and dyn== -1 and Dmur==1 and Dmuf==1 and pdf==self.orig_pdf:
2540+ if not misc.equal(tmp, onewgt.ref_wgt, sig_fig=4):
2541+ misc.sprint(tmp, onewgt.ref_wgt, (tmp-onewgt.ref_wgt)/tmp)
2542+ misc.sprint(onewgt)
2543+ misc.sprint(cevent)
2544+ raise Exception, 'not enough agreement between stored value and computed one'
2545+
2546+
2547+ return wgt
2548+
2549+
2550+def call_systematics(args, result=sys.stdout, running=True,
2551+ log=lambda x:sys.stdout.write(str(x)+'\n')):
2552+ """calling systematics from a list of arguments"""
2553+
2554+ input, output = args[0:2]
2555+ opts = {}
2556+ for arg in args[2:]:
2557+ if '=' in arg:
2558+ key,values= arg.split('=')
2559+ key = key.replace('-','')
2560+ values = values.split(',')
2561+ if key == 'together':
2562+ if key in opts:
2563+ opts[key].append(tuple(values))
2564+ else:
2565+ opts[key]=[tuple(values)]
2566+ elif key == 'result':
2567+ result = open(values[0],'w')
2568+ elif key in ['start_event', 'stop_event']:
2569+ opts[key] = int(values[0])
2570+ elif key == 'write_banner':
2571+ opts[key] = banner_mod.ConfigFile.format_variable(values[0], bool, 'write_banner')
2572+ else:
2573+ if key in opts:
2574+ opts[key] += values
2575+ else:
2576+ opts[key] = values
2577+ else:
2578+ raise SystematicsError, "unknow argument", arg
2579+
2580+ #load run_card and extract parameter if needed.
2581+ if 'from_card' in opts:
2582+ if opts['from_card'] != ['internal']:
2583+ card = banner.RunCard(opts['from_card'][0])
2584+ else:
2585+ lhe = lhe_parser.EventFile(input)
2586+ card = banner.RunCard(banner.Banner(lhe.banner)['mgruncard'])
2587+ lhe.close()
2588+
2589+ if isinstance(card, banner.RunCardLO):
2590+ # LO case
2591+ opts['mur'] = [float(x) for x in card['sys_scalefact'].split()]
2592+ opts['muf'] = opts['mur']
2593+ if card['sys_alpsfact'] != 'None':
2594+ opts['alps'] = [float(x) for x in card['sys_alpsfact'].split()]
2595+ else:
2596+ opts['alps'] = [1.0]
2597+ opts['together'] = [('mur','muf','alps','dyn')]
2598+ pdfs = card['sys_pdf'].split('&&')
2599+ opts['dyn'] = [-1,1,2,3,4]
2600+ opts['pdf'] = []
2601+ for pdf in pdfs:
2602+ split = pdf.split()
2603+ if len(split)==1:
2604+ opts['pdf'].append('%s' %pdf)
2605+ else:
2606+ pdf,nb = split
2607+ for i in range(int(nb)):
2608+ opts['pdf'].append('%s@%s' % (pdf, i))
2609+ if not opts['pdf']:
2610+ opts['pdf'] = 'central'
2611+ else:
2612+ #NLO case
2613+ raise Exception
2614+ del opts['from_card']
2615+
2616+
2617+ obj = Systematics(input, output, log=log,**opts)
2618+ if running:
2619+ obj.run(result)
2620+ return obj
2621+
2622+if __name__ == "__main__":
2623+ call_systematics(sys.argv[1:])
2624+
2625+
2626+
2627+
2628+
2629+
2630+
2631+
2632
2633=== modified file 'tests/acceptance_tests/test_cmd_madevent.py'
2634--- tests/acceptance_tests/test_cmd_madevent.py 2016-08-17 16:25:19 +0000
2635+++ tests/acceptance_tests/test_cmd_madevent.py 2016-08-24 22:15:21 +0000
2636@@ -117,9 +117,10 @@
2637 interface.exec_cmd('set madanalysis_path %s --no_save' % pjoin(MG5DIR, 'MadAnalysis'))
2638 interface.onecmd('output madevent %s -f' % self.run_dir)
2639
2640- if not os.path.exists(pjoin(interface.options['syscalc_path'],'sys_calc')):
2641- print "install SysCalc"
2642- interface.onecmd('install SysCalc')
2643+ if os.path.exists(pjoin(interface.options['syscalc_path'],'sys_calc')):
2644+ shutil.rmtree(interface.options['syscalc_path'])
2645+ #print "install SysCalc"
2646+ #interface.onecmd('install SysCalc')
2647
2648
2649 self.cmd_line = MECmd.MadEventCmdShell(me_dir=self.run_dir)
2650@@ -475,10 +476,10 @@
2651 # check that the html has the information
2652 self.assertTrue('syst' in data[0].parton)
2653 # check that the code was runned correctly
2654- fsock = open('%s/Events/%s/%s_parton_syscalc.log' % \
2655- (self.run_dir, data[0]['run_name'], data[0]['tag']),'r')
2656+ fsock = open('%s/Events/%s/parton_systematics.log' % \
2657+ (self.run_dir, data[0]['run_name']),'r')
2658 text = fsock.read()
2659- self.assertTrue(text.count('cross-section') >= 3)
2660+ self.assertTrue(text.count('dynamical scheme') >= 3)
2661
2662
2663 def check_pythia_output(self, run_name='run_01', syst=False):
2664@@ -488,9 +489,9 @@
2665 self.assertTrue('lhe' in data[0].pythia)
2666 self.assertTrue('log' in data[0].pythia)
2667
2668- if syst:
2669- # check that the html has the information
2670- self.assertTrue('rwt' in data[0].pythia)
2671+# if syst:
2672+# # check that the html has the information
2673+# self.assertTrue('rwt' in data[0].pythia)
2674
2675 def check_matched_plot(self, run_name='run_01', mintime=None, tag='fermi'):
2676 """ """
2677
2678=== modified file 'tests/input_files/run_card_matching.dat'
2679--- tests/input_files/run_card_matching.dat 2013-11-22 13:16:07 +0000
2680+++ tests/input_files/run_card_matching.dat 2016-08-24 22:15:21 +0000
2681@@ -234,9 +234,9 @@
2682 # will be use by SysCalc (if install)
2683 #**************************************
2684 #
2685-0.5 = sys_scalefact # Central scale factors
2686-0.5 = sys_alpsfact # \alpha_s emission scale factors
2687+0.5 1 2= sys_scalefact # Central scale factors
2688+1= sys_alpsfact # \alpha_s emission scale factors
2689 30 = sys_matchscale # variation of merging scale
2690 # PDF sets and number of members (0 or none for all members).
2691-CT10nlo.LHgrid 1 = sys_pdf # matching scales
2692+CT10nlo 1 = sys_pdf # matching scales
2693 # MSTW2008nlo68cl.LHgrid 1 = sys_pdf

Subscribers

People subscribed via source and target branches

to all changes: