Merge lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4

Proposed by marco zaro
Status: Merged
Merged at revision: 344
Proposed branch: lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory
Merge into: lp:~maddevelopers/mg5amcnlo/2.3.4
Diff against target: 20865 lines (+19155/-274)
106 files modified
UpdateNotes.txt (+6/-0)
madgraph/fks/fks_base.py (+79/-60)
madgraph/fks/fks_helas_objects.py (+449/-45)
madgraph/interface/amcatnlo_interface.py (+171/-49)
madgraph/interface/madgraph_interface.py (+6/-3)
madgraph/iolibs/export_fks.py (+38/-42)
madgraph/iolibs/template_files/loop_optimized/compute_color_flows.inc (+11/-11)
madgraph/loop/loop_exporters.py (+8/-3)
madgraph/various/banner.py (+1/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f (+278/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f (+760/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f (+423/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f (+540/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa.f (+756/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f (+524/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%coef_specs.inc (+6/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%improve_ps.f (+981/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f (+2267/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_num.f (+131/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%mp_compute_loop_coefs.f (+296/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%nexternal.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%ngraphs.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%nsqso_born.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%nsquaredSO.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%pmass.inc (+3/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%polynomial.f (+455/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%b_sf_001.f (+153/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born.f (+292/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born_conf.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born_hel.f (+151/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born_leshouche.inc (+8/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born_maxamps.inc (+3/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born_ngraphs.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%born_nhel.inc (+3/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%coloramps.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%config_subproc_map.inc (+1/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%configs_and_props_decl.inc (+12/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%fks_info.inc (+53/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%get_color.f (+48/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%leshouche_decl.inc (+7/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_1.f (+172/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_2.f (+172/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%matrix_3.f (+172/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%nFKSconfigs.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%ncombs.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%nexternal.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%ngraphs.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%parton_lum_1.f (+112/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%parton_lum_2.f (+110/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%parton_lum_3.f (+110/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%parton_lum_chooser.f (+23/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%pmass.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%real_from_born_configs.inc (+6/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%real_me_chooser.f (+22/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%sborn_sf.f (+17/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%sborn_sf_dum.f (+13/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%CT_interface.f (+278/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%GOLEM_interface.f (+760/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%TIR_interface.f (+423/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%born_matrix.f (+540/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%check_sa.f (+756/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%check_sa_born_splitOrders.f (+524/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%coef_specs.inc (+6/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%improve_ps.f (+981/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f (+2267/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_num.f (+131/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%mp_compute_loop_coefs.f (+296/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%nexternal.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%ngraphs.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%nsqso_born.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%nsquaredSO.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%pmass.inc (+3/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%polynomial.f (+455/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%b_sf_001.f (+153/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born.f (+292/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born_conf.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born_hel.f (+151/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born_leshouche.inc (+8/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born_maxamps.inc (+3/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born_ngraphs.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%born_nhel.inc (+3/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%coloramps.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%config_subproc_map.inc (+1/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%configs_and_props_decl.inc (+12/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%fks_info.inc (+53/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%get_color.f (+48/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%leshouche_decl.inc (+7/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_1.f (+172/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_2.f (+172/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%matrix_3.f (+172/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%nFKSconfigs.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%ncombs.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%nexternal.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%ngraphs.inc (+2/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%parton_lum_1.f (+112/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%parton_lum_2.f (+110/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%parton_lum_3.f (+110/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%parton_lum_chooser.f (+23/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%pmass.inc (+4/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%real_from_born_configs.inc (+6/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%real_me_chooser.f (+22/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%sborn_sf.f (+17/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%sborn_sf_dum.f (+13/-0)
tests/time_db (+59/-57)
tests/unit_tests/iolibs/test_export_fks.py (+93/-0)
tests/unit_tests/iolibs/test_export_fks_EW.py (+6/-3)
To merge this branch: bzr merge lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory
Reviewer Review Type Date Requested Status
Valentin Hirschi Approve
Olivier Mattelaer Approve
Rikkert Frederix Approve
Stefano Frixione Pending
MadDevelopers Pending
Review via email: mp+286337@code.launchpad.net

Description of the change

A new generation mode at NLO has been implemented, in order to spread the generation on all available CPUs and to use less memory (writing partial results to disk). The new generation mode can be triggered with
set new_nlo_generation True from the interface.
It can really be amazingly faster and more memory-efficient than before (Rikkert generated ttjj@NLO in 10 -ten- minutes on a 20-core cluster node), so that we can try generating higher multiplicities and look what the new bottlenecks are.
 The output should be identical (I have added an unit test to check this, yet for a very simple process : pp>w+)

As many threads start together, the logging is not as nice as the old one.

I have to thank Josh Bendavid for having provided the code for this.

Thanks for having a look and send comments...
Cheers,

Marco

To post a comment you must log in.
Revision history for this message
marco zaro (marco-zaro) wrote :

just to be clear, we are speaking about code generation (not event generation, unfortunately...)

Revision history for this message
Rikkert Frederix (frederix) wrote :

Hi Marco,

I would rename the option "new_nlo_generation" to something like "low_memory_multicore_nlo_generation" (or something shorter/better).

It seems like Ctrl-C doesn't work correctly during the 'output' of the process.

Cheers,
Rik

review: Needs Fixing
Revision history for this message
marco zaro (marco-zaro) wrote :

Hi Rik
On 18 Feb 2016, at 09:29, Rikkert Frederix <email address hidden> wrote:

> Review: Needs Fixing
>
> Hi Marco,
>
> I would rename the option "new_nlo_generation" to something like "low_memory_multicore_nlo_generation" (or something shorter/better).
>
I did not think much about the name

> It seems like Ctrl-C doesn't work correctly during the 'output' of the process.
>
argh, that was something we have already seen when we used the same multicore threading for compiling…
I will have a look and let you know…
Cheers,

Marco

> Cheers,
> Rik
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

311. By marco zaro

ctrl+c should be caught now

Revision history for this message
marco zaro (marco-zaro) wrote :

Ciao Rik,
catching ctrl+c should be ok now.
Another thing i have noticed, is that by default all available cores are used… The nb_core option should be used instead for that, i will fix also this issue
Ciao,

Marco
On 18 Feb 2016, at 10:33, marco zaro <email address hidden> wrote:

> Hi Rik
> On 18 Feb 2016, at 09:29, Rikkert Frederix <email address hidden> wrote:
>
>> Review: Needs Fixing
>>
>> Hi Marco,
>>
>> I would rename the option "new_nlo_generation" to something like "low_memory_multicore_nlo_generation" (or something shorter/better).
>>
> I did not think much about the name
>
>> It seems like Ctrl-C doesn't work correctly during the 'output' of the process.
>>
> argh, that was something we have already seen when we used the same multicore threading for compiling…
> I will have a look and let you know…
> Cheers,
>
> Marco
>
>> Cheers,
>> Rik
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
>> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

312. By marco zaro

changed the name of the option (new_nlo_generation -> low_mem_multicore_nlo_generation).
Now if nb_core is set, the number of used cores is determined by that

Revision history for this message
marco zaro (marco-zaro) wrote :

Hi,
I have changed the name of the option (low_mem_multicore_nlo_generation). Also the number of cores is dictated by nb_cores (if set)
cheers,

Marco
On 18 Feb 2016, at 12:24, marco zaro <email address hidden> wrote:

> Ciao Rik,
> catching ctrl+c should be ok now.
> Another thing i have noticed, is that by default all available cores are used… The nb_core option should be used instead for that, i will fix also this issue
> Ciao,
>
> Marco
> On 18 Feb 2016, at 10:33, marco zaro <email address hidden> wrote:
>
>> Hi Rik
>> On 18 Feb 2016, at 09:29, Rikkert Frederix <email address hidden> wrote:
>>
>>> Review: Needs Fixing
>>>
>>> Hi Marco,
>>>
>>> I would rename the option "new_nlo_generation" to something like "low_memory_multicore_nlo_generation" (or something shorter/better).
>>>
>> I did not think much about the name
>>
>>> It seems like Ctrl-C doesn't work correctly during the 'output' of the process.
>>>
>> argh, that was something we have already seen when we used the same multicore threading for compiling…
>> I will have a look and let you know…
>> Cheers,
>>
>> Marco
>>
>>> Cheers,
>>> Rik
>>>
>>> --
>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
>>> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.
>>
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
>> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

Revision history for this message
Rikkert Frederix (frederix) :
review: Approve
Revision history for this message
Olivier Mattelaer (olivier-mattelaer) wrote :

Hi Marco,

That's so cool, I missed the fact that a library was allowing to evade the GIL limitation.
This is so much nicer than using the Threading library (like the multi-core option of MG).
I guess that we should think to change the multi-core options of MG to move to this multiprocessing library.

1) Is there any reason why this mode is not set as default?
2) It should be indicated in the UpdateNotes how to use it.
3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
4)can you check why you set self['has_loops'] = False (line 466 of the diff) around line 220 of madgraph/fks/fks_helas_objects.py
6) Would not be better to not write the data on a file (via pickle), but add them in a multiprocessing.queue object that can be used by the main process? Or do you see a problem with that?
7) in madgraph/interface/amcatnlo_interface.py you define glob_directories_map
Defining a non static global is always annoying and creates troubles very often ( I have to apologize to have done it in aloha). Is it really necessary? and is it thread safe?

8) Is it possible to have a printout when the last process is generated and that you collect the info (I guess that you do that) since it hangs on one core on a line of type:
INFO: Generating born process: s~ s > t t t~ t~ [ all = QCD ]
I guess he is doing something else actually.

So this is very very nice. Nice work!!

Cheers,

Olivier

review: Needs Information
313. By Valentin Hirschi

1. Fixed function 'fix_coef_specs'

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

Hi Marco,

Awesome job, really.

On the testing side, I noticed the following:

a) Doing a full comparison of p p > t t~ j [QCD] generated in single and multicore, I noticed that the file 'SubProcesses/proc_characteristics' was incorrect in mutlicore, since it contains:
 ninitial = 0
 nexternal = 0

b) This implementation broke process generation with another OLP than MadLoop (i.e. Gosam).
When selecting Gosam as OLP, the generation crashes at:

FKSProcessError : Wrong type of argument for matrix_elements in function write_lh_order.

from

[...]
  File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/iolibs/export_fks.py", line 1355, in generate_virtuals_from_OLP
    self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)
  File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/iolibs/export_fks.py", line 1623, in write_lh_order
    'matrix_elements in function write_lh_order.')

Of course we don't want to have any sort of parallelization on behalf of GoSam, but normally here the only thing that is done is the writing of the leshouches order. The actual generation with GoSam should take place at the very end, basically after that the process is generated and the parallelization is done. So if this is not too difficult to fix, I think it would be worth doing so. Otherwise, explicitly forbid this mode when OLP is not MadLoop, but that would be weird because there is no reason for doing so.
Note that you don't need to have GoSam installed to test this I think, because it crashes when writing the leshouches order file, so before actually calling GoSam.

c) This is a problem which was not introduced by this branch but that I realized now and needs to be fixed. It showed up because I saw that you had to play around with the overall 'max_loop_coefs' value for different FKSHelasMatrixElements.
So imagine you do the following:

MG5_aMC>import model HC_NLO_X0_UFO-heft
MG5_aMC>generate g g > x0 [virt=QCD]
MG5_aMC>add process g g > z [virt=QCD]
MG5_aMC>output MyTest

Then you will see the message:
"""
ML5 has just output processes which do not share the same maximum loop wavefunction size or the same maximum loop vertex rank. This is potentially dangerous. Please prefer to output them separately.
"""

This is because if you look at the different coef_specs.inc, you will find:

> cat MyTest/SubProcesses/P0_gg_x0/coef_specs.inc
      INTEGER MAXLWFSIZE
      PARAMETER (MAXLWFSIZE=4)
      INTEGER LOOP_MAXCOEFS
      PARAMETER (LOOP_MAXCOEFS=70)
      INTEGER VERTEXMAXCOEFS
      PARAMETER (VERTEXMAXCOEFS=15)

> cat MyTest/SubProcesses/P1_gg_z/coef_specs.inc
      INTEGER MAXLWFSIZE
      PARAMETER (MAXLWFSIZE=4)
      INTEGER LOOP_MAXCOEFS
      PARAMETER (LOOP_MAXCOEFS=35)
      INTEGER VERTEXMAXCOEFS
      PARAMETER (VERTEXMAXCOEFS=5)

But then you see that in Source/DHELAS you have:

> cat MyTest/Source/DHELAS/coef_specs.inc
      INTEGER MAXLWFSIZE
      PARAMETER (MAXLWFSIZE=4)
      INTEGER VERTEXMAXCOEFS
      PARAMETER (VERTEXMAXCOEFS=15)

The tricky thing here is that DHELAS, which is used by both subprocesses, depends on VERTEXMAXCOEFS which is different for both subprocesses.
It is actually fine but *...

Read more...

review: Needs Fixing
Revision history for this message
marco zaro (marco-zaro) wrote :

Ciao Olivier,
thanks for the review
On 18 Feb 2016, at 20:44, Olivier Mattelaer <email address hidden> wrote:

> Review: Needs Information
>
> Hi Marco,
>
> That's so cool, I missed the fact that a library was allowing to evade the GIL limitation.
> This is so much nicer than using the Threading library (like the multi-core option of MG).
> I guess that we should think to change the multi-core options of MG to move to this multiprocessing library.
>
>
> 1) Is there any reason why this mode is not set as default?
for simple processes, the old mode works OK, and has a more detailed logging (number of diagrams, etc).
> 2) It should be indicated in the UpdateNotes how to use it.
done
> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
i don’t know… many of the interface options do not make sense for the amplitude generation, and will therefore have to be removed…

> 4)can you check why you set self['has_loops'] = False (line 466 of the diff) around line 220 of madgraph/fks/fks_helas_objects.py
it starts with false, but then, after the calls to
async_generate_born
it is updated
                has_loops = bornout[2]
                self['has_loops'] = self['has_loops'] or has_loops

> 6) Would not be better to not write the data on a file (via pickle), but add them in a multiprocessing.queue object that can be used by the main process? Or do you see a problem with that?
 The main reason to use pkl files is to keep the RAM as free as possible
> 7) in madgraph/interface/amcatnlo_interface.py you define glob_directories_map
> Defining a non static global is always annoying and creates troubles very often ( I have to apologize to have done it in aloha). Is it really necessary? and is it thread safe?
It can be put as global. Josh said that this is related to memory/cpu performance. I think it is thread safe
> 8) Is it possible to have a printout when the last process is generated and that you collect the info (I guess that you do that) since it hangs on one core on a line of type:
> INFO: Generating born process: s~ s > t t t~ t~ [ all = QCD ]
> I guess he is doing something else actually.
i have added some logging
>
> So this is very very nice. Nice work!!
>
I have pushed the changes, now I’ll look at valentin’s comments…

Cheers,

Marco

> Cheers,
>
> Olivier
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

314. By marco zaro

update with olivier's comments

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

Thanks,

>> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
> i don’t know… many of the interface options do not make sense for the amplitude generation, and will therefore have to be removed…

No just keep self.options and do not touch to it
such that this is just a pointer to the one of the interface
This will not use any memory (just a pointer in term of memory), actually less that what you do right now. (since you store two additional value if we count OLP)
So yes you will have access to pointless information but I do not see the problem.

Thanks a lot for the rest.

Cheers,

Olivier

> On Feb 19, 2016, at 09:33, marco zaro <email address hidden> wrote:
>
> Ciao Olivier,
> thanks for the review
> On 18 Feb 2016, at 20:44, Olivier Mattelaer <email address hidden> wrote:
>
>> Review: Needs Information
>>
>> Hi Marco,
>>
>> That's so cool, I missed the fact that a library was allowing to evade the GIL limitation.
>> This is so much nicer than using the Threading library (like the multi-core option of MG).
>> I guess that we should think to change the multi-core options of MG to move to this multiprocessing library.
>>
>>
>> 1) Is there any reason why this mode is not set as default?
> for simple processes, the old mode works OK, and has a more detailed logging (number of diagrams, etc).
>> 2) It should be indicated in the UpdateNotes how to use it.
> done
>> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
> i don’t know… many of the interface options do not make sense for the amplitude generation, and will therefore have to be removed…
>
>> 4)can you check why you set self['has_loops'] = False (line 466 of the diff) around line 220 of madgraph/fks/fks_helas_objects.py
> it starts with false, but then, after the calls to
> async_generate_born
> it is updated
> has_loops = bornout[2]
> self['has_loops'] = self['has_loops'] or has_loops
>
>> 6) Would not be better to not write the data on a file (via pickle), but add them in a multiprocessing.queue object that can be used by the main process? Or do you see a problem with that?
> The main reason to use pkl files is to keep the RAM as free as possible
>> 7) in madgraph/interface/amcatnlo_interface.py you define glob_directories_map
>> Defining a non static global is always annoying and creates troubles very often ( I have to apologize to have done it in aloha). Is it really necessary? and is it thread safe?
> It can be put as global. Josh said that this is related to memory/cpu performance. I think it is thread safe
>> 8) Is it possible to have a printout when the last process is generated and that you collect the info (I guess that you do that) sinc...

Read more...

315. By marco zaro

nexternal and ninitial are correctly written in proc_characteristics

316. By marco zaro

fixed bug occurring with write_lh_order for OLP!=MadLoop

Revision history for this message
marco zaro (marco-zaro) wrote :
Download full text (7.2 KiB)

Ciao Valentin,
thanks for your review, so:

On 18 Feb 2016, at 23:48, Valentin Hirschi <email address hidden> wrote:

> Review: Needs Fixing
>
> Hi Marco,
>
> Awesome job, really.
>
> On the testing side, I noticed the following:
>
> a) Doing a full comparison of p p > t t~ j [QCD] generated in single and multicore, I noticed that the file 'SubProcesses/proc_characteristics' was incorrect in mutlicore, since it contains:
> ninitial = 0
> nexternal = 0
hmm right, i was not checking this file with the unit test. I have added it to the test, it should be fixed now (rev 315)
>
> b) This implementation broke process generation with another OLP than MadLoop (i.e. Gosam).
> When selecting Gosam as OLP, the generation crashes at:
>
> FKSProcessError : Wrong type of argument for matrix_elements in function write_lh_order.
>
> from
>
> [...]
> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/iolibs/export_fks.py", line 1355, in generate_virtuals_from_OLP
> self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)
> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/iolibs/export_fks.py", line 1623, in write_lh_order
> 'matrix_elements in function write_lh_order.’)
>
> Of course we don't want to have any sort of parallelization on behalf of GoSam, but normally here the only thing that is done is the writing of the leshouches order. The actual generation with GoSam should take place at the very end, basically after that the process is generated and the parallelization is done. So if this is not too difficult to fix, I think it would be worth doing so. Otherwise, explicitly forbid this mode when OLP is not MadLoop, but that would be weird because there is no reason for doing so.
> Note that you don't need to have GoSam installed to test this I think, because it crashes when writing the leshouches order file, so before actually calling GoSam.
>
I have added a unit test also for this and reproduced your error. It should be fixed now (rev 316), please check that i did not screw anything, as some modifications to the code were needed (the unit tests pass after trival fixes)

> c) This is a problem which was not introduced by this branch but that I realized now and needs to be fixed. It showed up because I saw that you had to play around with the overall 'max_loop_coefs' value for different FKSHelasMatrixElements.
> So imagine you do the following:
>
> MG5_aMC>import model HC_NLO_X0_UFO-heft
> MG5_aMC>generate g g > x0 [virt=QCD]
> MG5_aMC>add process g g > z [virt=QCD]
> MG5_aMC>output MyTest
>
> Then you will see the message:
> """
> ML5 has just output processes which do not share the same maximum loop wavefunction size or the same maximum loop vertex rank. This is potentially dangerous. Please prefer to output them separately.
> """
>
> This is because if you look at the different coef_specs.inc, you will find:
>
>> cat MyTest/SubProcesses/P0_gg_x0/coef_specs.inc
> INTEGER MAXLWFSIZE
> PARAMETER (MAXLWFSIZE=4)
> INTEGER LOOP_MAXCOEFS
> PARAMETER (LOOP_MAXCOEFS=70)
> INTEGER VERTEXMAXCOEFS
> PARAMETER (VERTEXMAXCOEFS=15)
>
>> cat MyTest/SubProcess...

Read more...

317. By Valentin Hirschi

1. Fixed an issue of the FKS output when the MG option loop_color_flow is set to True.

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

Hi Marco,

One more thing, please check how you defined the possible values for the option 'low_mem_multicore_nlo_generation' that you have introduced because the auto completion to 'True False default' does not work.

> Ciao Valentin,
> thanks for your review, so:

> On 18 Feb 2016, at 23:48, Valentin Hirschi <email address hidden> wrote:
>
> > Review: Needs Fixing
> >
> > Hi Marco,
> >
> > Awesome job, really.
> >
> > On the testing side, I noticed the following:
> >
> > a) Doing a full comparison of p p > t t~ j [QCD] generated in single and
> multicore, I noticed that the file 'SubProcesses/proc_characteristics' was
> incorrect in mutlicore, since it contains:
> > ninitial = 0
> > nexternal = 0
> hmm right, i was not checking this file with the unit test. I have added it to
> the test, it should be fixed now (rev 315)
> >
> > b) This implementation broke process generation with another OLP than
> MadLoop (i.e. Gosam).
> > When selecting Gosam as OLP, the generation crashes at:
> >
> > FKSProcessError : Wrong type of argument for matrix_elements in function
> write_lh_order.
> >
> > from
> >
> > [...]
> > File "/Users/valentin/Documents/Work/MG5/2.3.4-low-
> memory/madgraph/iolibs/export_fks.py", line 1355, in
> generate_virtuals_from_OLP
> > self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)
> > File "/Users/valentin/Documents/Work/MG5/2.3.4-low-
> memory/madgraph/iolibs/export_fks.py", line 1623, in write_lh_order
> > 'matrix_elements in function write_lh_order.’)
> >
> > Of course we don't want to have any sort of parallelization on behalf of
> GoSam, but normally here the only thing that is done is the writing of the
> leshouches order. The actual generation with GoSam should take place at the
> very end, basically after that the process is generated and the
> parallelization is done. So if this is not too difficult to fix, I think it
> would be worth doing so. Otherwise, explicitly forbid this mode when OLP is
> not MadLoop, but that would be weird because there is no reason for doing so.
> > Note that you don't need to have GoSam installed to test this I think,
> because it crashes when writing the leshouches order file, so before actually
> calling GoSam.
> >
> I have added a unit test also for this and reproduced your error. It should be
> fixed now (rev 316), please check that i did not screw anything, as some
> modifications to the code were needed (the unit tests pass after trival fixes)

So when not using the low-mem option, now OLP = GoSam crashes with:

AttributeError : 'MasterCmd' object has no attribute 'born_processes_for_olp'

at

[...]
  File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/interface/amcatnlo_interface.py", line 557, in do_output
    self.born_processes_for_olp,self._export_dir,self.options['OLP'])
AttributeError: 'MasterCmd' object has no attribute 'born_processes_for_olp'

And when doing it with the new option turned on, it goes through but crashes at the compilation because apparently you still linking the version of 'BinothLHA.f' for MadLoop and not the one for GoSam as it is done when not using the new option.

>
> > c) This is a problem which was n...

Read more...

318. By marco zaro

fixed autocompletion, behaviour with gosam and [real=QCD]

Revision history for this message
marco zaro (marco-zaro) wrote :
Download full text (10.0 KiB)

Ciao Valentin,
Thanks for the review…
I should have addressed all points now, see below…
cheers,

Marco
On 19 Feb 2016, at 21:04, Valentin Hirschi <email address hidden> wrote:

> Hi Marco,
>
> One more thing, please check how you defined the possible values for the option 'low_mem_multicore_nlo_generation' that you have introduced because the auto completion to 'True False default' does not work.

fixed
>
>> Ciao Valentin,
>> thanks for your review, so:
>
>> On 18 Feb 2016, at 23:48, Valentin Hirschi <email address hidden> wrote:
>>
>>> Review: Needs Fixing
>>>
>>> Hi Marco,
>>>
>>> Awesome job, really.
>>>
>>> On the testing side, I noticed the following:
>>>
>>> a) Doing a full comparison of p p > t t~ j [QCD] generated in single and
>> multicore, I noticed that the file 'SubProcesses/proc_characteristics' was
>> incorrect in mutlicore, since it contains:
>>> ninitial = 0
>>> nexternal = 0
>> hmm right, i was not checking this file with the unit test. I have added it to
>> the test, it should be fixed now (rev 315)
>>>
>>> b) This implementation broke process generation with another OLP than
>> MadLoop (i.e. Gosam).
>>> When selecting Gosam as OLP, the generation crashes at:
>>>
>>> FKSProcessError : Wrong type of argument for matrix_elements in function
>> write_lh_order.
>>>
>>> from
>>>
>>> [...]
>>> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-
>> memory/madgraph/iolibs/export_fks.py", line 1355, in
>> generate_virtuals_from_OLP
>>> self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)
>>> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-
>> memory/madgraph/iolibs/export_fks.py", line 1623, in write_lh_order
>>> 'matrix_elements in function write_lh_order.’)
>>>
>>> Of course we don't want to have any sort of parallelization on behalf of
>> GoSam, but normally here the only thing that is done is the writing of the
>> leshouches order. The actual generation with GoSam should take place at the
>> very end, basically after that the process is generated and the
>> parallelization is done. So if this is not too difficult to fix, I think it
>> would be worth doing so. Otherwise, explicitly forbid this mode when OLP is
>> not MadLoop, but that would be weird because there is no reason for doing so.
>>> Note that you don't need to have GoSam installed to test this I think,
>> because it crashes when writing the leshouches order file, so before actually
>> calling GoSam.
>>>
>> I have added a unit test also for this and reproduced your error. It should be
>> fixed now (rev 316), please check that i did not screw anything, as some
>> modifications to the code were needed (the unit tests pass after trival fixes)
>
> So when not using the low-mem option, now OLP = GoSam crashes with:
>
> AttributeError : 'MasterCmd' object has no attribute ‘born_processes_for_olp'

>
> at
>
> [...]
> File "/Users/valentin/Documents/Work/MG5/2.3.4-low-memory/madgraph/interface/amcatnlo_interface.py", line 557, in do_output
> self.born_processes_for_olp,self._export_dir,self.options['OLP'])
> AttributeError: 'MasterCmd' object has no attribute ‘born_processes_for_olp’
that should be ok now, please chec...

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

Hi Marco,

I will accept the merge even without the fix associate to this:

>> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)

But please do like that next time, since it is easier if we have access to the full user options rather that duplicate the information.

Cheers,

Olivier

review: Approve
319. By marco zaro

changed the arguments passed to FKSMultiProcess to options

Revision history for this message
marco zaro (marco-zaro) wrote :

Hi Olivier,
I have changed the arguments to options, please have a look. I have also removed a chunk of code from fks_base (10 lines starting with #check limitations of FKS) which did not make any sense and was never used.

Let me know if this is fine…
Cheers,

Marco
On 24 Feb 2016, at 13:31, Olivier Mattelaer <email address hidden> wrote:

> Review: Approve
>
> Hi Marco,
>
> I will accept the merge even without the fix associate to this:
>
>>> 3) concerning the ncores options, would not be better to pass to the class the self.options of the interface? Such that we do not have to multiply the optional argument when we to have access to any of the options define in that dictionary (this should have been done for OLP)
>
> But please do like that next time, since it is easier if we have access to the full user options rather that duplicate the information.
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

320. By marco zaro

imerged with latest 2.3.4

Revision history for this message
Valentin Hirschi (valentin-hirschi) wrote :

Hi Marco,

I tested the OLP=GoSam behavior and a couple other things and everything went smooth.

So please just confirm that you checked that there couldn't be any border effect on the modification of the function 'RunCardNLO.create_default_for_process(' from other parts of the code that might use it and also the fact that no deepcopy at all is used in the new low-mem mode.

Anyway, I assume that the above is indeed the case so that I already give my green light.

Great work,

Cheers

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

> Hi Olivier,
> I have changed the arguments to options, please have a look. I have also
> removed a chunk of code from fks_base (10 lines starting with #check
> limitations of FKS) which did not make any sense and was never used.

That's already much nicer. Thanks.

Actually my argument was to avoid to use the full options of MadGraph:
and therefore have instead of
self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,fks_options)
the following line:
self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,self.options)

Obviously, this force to have a slightly different behavior for
ncores_for_proc_gen
since that one is not part of the self.options dictionary.
The two advantages: all the configuration options are available and less memory consumption (since you have one less dictionary in memory).

But this is already so much nicer that I'm more than happy with your nice change.
So this comment was just to clarify my point and not to encourage you to perform any changes.

Thanks,

Olivier

Olivier

Revision history for this message
marco zaro (marco-zaro) wrote :

Ciao Olivier,
thanks for the comment
On 29 Feb 2016, at 11:50, Olivier Mattelaer <email address hidden> wrote:

>> Hi Olivier,
>> I have changed the arguments to options, please have a look. I have also
>> removed a chunk of code from fks_base (10 lines starting with #check
>> limitations of FKS) which did not make any sense and was never used.
>
> That's already much nicer. Thanks.
>
> Actually my argument was to avoid to use the full options of MadGraph:
> and therefore have instead of
> self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,fks_options)
> the following line:
> self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,self.options)
anyway the same options will have to be passed to the MultiProcess Instance, and many of them are not valid.
So sooner or later one has to create a new dictionary…
I will proceed with the merge

Cheers,

Marco
>
> Obviously, this force to have a slightly different behavior for
> ncores_for_proc_gen
> since that one is not part of the self.options dictionary.
> The two advantages: all the configuration options are available and less memory consumption (since you have one less dictionary in memory).
>
> But this is already so much nicer that I'm more than happy with your nice change.
> So this comment was just to clarify my point and not to encourage you to perform any changes.
>
> Thanks,
>
> Olivier
>
> Olivier
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.3.4-low-memory/+merge/286337
> Your team MadDevelopers is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.3.4-low-memory into lp:~maddevelopers/mg5amcnlo/2.3.4.

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt 2016-02-23 22:55:27 +0000
+++ UpdateNotes.txt 2016-02-24 13:54:50 +0000
@@ -1,6 +1,12 @@
1Update notes for MadGraph5_aMC@NLO (in reverse time order)1Update notes for MadGraph5_aMC@NLO (in reverse time order)
22
32.4.0(XX/XX/XX)32.4.0(XX/XX/XX)
4 MZ: new NLO generation mode. It is more efficient from the memory and CPU point of
5 view, in particular for high-multiplicity processes.
6 Many thanks to Josh Bendavid for his fundamental contribution for this.
7 The mode can be enabled with
8 > set low_mem_multicore_nlo_generation True
9 before generating the process.
4 VH: Interfaced MadLoop to Samurai and Ninja (the latter is now the default)10 VH: Interfaced MadLoop to Samurai and Ninja (the latter is now the default)
5 OM: Adding the possibility to use new syntax for tree-level processes:11 OM: Adding the possibility to use new syntax for tree-level processes:
6 QED==2 and QCD>2: The first allows to select exactly a power of the coupling (at amplitude level12 QED==2 and QCD>2: The first allows to select exactly a power of the coupling (at amplitude level
713
=== modified file 'madgraph/fks/fks_base.py'
--- madgraph/fks/fks_base.py 2015-10-20 10:56:04 +0000
+++ madgraph/fks/fks_base.py 2016-02-24 13:54:50 +0000
@@ -44,15 +44,18 @@
44 def default_setup(self):44 def default_setup(self):
45 """Default values for all properties"""45 """Default values for all properties"""
46 super(FKSMultiProcess, self).default_setup()46 super(FKSMultiProcess, self).default_setup()
47 self['real_amplitudes'] = diagram_generation.AmplitudeList()
48 self['pdgs'] = []
47 self['born_processes'] = FKSProcessList()49 self['born_processes'] = FKSProcessList()
48 if not 'OLP' in self.keys():50 if not 'OLP' in self.keys():
49 self['OLP'] = 'MadLoop'51 self['OLP'] = 'MadLoop'
52 self['ncores_for_proc_gen'] = 0
50 53
51 def get_sorted_keys(self):54 def get_sorted_keys(self):
52 """Return particle property names as a nicely sorted list."""55 """Return particle property names as a nicely sorted list."""
53 keys = super(FKSMultiProcess, self).get_sorted_keys()56 keys = super(FKSMultiProcess, self).get_sorted_keys()
54 keys += ['born_processes', 'real_amplitudes', 'real_pdgs', 'has_isr', 57 keys += ['born_processes', 'real_amplitudes', 'real_pdgs', 'has_isr',
55 'has_fsr', 'OLP']58 'has_fsr', 'OLP', 'ncores_for_proc_gen']
56 return keys59 return keys
5760
58 def filter(self, name, value):61 def filter(self, name, value):
@@ -77,10 +80,15 @@
77 if not isinstance(value,str):80 if not isinstance(value,str):
78 raise self.PhysicsObjectError, \81 raise self.PhysicsObjectError, \
79 "%s is not a valid string for OLP " % str(value)82 "%s is not a valid string for OLP " % str(value)
83
84 if name == 'ncores_for_proc_gen':
85 if not isinstance(value,int):
86 raise self.PhysicsObjectError, \
87 "%s is not a valid value for ncores_for_proc_gen " % str(value)
80 88
81 return super(FKSMultiProcess,self).filter(name, value)89 return super(FKSMultiProcess,self).filter(name, value)
82 90
83 def __init__(self, *arguments, **options):91 def __init__(self, procdef=None, options={}):
84 """Initializes the original multiprocess, then generates the amps for the 92 """Initializes the original multiprocess, then generates the amps for the
85 borns, then generate the born processes and the reals.93 borns, then generate the born processes and the reals.
86 Real amplitudes are stored in real_amplitudes according on the pdgs of their94 Real amplitudes are stored in real_amplitudes according on the pdgs of their
@@ -94,17 +102,25 @@
94 for logg in loggers_off:102 for logg in loggers_off:
95 logg.setLevel(logging.WARNING)103 logg.setLevel(logging.WARNING)
96 104
97 self['real_amplitudes'] = diagram_generation.AmplitudeList()105 # OLP option
98 self['pdgs'] = []106 olp='MadLoop'
99
100 if 'OLP' in options.keys():107 if 'OLP' in options.keys():
101 self['OLP']=options['OLP']108 olp = options['OLP']
102 del options['OLP']109 del options['OLP']
103110
111 ncores_for_proc_gen = 0
112 # ncores_for_proc_gen has the following meaning
113 # 0 : do things the old way
114 # > 0 use ncores_for_proc_gen
115 # -1 : use all cores
116 if 'ncores_for_proc_gen' in options.keys():
117 ncores_for_proc_gen = options['ncores_for_proc_gen']
118 del options['ncores_for_proc_gen']
104119
105 try:120 try:
106 # Now generating the borns for the first time.121 # Now generating the borns for the first time.
107 super(FKSMultiProcess, self).__init__(*arguments,**options)122 super(FKSMultiProcess, self).__init__(procdef, **options)
123
108 except diagram_generation.NoDiagramException as error:124 except diagram_generation.NoDiagramException as error:
109 # If no born, then this process most likely does not have any.125 # If no born, then this process most likely does not have any.
110 raise NoBornException, "Born diagrams could not be generated for the "+\126 raise NoBornException, "Born diagrams could not be generated for the "+\
@@ -113,17 +129,9 @@
113 " processes yet, but you can still use MadLoop if you want to "+\129 " processes yet, but you can still use MadLoop if you want to "+\
114 "only generate them."+\130 "only generate them."+\
115 " For this, use the 'virt=' mode, without multiparticle labels."131 " For this, use the 'virt=' mode, without multiparticle labels."
116 132
117 #check limitation of FKS133 self['OLP'] = olp
118 if arguments and isinstance(arguments, MG.Process):134 self['ncores_for_proc_gen'] = ncores_for_proc_gen
119 myprocdef = arguments[0]
120 misc.sprint( myprocdef.keys())
121 if myprocdef['perturbation_couplings']!=['QCD']:
122 raise InvalidCmd("FKS for reals only available in QCD for now, you asked %s" \
123 % ', '.join(myprocdef['perturbation_couplings']))
124 elif myprocdef.get_ninitial()==1:
125 raise InvalidCmd("At this stage aMC@NLO cannot handle decay process.\n"+\
126 " Only Leading Order (loop-induced and tree level) decay are supported.")
127 135
128 #check process definition(s):136 #check process definition(s):
129 # a process such as g g > g g will lead to real emissions 137 # a process such as g g > g g will lead to real emissions
@@ -164,58 +172,63 @@
164 'Process', ''),172 'Process', ''),
165 i + 1, len(amps)))173 i + 1, len(amps)))
166174
167 born = FKSProcess(amp)175 born = FKSProcess(amp, ncores_for_proc_gen = self['ncores_for_proc_gen'])
168 self['born_processes'].append(born)176 self['born_processes'].append(born)
169 born.generate_reals(self['pdgs'], self['real_amplitudes'])177 born.generate_reals(self['pdgs'], self['real_amplitudes'])
170178
171 born_pdg_list = [[l['id'] for l in born.born_proc['legs']] \179 if not self['ncores_for_proc_gen']:
172 for born in self['born_processes'] ]180 # old generation mode
173181
174 for born in self['born_processes']:182 born_pdg_list = [[l['id'] for l in born.born_proc['legs']] \
175 for real in born.real_amps:183 for born in self['born_processes'] ]
176 real.find_fks_j_from_i(born_pdg_list)184
177185 for born in self['born_processes']:
178 if amps:186 for real in born.real_amps:
179 if self['process_definitions'][0].get('NLO_mode') == 'all':187 real.find_fks_j_from_i(born_pdg_list)
180 self.generate_virtuals()188 if amps:
181 189 if self['process_definitions'][0].get('NLO_mode') == 'all':
182 elif not self['process_definitions'][0].get('NLO_mode') in ['all', 'real','LOonly']:190 self.generate_virtuals()
183 raise fks_common.FKSProcessError(\191
184 "Not a valid NLO_mode for a FKSMultiProcess: %s" % \192 elif not self['process_definitions'][0].get('NLO_mode') in ['all', 'real','LOonly']:
185 self['process_definitions'][0].get('NLO_mode'))193 raise fks_common.FKSProcessError(\
186194 "Not a valid NLO_mode for a FKSMultiProcess: %s" % \
187 # now get the total number of diagrams195 self['process_definitions'][0].get('NLO_mode'))
188 n_diag_born = sum([len(amp.get('diagrams')) 196
189 for amp in self.get_born_amplitudes()])197 # now get the total number of diagrams
190 n_diag_real = sum([len(amp.get('diagrams')) 198 n_diag_born = sum([len(amp.get('diagrams'))
191 for amp in self.get_real_amplitudes()])199 for amp in self.get_born_amplitudes()])
192 n_diag_virt = sum([len(amp.get('loop_diagrams')) 200 n_diag_real = sum([len(amp.get('diagrams'))
193 for amp in self.get_virt_amplitudes()])201 for amp in self.get_real_amplitudes()])
194202 n_diag_virt = sum([len(amp.get('loop_diagrams'))
195 if n_diag_virt == 0 and n_diag_real ==0 and \203 for amp in self.get_virt_amplitudes()])
196 not self['process_definitions'][0].get('NLO_mode') == 'LOonly':204
197 raise fks_common.FKSProcessError(205 if n_diag_virt == 0 and n_diag_real ==0 and \
198 'This process does not have any correction up to NLO in %s'\206 not self['process_definitions'][0].get('NLO_mode') == 'LOonly':
199 %','.join(perturbation))207 raise fks_common.FKSProcessError(
200208 'This process does not have any correction up to NLO in %s'\
201 logger.info(('Generated %d subprocesses with %d real emission diagrams, ' + \209 %','.join(perturbation))
202 '%d born diagrams and %d virtual diagrams') % \210
203 (len(self['born_processes']), n_diag_real, n_diag_born, n_diag_virt))211 logger.info(('Generated %d subprocesses with %d real emission diagrams, ' + \
204212 '%d born diagrams and %d virtual diagrams') % \
205 for i, logg in enumerate(loggers_off):213 (len(self['born_processes']), n_diag_real, n_diag_born, n_diag_virt))
206 logg.setLevel(old_levels[i])214
215 for i, logg in enumerate(loggers_off):
216 logg.setLevel(old_levels[i])
207217
208 self['has_isr'] = any([proc.isr for proc in self['born_processes']])218 self['has_isr'] = any([proc.isr for proc in self['born_processes']])
209 self['has_fsr'] = any([proc.fsr for proc in self['born_processes']])219 self['has_fsr'] = any([proc.fsr for proc in self['born_processes']])
210220
211 def add(self, other):221 def add(self, other):
212 """combines self and other, extending the lists of born/real amplitudes"""222 """combines self and other, extending the lists of born/real amplitudes"""
223 self['process_definitions'].extend(other['process_definitions'])
224 self['amplitudes'].extend(other['amplitudes'])
213 self['born_processes'].extend(other['born_processes'])225 self['born_processes'].extend(other['born_processes'])
214 self['real_amplitudes'].extend(other['real_amplitudes'])226 self['real_amplitudes'].extend(other['real_amplitudes'])
215 self['pdgs'].extend(other['pdgs'])227 self['pdgs'].extend(other['pdgs'])
216 self['has_isr'] = self['has_isr'] or other['has_isr']228 self['has_isr'] = self['has_isr'] or other['has_isr']
217 self['has_fsr'] = self['has_fsr'] or other['has_fsr']229 self['has_fsr'] = self['has_fsr'] or other['has_fsr']
218 self['OLP'] = other['OLP']230 self['OLP'] = other['OLP']
231 self['ncores_for_proc_gen'] = other['ncores_for_proc_gen']
219232
220 def get_born_amplitudes(self):233 def get_born_amplitudes(self):
221 """return an amplitudelist with the born amplitudes"""234 """return an amplitudelist with the born amplitudes"""
@@ -394,11 +407,16 @@
394 """The class for a FKS process. Starts from the born process and finds407 """The class for a FKS process. Starts from the born process and finds
395 all the possible splittings.""" 408 all the possible splittings."""
396 409
397 def __init__(self, start_proc = None, remove_reals = True):410 def __init__(self, start_proc = None, remove_reals = True, ncores_for_proc_gen=0):
398 """initialization: starts either from an amplitude or a process,411 """initialization: starts either from an amplitude or a process,
399 then init the needed variables.412 then init the needed variables.
400 remove_borns tells if the borns not needed for integration will be removed413 remove_borns tells if the borns not needed for integration will be removed
401 from the born list (mainly used for testing)"""414 from the born list (mainly used for testing)
415 ncores_for_proc_gen has the following meaning
416 0 : do things the old way
417 > 0 use ncores_for_proc_gen
418 -1 : use all cores
419 """
402 420
403 self.splittings = {}421 self.splittings = {}
404 self.reals = []422 self.reals = []
@@ -416,6 +434,7 @@
416 self.nincoming = 0434 self.nincoming = 0
417 self.virt_amp = None435 self.virt_amp = None
418 self.perturbation = 'QCD'436 self.perturbation = 'QCD'
437 self.ncores_for_proc_gen = ncores_for_proc_gen
419438
420 if not remove_reals in [True, False]:439 if not remove_reals in [True, False]:
421 raise fks_common.FKSProcessError(\440 raise fks_common.FKSProcessError(\
@@ -509,7 +528,6 @@
509 self.real_amps = real_amps528 self.real_amps = real_amps
510529
511530
512
513 def generate_reals(self, pdg_list, real_amp_list, combine=True): #test written531 def generate_reals(self, pdg_list, real_amp_list, combine=True): #test written
514 """For all the possible splittings, creates an FKSRealProcess.532 """For all the possible splittings, creates an FKSRealProcess.
515 It removes double counted configorations from the ones to integrates and533 It removes double counted configorations from the ones to integrates and
@@ -534,8 +552,9 @@
534 self.find_reals_to_integrate()552 self.find_reals_to_integrate()
535 if combine:553 if combine:
536 self.combine_real_amplitudes()554 self.combine_real_amplitudes()
537 self.generate_real_amplitudes(pdg_list, real_amp_list)555 if not self.ncores_for_proc_gen:
538 self.link_born_reals()556 self.generate_real_amplitudes(pdg_list, real_amp_list)
557 self.link_born_reals()
539558
540559
541 def link_born_reals(self):560 def link_born_reals(self):
542561
=== modified file 'madgraph/fks/fks_helas_objects.py'
--- madgraph/fks/fks_helas_objects.py 2015-10-01 16:00:08 +0000
+++ madgraph/fks/fks_helas_objects.py 2016-02-24 13:54:50 +0000
@@ -25,20 +25,183 @@
25import madgraph.fks.fks_base as fks_base25import madgraph.fks.fks_base as fks_base
26import madgraph.fks.fks_common as fks_common26import madgraph.fks.fks_common as fks_common
27import madgraph.loop.loop_helas_objects as loop_helas_objects27import madgraph.loop.loop_helas_objects as loop_helas_objects
28import madgraph.loop.loop_diagram_generation as loop_diagram_generation
28import copy29import copy
29import logging30import logging
30import array31import array
32import multiprocessing
33import signal
34import tempfile
35import cPickle
36import itertools
37import os
3138
32logger = logging.getLogger('madgraph.fks_helas_objects')39logger = logging.getLogger('madgraph.fks_helas_objects')
3340
3441
42#functions to be used in the ncores_for_proc_gen mode
43def async_generate_real(args):
44 i = args[0]
45 real_amp = args[1]
46
47 #amplitude generation
48 amplitude = real_amp.generate_real_amplitude()
49 helasreal = helas_objects.HelasMatrixElement(amplitude)
50 logger.info('Generating real %s' % \
51 real_amp.process.nice_string(print_weighted=False).replace('Process', 'process'))
52
53 # Keep track of already generated color objects, to reuse as
54 # much as possible
55 list_colorize = []
56 list_color_basis = []
57 list_color_matrices = []
58
59 # Now this keeps track of the color matrices created from the loop-born
60 # color basis. Keys are 2-tuple with the index of the loop and born basis
61 # in the list above and the value is the resulting matrix.
62 dict_loopborn_matrices = {}
63 # The dictionary below is simply a container for convenience to be
64 # passed to the function process_color.
65 color_information = { 'list_colorize' : list_colorize,
66 'list_color_basis' : list_color_basis,
67 'list_color_matrices' : list_color_matrices,
68 'dict_loopborn_matrices' : dict_loopborn_matrices}
69
70 helas_objects.HelasMultiProcess.process_color(helasreal,color_information)
71
72 outdata = [amplitude,helasreal]
73
74 output = tempfile.NamedTemporaryFile(delete = False)
75 cPickle.dump(outdata,output,protocol=2)
76 output.close()
77
78 return [output.name,helasreal.get_num_configs(),helasreal.get_nexternal_ninitial()[0]]
79
80
81def async_generate_born(args):
82 i = args[0]
83 born = args[1]
84 born_pdg_list = args[2]
85 loop_orders = args[3]
86 pdg_list = args[4]
87 loop_optimized = args[5]
88 OLP = args[6]
89 realmapout = args[7]
90
91 logger.info('Generating born %s' % \
92 born.born_proc.nice_string(print_weighted=False).replace('Process', 'process'))
93
94 #load informations on reals from temp files
95 helasreal_list = []
96 for amp in born.real_amps:
97 idx = pdg_list.index(amp.pdgs)
98 infilename = realmapout[idx]
99 infile = open(infilename,'rb')
100 realdata = cPickle.load(infile)
101 infile.close()
102 amp.amplitude = realdata[0]
103 helasreal_list.append(realdata[1])
104
105 born.link_born_reals()
106
107 for amp in born.real_amps:
108 amp.find_fks_j_from_i(born_pdg_list)
109
110 # generate the virtuals if needed
111 has_loops = False
112 if born.born_proc.get('NLO_mode') == 'all' and OLP == 'MadLoop':
113 myproc = copy.copy(born.born_proc)
114 # take the orders that are actually used by the matrix element
115 myproc['orders'] = loop_orders
116 myproc['legs'] = fks_common.to_legs(copy.copy(myproc['legs']))
117 myamp = loop_diagram_generation.LoopAmplitude(myproc)
118 if myamp.get('diagrams'):
119 has_loops = True
120 born.virt_amp = myamp
121
122 helasfull = FKSHelasProcess(born, helasreal_list,
123 loop_optimized = loop_optimized,
124 decay_ids=[],
125 gen_color=False)
126
127 processes = helasfull.born_matrix_element.get('processes')
128
129 metag = helas_objects.IdentifyMETag.create_tag(helasfull.born_matrix_element.get('base_amplitude'))
130
131 outdata = helasfull
132
133 output = tempfile.NamedTemporaryFile(delete = False)
134 cPickle.dump(outdata,output,protocol=2)
135 output.close()
136
137 return [output.name,metag,has_loops,processes]
138
139
140def async_finalize_matrix_elements(args):
141
142 i = args[0]
143 mefile = args[1]
144 duplist = args[2]
145
146 infile = open(mefile,'rb')
147 me = cPickle.load(infile)
148 infile.close()
149
150 #set unique id based on position in unique me list
151 me.get('processes')[0].set('uid', i)
152
153 # Always create an empty color basis, and the
154 # list of raw colorize objects (before
155 # simplification) associated with amplitude
156 col_basis = color_amp.ColorBasis()
157 new_amp = me.born_matrix_element.get_base_amplitude()
158 me.born_matrix_element.set('base_amplitude', new_amp)
159 colorize_obj = col_basis.create_color_dict_list(new_amp)
160
161 col_basis.build()
162 col_matrix = color_amp.ColorMatrix(col_basis)
163
164 me.born_matrix_element.set('color_basis',col_basis)
165 me.born_matrix_element.set('color_matrix',col_matrix)
166
167 for iother,othermefile in enumerate(duplist):
168 infileother = open(othermefile,'rb')
169 otherme = cPickle.load(infileother)
170 infileother.close()
171 me.add_process(otherme)
172
173 me.set_color_links()
174
175 initial_states=[]
176 for fksreal in me.real_processes:
177 # Pick out all initial state particles for the two beams
178 initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \
179 p in fksreal.matrix_element.get('processes')))))
180
181 if me.virt_matrix_element:
182 has_virtual = True
183 else:
184 has_virtual = False
185
186 #data to write to file
187 outdata = me
188
189 output = tempfile.NamedTemporaryFile(delete = False)
190 cPickle.dump(outdata,output,protocol=2)
191 output.close()
192
193 #data to be returned to parent process (filename plus small objects only)
194 return [output.name,initial_states,me.get_used_lorentz(),me.get_used_couplings(),has_virtual]
195
196
35class FKSHelasMultiProcess(helas_objects.HelasMultiProcess):197class FKSHelasMultiProcess(helas_objects.HelasMultiProcess):
36 """class to generate the helas calls for a FKSMultiProcess"""198 """class to generate the helas calls for a FKSMultiProcess"""
37199
38 def get_sorted_keys(self):200 def get_sorted_keys(self):
39 """Return particle property names as a nicely sorted list."""201 """Return particle property names as a nicely sorted list."""
40 keys = super(FKSHelasMultiProcess, self).get_sorted_keys()202 keys = super(FKSHelasMultiProcess, self).get_sorted_keys()
41 keys += ['real_matrix_elements', ['has_isr'], ['has_fsr']]203 keys += ['real_matrix_elements', ['has_isr'], ['has_fsr'],
204 'used_lorentz', 'used_couplings', 'max_configs', 'max_particles', 'processes']
42 return keys205 return keys
43206
44 def filter(self, name, value):207 def filter(self, name, value):
@@ -61,23 +224,190 @@
61224
62 self.loop_optimized = loop_optimized225 self.loop_optimized = loop_optimized
63226
64 # generate the real ME's if they are needed.227 self['used_lorentz'] = []
65 # note that it may not be always the case, e.g. it the NLO_mode is LOonly228 self['used_couplings'] = []
66 if fksmulti['real_amplitudes']:229 self['processes'] = []
67 logger.info('Generating real emission matrix-elements...')230
68 self['real_matrix_elements'] = self.generate_matrix_elements(231 self['max_particles'] = -1
69 copy.copy(fksmulti['real_amplitudes']), combine_matrix_elements = False)232 self['max_configs'] = -1
70 else:233
71 self['real_matrix_elements'] = helas_objects.HelasMatrixElementList()234 if not fksmulti['ncores_for_proc_gen']:
72235 # generate the real ME's if they are needed.
73 self['matrix_elements'] = self.generate_matrix_elements_fks(236 # note that it may not be always the case, e.g. it the NLO_mode is LOonly
74 fksmulti, 237 if fksmulti['real_amplitudes']:
75 gen_color, decay_ids)238 logger.info('Generating real emission matrix-elements...')
76 self['initial_states']=[]239 self['real_matrix_elements'] = self.generate_matrix_elements(
240 copy.copy(fksmulti['real_amplitudes']), combine_matrix_elements = False)
241 else:
242 self['real_matrix_elements'] = helas_objects.HelasMatrixElementList()
243
244 self['matrix_elements'] = self.generate_matrix_elements_fks(
245 fksmulti,
246 gen_color, decay_ids)
247 self['initial_states']=[]
248 self['has_loops'] = len(self.get_virt_matrix_elements()) > 0
249
250 else:
251 self['has_loops'] = False
252 #more efficient generation
253 born_procs = fksmulti.get('born_processes')
254 born_pdg_list = [[l['id'] for l in born.born_proc['legs']] \
255 for born in born_procs ]
256 loop_orders = {}
257 for born in born_procs:
258 for coup, val in fks_common.find_orders(born.born_amp).items():
259 try:
260 loop_orders[coup] = max([loop_orders[coup], val])
261 except KeyError:
262 loop_orders[coup] = val
263 pdg_list = []
264 real_amp_list = []
265 for born in born_procs:
266 for amp in born.real_amps:
267 if not pdg_list.count(amp.pdgs):
268 pdg_list.append(amp.pdgs)
269 real_amp_list.append(amp)
270
271 #generating and store in tmp files all output corresponding to each real_amplitude
272 real_out_list = []
273 realmapin = []
274 for i,real_amp in enumerate(real_amp_list):
275 realmapin.append([i,real_amp])
276
277 # start the pool instance with a signal instance to catch ctr+c
278 original_sigint_handler = signal.signal(signal.SIGINT, signal.SIG_IGN)
279 if fksmulti['ncores_for_proc_gen'] < 0: # use all cores
280 pool = multiprocessing.Pool(maxtasksperchild=1)
281 else:
282 pool = multiprocessing.Pool(processes=fksmulti['ncores_for_proc_gen'],maxtasksperchild=1)
283 signal.signal(signal.SIGINT, original_sigint_handler)
284
285 logger.info('Generating real matrix elements...')
286 try:
287 # the very large timeout passed to get is to be able to catch
288 # KeyboardInterrupts
289 realmapout = pool.map_async(async_generate_real,realmapin).get(9999999)
290 except KeyboardInterrupt:
291 pool.terminate()
292 raise KeyboardInterrupt
293
294 realmapfiles = []
295 for realout in realmapout:
296 realmapfiles.append(realout[0])
297
298 logger.info('Generating born and virtual matrix elements...')
299 #now loop over born and consume reals, generate virtuals
300 bornmapin = []
301 OLP=fksmulti['OLP']
302 for i,born in enumerate(born_procs):
303 bornmapin.append([i,born,born_pdg_list,loop_orders,pdg_list,loop_optimized,OLP,realmapfiles])
304
305 try:
306 bornmapout = pool.map_async(async_generate_born,bornmapin).get(9999999)
307 except KeyboardInterrupt:
308 pool.terminate()
309 raise KeyboardInterrupt
310
311 #remove real temp files
312 for realtmp in realmapout:
313 os.remove(realtmp[0])
314
315 logger.info('Collecting infos and finalizing matrix elements...')
316 unique_me_list = []
317 duplicate_me_lists = []
318 for bornout in bornmapout:
319 mefile = bornout[0]
320 metag = bornout[1]
321 has_loops = bornout[2]
322 self['has_loops'] = self['has_loops'] or has_loops
323 processes = bornout[3]
324 self['processes'].extend(processes)
325 unique = True
326 for ime2,bornout2 in enumerate(unique_me_list):
327 mefile2 = bornout2[0]
328 metag2 = bornout2[1]
329 if metag==metag2:
330 duplicate_me_lists[ime2].append(mefile)
331 unique = False
332 break;
333 if unique:
334 unique_me_list.append(bornout)
335 duplicate_me_lists.append([])
336
337 memapin = []
338 for i,bornout in enumerate(unique_me_list):
339 mefile = bornout[0]
340 memapin.append([i,mefile, duplicate_me_lists[i]])
341
342 try:
343 memapout = pool.map_async(async_finalize_matrix_elements,memapin).get(9999999)
344 except KeyboardInterrupt:
345 pool.terminate()
346 raise KeyboardInterrupt
347
348 #remove born+virtual temp files
349 for bornout in bornmapout:
350 mefile = bornout[0]
351 os.remove(mefile)
352
353 pool.close()
354 pool.join()
355
356 #set final list of matrix elements (paths to temp files)
357 matrix_elements = []
358 for meout in memapout:
359 matrix_elements.append(meout[0])
360
361 self['matrix_elements']=matrix_elements
362
363 #cache information needed for output which will not be available from
364 #the matrix elements later
365 initial_states = []
366 for meout in memapout:
367 me_initial_states = meout[1]
368 for state in me_initial_states:
369 initial_states.append(state)
370
371 # remove doubles from the list
372 checked = []
373 for e in initial_states:
374 if e not in checked:
375 checked.append(e)
376 initial_states=checked
377
378 self['initial_states']=initial_states
379
380 helas_list = []
381 for meout in memapout:
382 helas_list.extend(meout[2])
383 self['used_lorentz']=list(set(helas_list))
384
385 coupling_list = []
386 for meout in memapout:
387 coupling_list.extend([c for l in meout[3] for c in l])
388 self['used_couplings'] = list(set(coupling_list))
389
390 has_virtuals = False
391 for meout in memapout:
392 if meout[4]:
393 has_virtuals = True
394 break
395 self['has_virtuals'] = has_virtuals
396
397 configs_list = []
398 for meout in realmapout:
399 configs_list.append(meout[1])
400 self['max_configs'] = max(configs_list)
401
402 nparticles_list = []
403 for meout in realmapout:
404 nparticles_list.append(meout[2])
405 self['max_particles'] = max(nparticles_list)
77406
78 self['has_isr'] = fksmulti['has_isr']407 self['has_isr'] = fksmulti['has_isr']
79 self['has_fsr'] = fksmulti['has_fsr']408 self['has_fsr'] = fksmulti['has_fsr']
80 self['has_loops'] = len(self.get_virt_matrix_elements()) > 0 409
410 logger.info('... Done')
81411
82 for i, logg in enumerate(loggers_off):412 for i, logg in enumerate(loggers_off):
83 logg.setLevel(old_levels[i])413 logg.setLevel(old_levels[i])
@@ -85,19 +415,66 @@
85 def get_used_lorentz(self):415 def get_used_lorentz(self):
86 """Return a list of (lorentz_name, conjugate, outgoing) with416 """Return a list of (lorentz_name, conjugate, outgoing) with
87 all lorentz structures used by this HelasMultiProcess."""417 all lorentz structures used by this HelasMultiProcess."""
88 helas_list = []418
89 for me in self.get('matrix_elements'):419 if not self['used_lorentz']:
90 helas_list.extend(me.get_used_lorentz())420 helas_list = []
91 return list(set(helas_list))421 for me in self.get('matrix_elements'):
422 helas_list.extend(me.get_used_lorentz())
423 self['used_lorentz'] = list(set(helas_list))
424
425 return self['used_lorentz']
426
92427
93 def get_used_couplings(self):428 def get_used_couplings(self):
94 """Return a list with all couplings used by this429 """Return a list with all couplings used by this
95 HelasMatrixElement."""430 HelasMatrixElement."""
96 coupling_list = []431
97 for me in self.get('matrix_elements'):432 if not self['used_couplings']:
98 coupling_list.extend([c for l in me.get_used_couplings() for c in l])433 coupling_list = []
99 return list(set(coupling_list))434 for me in self.get('matrix_elements'):
435 coupling_list.extend([c for l in me.get_used_couplings() for c in l])
436 self['used_couplings'] = list(set(coupling_list))
437
438 return self['used_couplings']
439
440
441 def get_processes(self):
442 """Return a list with all couplings used by this
443 HelasMatrixElement."""
444
445 if not self['processes']:
446 process_list = []
447 for me in self.get('matrix_elements'):
448 process_list.extend(me.born_matrix_element.get('processes'))
449 self['processes'] = process_list
450
451 return self['processes']
452
453
454 def get_max_configs(self):
455 """Return max_configs"""
456
457 if self['max_configs'] < 0:
458 try:
459 self['max_configs'] = max([me.get_num_configs() \
460 for me in self['real_matrix_elements']])
461 except ValueError:
462 self['max_configs'] = max([me.born_matrix_element.get_num_configs() \
463 for me in self['matrix_elements']])
464
465 return self['max_configs']
466
467
468 def get_max_particles(self):
469 """Return max_paricles"""
470
471 if self['max_particles'] < 0:
472 self['max_particles'] = max([me.get_nexternal_ninitial()[0] \
473 for me in self['matrix_elements']])
474
475 return self['max_particles']
100 476
477
101 def get_matrix_elements(self):478 def get_matrix_elements(self):
102 """Extract the list of matrix elements"""479 """Extract the list of matrix elements"""
103 return self.get('matrix_elements') 480 return self.get('matrix_elements')
@@ -240,17 +617,32 @@
240 self.perturbation = fksproc.perturbation617 self.perturbation = fksproc.perturbation
241 real_amps_new = []618 real_amps_new = []
242 # combine for example u u~ > t t~ and d d~ > t t~619 # combine for example u u~ > t t~ and d d~ > t t~
243 for proc in fksproc.real_amps:620 if fksproc.ncores_for_proc_gen:
244 fksreal_me = FKSHelasRealProcess(proc, real_me_list, real_amp_list, **opts)621 # new NLO (multicore) generation mode
245 try:622 for real_me, proc in itertools.izip(real_me_list,fksproc.real_amps):
246 other = self.real_processes[self.real_processes.index(fksreal_me)]623 fksreal_me = FKSHelasRealProcess(proc, real_me, **opts)
247 other.matrix_element.get('processes').extend(\624 try:
248 fksreal_me.matrix_element.get('processes') )625 other = self.real_processes[self.real_processes.index(fksreal_me)]
249 except ValueError:626 other.matrix_element.get('processes').extend(\
250 if fksreal_me.matrix_element.get('processes') and \627 fksreal_me.matrix_element.get('processes') )
251 fksreal_me.matrix_element.get('diagrams'):628 except ValueError:
252 self.real_processes.append(fksreal_me)629 if fksreal_me.matrix_element.get('processes') and \
253 real_amps_new.append(proc)630 fksreal_me.matrix_element.get('diagrams'):
631 self.real_processes.append(fksreal_me)
632 real_amps_new.append(proc)
633 else:
634 #old mode
635 for proc in fksproc.real_amps:
636 fksreal_me = FKSHelasRealProcess(proc, real_me_list, real_amp_list, **opts)
637 try:
638 other = self.real_processes[self.real_processes.index(fksreal_me)]
639 other.matrix_element.get('processes').extend(\
640 fksreal_me.matrix_element.get('processes') )
641 except ValueError:
642 if fksreal_me.matrix_element.get('processes') and \
643 fksreal_me.matrix_element.get('diagrams'):
644 self.real_processes.append(fksreal_me)
645 real_amps_new.append(proc)
254 fksproc.real_amps = real_amps_new646 fksproc.real_amps = real_amps_new
255 if fksproc.virt_amp:647 if fksproc.virt_amp:
256 self.virt_matrix_element = \648 self.virt_matrix_element = \
@@ -410,23 +802,35 @@
410 self.fks_infos = fksrealproc.fks_infos802 self.fks_infos = fksrealproc.fks_infos
411 self.is_to_integrate = fksrealproc.is_to_integrate803 self.is_to_integrate = fksrealproc.is_to_integrate
412804
413 if len(real_me_list) != len(real_amp_list):805 # real_me_list is a list in the old NLO generation mode;
806 # in the new one it is a matrix element
807 if type(real_me_list) == list and len(real_me_list) != len(real_amp_list):
414 raise fks_common.FKSProcessError(808 raise fks_common.FKSProcessError(
415 'not same number of amplitudes and matrix elements: %d, %d' % \809 'not same number of amplitudes and matrix elements: %d, %d' % \
416 (len(real_amp_list), len(real_me_list)))810 (len(real_amp_list), len(real_me_list)))
417 if real_me_list and real_amp_list:811 if type(real_me_list) == list and real_me_list and real_amp_list:
418 self.matrix_element = copy.deepcopy(real_me_list[real_amp_list.index(fksrealproc.amplitude)])812 self.matrix_element = copy.deepcopy(real_me_list[real_amp_list.index(fksrealproc.amplitude)])
419 self.matrix_element['processes'] = copy.deepcopy(self.matrix_element['processes'])813 self.matrix_element['processes'] = copy.deepcopy(self.matrix_element['processes'])
814
815 elif type(real_me_list) == helas_objects.HelasMatrixElement:
816 #new NLO generation mode
817 self.matrix_element = real_me_list
818
420 else:819 else:
421 logger.info('generating matrix element...')820
422 self.matrix_element = helas_objects.HelasMatrixElement(821 if real_me_list and real_amp_list:
423 fksrealproc.amplitude, **opts)822 self.matrix_element = copy.deepcopy(real_me_list[real_amp_list.index(fksrealproc.amplitude)])
424 #generate the color for the real823 self.matrix_element['processes'] = copy.deepcopy(self.matrix_element['processes'])
425 self.matrix_element.get('color_basis').build(824 else:
426 self.matrix_element.get('base_amplitude'))825 logger.info('generating matrix element...')
427 self.matrix_element.set('color_matrix',826 self.matrix_element = helas_objects.HelasMatrixElement(
428 color_amp.ColorMatrix(827 fksrealproc.amplitude, **opts)
429 self.matrix_element.get('color_basis')))828 #generate the color for the real
829 self.matrix_element.get('color_basis').build(
830 self.matrix_element.get('base_amplitude'))
831 self.matrix_element.set('color_matrix',
832 color_amp.ColorMatrix(
833 self.matrix_element.get('color_basis')))
430 #self.fks_j_from_i = fksrealproc.find_fks_j_from_i()834 #self.fks_j_from_i = fksrealproc.find_fks_j_from_i()
431 self.fks_j_from_i = fksrealproc.fks_j_from_i835 self.fks_j_from_i = fksrealproc.fks_j_from_i
432836
433837
=== modified file 'madgraph/interface/amcatnlo_interface.py'
--- madgraph/interface/amcatnlo_interface.py 2015-10-20 14:39:50 +0000
+++ madgraph/interface/amcatnlo_interface.py 2016-02-24 13:54:50 +0000
@@ -24,6 +24,13 @@
24import optparse24import optparse
25import subprocess25import subprocess
26import shutil26import shutil
27import multiprocessing
28import signal
29import tempfile
30import itertools
31import os
32import cPickle
33
2734
28import madgraph35import madgraph
29from madgraph import MG4DIR, MG5DIR, MadGraph5Error36from madgraph import MG4DIR, MG5DIR, MadGraph5Error
@@ -52,6 +59,37 @@
52logger = logging.getLogger('cmdprint') # -> stdout59logger = logging.getLogger('cmdprint') # -> stdout
53logger_stderr = logging.getLogger('fatalerror') # ->stderr60logger_stderr = logging.getLogger('fatalerror') # ->stderr
5461
62# a new function for the improved NLO generation
63glob_directories_map = []
64def generate_directories_fks_async(i):
65
66 arglist = glob_directories_map[i]
67
68 curr_exporter = arglist[0]
69 mefile = arglist[1]
70 curr_fortran_model = arglist[2]
71 ime = arglist[3]
72 nme = arglist[4]
73 path = arglist[5]
74 olpopts = arglist[6]
75
76 infile = open(mefile,'rb')
77 me = cPickle.load(infile)
78 infile.close()
79
80 calls = curr_exporter.generate_directories_fks(me, curr_fortran_model, ime, nme, path, olpopts)
81 nexternal = curr_exporter.proc_characteristic['nexternal']
82 ninitial = curr_exporter.proc_characteristic['ninitial']
83 processes = me.born_matrix_element.get('processes')
84
85 #only available after export has been done, so has to be returned from here
86 max_loop_vertex_rank = -99
87 if me.virt_matrix_element:
88 max_loop_vertex_rank = me.virt_matrix_element.get_max_loop_vertex_rank()
89
90 return [calls, curr_exporter.fksdirs, max_loop_vertex_rank, ninitial, nexternal, processes]
91
92
55class CheckFKS(mg_interface.CheckValidForCmd):93class CheckFKS(mg_interface.CheckValidForCmd):
5694
5795
@@ -431,16 +469,30 @@
431# 469#
432# raise self.InvalidCmd("FKS for reals only available in QCD for now, you asked %s" \470# raise self.InvalidCmd("FKS for reals only available in QCD for now, you asked %s" \
433# % ', '.join(myprocdef['perturbation_couplings']))471# % ', '.join(myprocdef['perturbation_couplings']))
472 ##
473
474 # if the new nlo process generation mode is enabled, the number of cores to be
475 # used has to be passed
476 # ncores_for_proc_gen has the following meaning
477 # 0 : do things the old way
478 # > 0 use ncores_for_proc_gen
479 # -1 : use all cores
480 if self.options['low_mem_multicore_nlo_generation']:
481 if self.options['nb_core']:
482 self.ncores_for_proc_gen = int(self.options['nb_core'])
483 else:
484 self.ncores_for_proc_gen = -1
485 else:
486 self.ncores_for_proc_gen = 0
487
488 # this is the options dictionary to pass to the FKSMultiProcess
489 fks_options = {'OLP': self.options['OLP'],
490 'ignore_six_quark_processes': self.options['ignore_six_quark_processes'],
491 'ncores_for_proc_gen': self.ncores_for_proc_gen}
434 try:492 try:
435 self._fks_multi_proc.add(fks_base.FKSMultiProcess(myprocdef,493 self._fks_multi_proc.add(fks_base.FKSMultiProcess(myprocdef,fks_options))
436 collect_mirror_procs,
437 ignore_six_quark_processes,
438 OLP=self.options['OLP']))
439 except AttributeError: 494 except AttributeError:
440 self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,495 self._fks_multi_proc = fks_base.FKSMultiProcess(myprocdef,fks_options)
441 collect_mirror_procs,
442 ignore_six_quark_processes,
443 OLP=self.options['OLP'])
444496
445497
446 def do_output(self, line):498 def do_output(self, line):
@@ -498,7 +550,7 @@
498 # Generate the virtuals if from OLP550 # Generate the virtuals if from OLP
499 if self.options['OLP']!='MadLoop':551 if self.options['OLP']!='MadLoop':
500 self._curr_exporter.generate_virtuals_from_OLP(552 self._curr_exporter.generate_virtuals_from_OLP(
501 self._curr_matrix_elements,self._export_dir,self.options['OLP'])553 self.born_processes_for_olp,self._export_dir,self.options['OLP'])
502 554
503 # Remember that we have done export555 # Remember that we have done export
504 self._done_export = (self._export_dir, self._export_format)556 self._done_export = (self._export_dir, self._export_format)
@@ -531,41 +583,50 @@
531 self._fks_multi_proc, 583 self._fks_multi_proc,
532 loop_optimized= self.options['loop_optimized_output'])584 loop_optimized= self.options['loop_optimized_output'])
533 585
534 ndiags = sum([len(me.get('diagrams')) for \586 if not self.options['low_mem_multicore_nlo_generation']:
535 me in self._curr_matrix_elements.\587 # generate the code the old way
536 get_matrix_elements()])588 ndiags = sum([len(me.get('diagrams')) for \
537 # assign a unique id number to all process and589 me in self._curr_matrix_elements.\
538 # generate a list of possible PDF combinations590 get_matrix_elements()])
539 uid = 0 591 # assign a unique id number to all process and
540 initial_states=[]592 # generate a list of possible PDF combinations
541 for me in self._curr_matrix_elements.get_matrix_elements():593 uid = 0
542 uid += 1 # update the identification number594 initial_states=[]
543 me.get('processes')[0].set('uid', uid)595 for me in self._curr_matrix_elements.get_matrix_elements():
544 try:596 uid += 1 # update the identification number
545 initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \597 me.get('processes')[0].set('uid', uid)
546 p in me.born_matrix_element.get('processes')))))
547 except IndexError:
548 initial_states.append(sorted(list(set((p.get_initial_pdg(1)) for \
549 p in me.born_matrix_element.get('processes')))))
550
551 for fksreal in me.real_processes:
552 # Pick out all initial state particles for the two beams
553 try:598 try:
554 initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \599 initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \
555 p in fksreal.matrix_element.get('processes')))))600 p in me.born_matrix_element.get('processes')))))
556 except IndexError:601 except IndexError:
557 initial_states.append(sorted(list(set((p.get_initial_pdg(1)) for \602 initial_states.append(sorted(list(set((p.get_initial_pdg(1)) for \
558 p in fksreal.matrix_element.get('processes')))))603 p in me.born_matrix_element.get('processes')))))
559
560 604
561 # remove doubles from the list605 for fksreal in me.real_processes:
562 checked = []606 # Pick out all initial state particles for the two beams
563 for e in initial_states:607 try:
564 if e not in checked:608 initial_states.append(sorted(list(set((p.get_initial_pdg(1),p.get_initial_pdg(2)) for \
565 checked.append(e)609 p in fksreal.matrix_element.get('processes')))))
566 initial_states=checked610 except IndexError:
567611 initial_states.append(sorted(list(set((p.get_initial_pdg(1)) for \
568 self._curr_matrix_elements.set('initial_states',initial_states)612 p in fksreal.matrix_element.get('processes')))))
613
614
615 # remove doubles from the list
616 checked = []
617 for e in initial_states:
618 if e not in checked:
619 checked.append(e)
620 initial_states=checked
621
622 self._curr_matrix_elements.set('initial_states',initial_states)
623
624 else:
625 #new NLO generation
626 if self._curr_matrix_elements['has_loops']:
627 self._curr_exporter.opt['mp'] = True
628 self._curr_exporter.model = self._curr_model
629 ndiags = 0
569630
570 cpu_time2 = time.time()631 cpu_time2 = time.time()
571 return ndiags, cpu_time2 - cpu_time1632 return ndiags, cpu_time2 - cpu_time1
@@ -586,23 +647,84 @@
586 for charac in ['has_isr', 'has_fsr', 'has_loops']:647 for charac in ['has_isr', 'has_fsr', 'has_loops']:
587 proc_charac[charac] = self._curr_matrix_elements[charac]648 proc_charac[charac] = self._curr_matrix_elements[charac]
588649
650 # prepare for the generation
651 # glob_directories_map is for the new NLO generation
652 global glob_directories_map
653 glob_directories_map = []
589654
655 self.born_processes_for_olp = []
590 for ime, me in \656 for ime, me in \
591 enumerate(self._curr_matrix_elements.get('matrix_elements')):657 enumerate(self._curr_matrix_elements.get('matrix_elements')):
592 #me is a FKSHelasProcessFromReals658 if not self.options['low_mem_multicore_nlo_generation']:
593 calls = calls + \659 #me is a FKSHelasProcessFromReals
594 self._curr_exporter.generate_directories_fks(me, 660 calls = calls + \
595 self._curr_fortran_model, 661 self._curr_exporter.generate_directories_fks(me,
596 ime, len(self._curr_matrix_elements.get('matrix_elements')), 662 self._curr_fortran_model,
597 path,self.options['OLP'])663 ime, len(self._curr_matrix_elements.get('matrix_elements')),
598 self._fks_directories.extend(self._curr_exporter.fksdirs)664 path,self.options['OLP'])
665 self._fks_directories.extend(self._curr_exporter.fksdirs)
666 self.born_processes_for_olp.append(me.born_matrix_element.get('processes')[0])
667 else:
668 glob_directories_map.append(\
669 [self._curr_exporter, me, self._curr_fortran_model,
670 ime, len(self._curr_matrix_elements.get('matrix_elements')),
671 path, self.options['OLP']])
672
673 if self.options['low_mem_multicore_nlo_generation']:
674 # start the pool instance with a signal instance to catch ctr+c
675 logger.info('Writing directories...')
676 original_sigint_handler = signal.signal(signal.SIGINT, signal.SIG_IGN)
677 if self.ncores_for_proc_gen < 0: # use all cores
678 pool = multiprocessing.Pool(maxtasksperchild=1)
679 else:
680 pool = multiprocessing.Pool(processes=self.ncores_for_proc_gen,maxtasksperchild=1)
681 signal.signal(signal.SIGINT, original_sigint_handler)
682 try:
683 # the very large timeout passed to get is to be able to catch
684 # KeyboardInterrupts
685 diroutputmap = pool.map_async(generate_directories_fks_async,
686 range(len(glob_directories_map))).get(9999999)
687 except KeyboardInterrupt:
688 pool.terminate()
689 raise KeyboardInterrupt
690
691 pool.close()
692 pool.join()
693
694 #clean up tmp files containing final matrix elements
695 for mefile in self._curr_matrix_elements.get('matrix_elements'):
696 os.remove(mefile)
697
698 for charac in ['nexternal', 'ninitial']:
699 proc_charac[charac] = self._curr_exporter.proc_characteristic[charac]
700 # ninitial and nexternal
701 proc_charac['nexternal'] = max([diroutput[4] for diroutput in diroutputmap])
702 ninitial_set = set([diroutput[3] for diroutput in diroutputmap])
703 if len(ninitial_set) != 1:
704 raise MadGraph5Error, ("Invalid ninitial values: %s" % ' ,'.join(list(ninitial_set)))
705 proc_charac['ninitial'] = list(ninitial_set)[0]
706
707 self.born_processes = []
708 self.born_processes_for_olp = []
709 max_loop_vertex_ranks = []
710
711 for diroutput in diroutputmap:
712 calls = calls + diroutput[0]
713 self._fks_directories.extend(diroutput[1])
714 max_loop_vertex_ranks.append(diroutput[2])
715 self.born_processes.extend(diroutput[5])
716 self.born_processes_for_olp.append(diroutput[5][0])
717
718 else:
719 max_loop_vertex_ranks = [me.get_max_loop_vertex_rank() for \
720 me in self._curr_matrix_elements.get_virt_matrix_elements()]
721
599 card_path = os.path.join(path, os.path.pardir, 'SubProcesses', \722 card_path = os.path.join(path, os.path.pardir, 'SubProcesses', \
600 'procdef_mg5.dat')723 'procdef_mg5.dat')
601 724
602 if self.options['loop_optimized_output'] and \725 if self.options['loop_optimized_output'] and \
603 len(self._curr_matrix_elements.get_virt_matrix_elements()) > 0:726 len(max_loop_vertex_ranks) > 0:
604 self._curr_exporter.write_coef_specs_file(\727 self._curr_exporter.write_coef_specs_file(max_loop_vertex_ranks)
605 self._curr_matrix_elements.get_virt_matrix_elements())
606 if self._generate_info:728 if self._generate_info:
607 self._curr_exporter.write_procdef_mg5(card_path, #729 self._curr_exporter.write_procdef_mg5(card_path, #
608 self._curr_model['name'],730 self._curr_model['name'],
609731
=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py 2016-02-23 22:55:27 +0000
+++ madgraph/interface/madgraph_interface.py 2016-02-24 13:54:50 +0000
@@ -1332,7 +1332,8 @@
13321332
1333 if len(args) == 1 and args[0] in ['complex_mass_scheme',\1333 if len(args) == 1 and args[0] in ['complex_mass_scheme',\
1334 'loop_optimized_output',\1334 'loop_optimized_output',\
1335 'loop_color_flows']:1335 'loop_color_flows',\
1336 'low_mem_multicore_nlo_generation']:
1336 args.append('True')1337 args.append('True')
13371338
1338 if len(args) > 2 and '=' == args[1]:1339 if len(args) > 2 and '=' == args[1]:
@@ -1381,7 +1382,7 @@
1381 if not args[1].isdigit():1382 if not args[1].isdigit():
1382 raise self.InvalidCmd('%s values should be a integer' % args[0])1383 raise self.InvalidCmd('%s values should be a integer' % args[0])
13831384
1384 if args[0] in ['loop_optimized_output', 'loop_color_flows']:1385 if args[0] in ['loop_optimized_output', 'loop_color_flows', 'low_mem_multicore_nlo_generation']:
1385 try:1386 try:
1386 args[1] = banner_module.ConfigFile.format_variable(args[1], bool, args[0])1387 args[1] = banner_module.ConfigFile.format_variable(args[1], bool, args[0])
1387 except Exception:1388 except Exception:
@@ -2381,7 +2382,8 @@
23812382
2382 if len(args) == 2:2383 if len(args) == 2:
2383 if args[1] in ['group_subprocesses', 'complex_mass_scheme',\2384 if args[1] in ['group_subprocesses', 'complex_mass_scheme',\
2384 'loop_optimized_output', 'loop_color_flows']:2385 'loop_optimized_output', 'loop_color_flows',\
2386 'low_mem_multicore_nlo_generation']:
2385 return self.list_completion(text, ['False', 'True', 'default'])2387 return self.list_completion(text, ['False', 'True', 'default'])
2386 elif args[1] in ['ignore_six_quark_processes']:2388 elif args[1] in ['ignore_six_quark_processes']:
2387 return self.list_completion(text, self._multiparticles.keys())2389 return self.list_completion(text, self._multiparticles.keys())
@@ -2698,6 +2700,7 @@
26982700
2699 options_madgraph= {'group_subprocesses': 'Auto',2701 options_madgraph= {'group_subprocesses': 'Auto',
2700 'ignore_six_quark_processes': False,2702 'ignore_six_quark_processes': False,
2703 'low_mem_multicore_nlo_generation': False,
2701 'complex_mass_scheme': False,2704 'complex_mass_scheme': False,
2702 'gauge':'unitary',2705 'gauge':'unitary',
2703 'stdout_level':None,2706 'stdout_level':None,
27042707
=== modified file 'madgraph/iolibs/export_fks.py'
--- madgraph/iolibs/export_fks.py 2016-02-23 19:44:10 +0000
+++ madgraph/iolibs/export_fks.py 2016-02-24 13:54:50 +0000
@@ -55,6 +55,20 @@
55_file_path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0] + '/'55_file_path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0] + '/'
56logger = logging.getLogger('madgraph.export_fks')56logger = logging.getLogger('madgraph.export_fks')
5757
58
59def make_jpeg_async(args):
60 Pdir = args[0]
61 old_pos = args[1]
62 dir_path = args[2]
63
64 devnull = os.open(os.devnull, os.O_RDWR)
65
66 os.chdir(Pdir)
67 subprocess.call([os.path.join(old_pos, dir_path, 'bin', 'internal', 'gen_jpeg-pl')],
68 stdout = devnull)
69 os.chdir(os.path.pardir)
70
71
58#=================================================================================72#=================================================================================
59# Class for used of the (non-optimized) Loop process73# Class for used of the (non-optimized) Loop process
60#=================================================================================74#=================================================================================
@@ -294,12 +308,9 @@
294 #===========================================================================308 #===========================================================================
295 # write_maxparticles_file309 # write_maxparticles_file
296 #===========================================================================310 #===========================================================================
297 def write_maxparticles_file(self, writer, matrix_elements):311 def write_maxparticles_file(self, writer, maxparticles):
298 """Write the maxparticles.inc file for MadEvent"""312 """Write the maxparticles.inc file for MadEvent"""
299313
300 maxparticles = max([me.get_nexternal_ninitial()[0] \
301 for me in matrix_elements['matrix_elements']])
302
303 lines = "integer max_particles, max_branch\n"314 lines = "integer max_particles, max_branch\n"
304 lines += "parameter (max_particles=%d) \n" % maxparticles315 lines += "parameter (max_particles=%d) \n" % maxparticles
305 lines += "parameter (max_branch=max_particles-1)"316 lines += "parameter (max_branch=max_particles-1)"
@@ -313,16 +324,9 @@
313 #===========================================================================324 #===========================================================================
314 # write_maxconfigs_file325 # write_maxconfigs_file
315 #===========================================================================326 #===========================================================================
316 def write_maxconfigs_file(self, writer, matrix_elements):327 def write_maxconfigs_file(self, writer, maxconfigs):
317 """Write the maxconfigs.inc file for MadEvent"""328 """Write the maxconfigs.inc file for MadEvent"""
318329
319 try:
320 maxconfigs = max([me.get_num_configs() \
321 for me in matrix_elements['real_matrix_elements']])
322 except ValueError:
323 maxconfigs = max([me.born_matrix_element.get_num_configs() \
324 for me in matrix_elements['matrix_elements']])
325
326 lines = "integer lmaxconfigs\n"330 lines = "integer lmaxconfigs\n"
327 lines += "parameter (lmaxconfigs=%d)" % maxconfigs331 lines += "parameter (lmaxconfigs=%d)" % maxconfigs
328332
@@ -455,7 +459,7 @@
455 # likely also generate it for each subproc.459 # likely also generate it for each subproc.
456 if OLP=='NJET':460 if OLP=='NJET':
457 filename = 'OLE_order.lh'461 filename = 'OLE_order.lh'
458 self.write_lh_order(filename, matrix_element, OLP)462 self.write_lh_order(filename, [matrix_element.born_matrix_element.get('processes')[0]], OLP)
459 463
460 if matrix_element.virt_matrix_element:464 if matrix_element.virt_matrix_element:
461 calls += self.generate_virt_directory( \465 calls += self.generate_virt_directory( \
@@ -683,14 +687,11 @@
683 #===========================================================================687 #===========================================================================
684 # create the run_card 688 # create the run_card
685 #===========================================================================689 #===========================================================================
686 def create_run_card(self, matrix_elements, history):690 def create_run_card(self, processes, history):
687 """ """691 """ """
688 692
689 run_card = banner_mod.RunCardNLO()693 run_card = banner_mod.RunCardNLO()
690 694
691 processes = [me.get('processes')
692 for me in matrix_elements['matrix_elements']]
693
694 run_card.create_default_for_process(self.proc_characteristic, 695 run_card.create_default_for_process(self.proc_characteristic,
695 history,696 history,
696 processes)697 processes)
@@ -709,7 +710,7 @@
709 self.proc_characteristic['grouped_matrix'] = False710 self.proc_characteristic['grouped_matrix'] = False
710 self.create_proc_charac()711 self.create_proc_charac()
711712
712 self.create_run_card(matrix_elements, history)713 self.create_run_card(matrix_elements.get_processes(), history)
713# modelname = self.model.get('name')714# modelname = self.model.get('name')
714# if modelname == 'mssm' or modelname.startswith('mssm-'):715# if modelname == 'mssm' or modelname.startswith('mssm-'):
715# param_card = os.path.join(self.dir_path, 'Cards','param_card.dat')716# param_card = os.path.join(self.dir_path, 'Cards','param_card.dat')
@@ -723,14 +724,15 @@
723 self.write_get_mass_width_file(writers.FortranWriter(filename), makeinc, self.model)724 self.write_get_mass_width_file(writers.FortranWriter(filename), makeinc, self.model)
724725
725# # Write maxconfigs.inc based on max of ME's/subprocess groups726# # Write maxconfigs.inc based on max of ME's/subprocess groups
727
726 filename = os.path.join(self.dir_path,'Source','maxconfigs.inc')728 filename = os.path.join(self.dir_path,'Source','maxconfigs.inc')
727 self.write_maxconfigs_file(writers.FortranWriter(filename),729 self.write_maxconfigs_file(writers.FortranWriter(filename),
728 matrix_elements)730 matrix_elements.get_max_configs())
729 731
730# # Write maxparticles.inc based on max of ME's/subprocess groups732# # Write maxparticles.inc based on max of ME's/subprocess groups
731 filename = os.path.join(self.dir_path,'Source','maxparticles.inc')733 filename = os.path.join(self.dir_path,'Source','maxparticles.inc')
732 self.write_maxparticles_file(writers.FortranWriter(filename),734 self.write_maxparticles_file(writers.FortranWriter(filename),
733 matrix_elements)735 matrix_elements.get_max_particles())
734736
735 # Touch "done" file737 # Touch "done" file
736 os.system('touch %s/done' % os.path.join(self.dir_path,'SubProcesses'))738 os.system('touch %s/done' % os.path.join(self.dir_path,'SubProcesses'))
@@ -1341,7 +1343,8 @@
1341 matrix_element, i,1343 matrix_element, i,
1342 fortran_model)1344 fortran_model)
13431345
1344 def generate_virtuals_from_OLP(self,FKSHMultiproc,export_path, OLP):1346
1347 def generate_virtuals_from_OLP(self,process_list,export_path, OLP):
1345 """Generates the library for computing the loop matrix elements1348 """Generates the library for computing the loop matrix elements
1346 necessary for this process using the OLP specified."""1349 necessary for this process using the OLP specified."""
1347 1350
@@ -1350,7 +1353,7 @@
1350 if not os.path.exists(virtual_path):1353 if not os.path.exists(virtual_path):
1351 os.makedirs(virtual_path)1354 os.makedirs(virtual_path)
1352 filename = os.path.join(virtual_path,'OLE_order.lh')1355 filename = os.path.join(virtual_path,'OLE_order.lh')
1353 self.write_lh_order(filename, FKSHMultiproc.get('matrix_elements'),OLP)1356 self.write_lh_order(filename, process_list, OLP)
13541357
1355 fail_msg='Generation of the virtuals with %s failed.\n'%OLP+\1358 fail_msg='Generation of the virtuals with %s failed.\n'%OLP+\
1356 'Please check the virt_generation.log file in %s.'\1359 'Please check the virt_generation.log file in %s.'\
@@ -1419,20 +1422,19 @@
1419 proc_to_label = self.parse_contract_file(1422 proc_to_label = self.parse_contract_file(
1420 pjoin(virtual_path,'OLE_order.olc'))1423 pjoin(virtual_path,'OLE_order.olc'))
14211424
1422 self.write_BinothLHA_inc(FKSHMultiproc,proc_to_label,\1425 self.write_BinothLHA_inc(process_list,proc_to_label,\
1423 pjoin(export_path,'SubProcesses'))1426 pjoin(export_path,'SubProcesses'))
1424 1427
1425 # Link the contract file to within the SubProcess directory1428 # Link the contract file to within the SubProcess directory
1426 ln(pjoin(virtual_path,'OLE_order.olc'),pjoin(export_path,'SubProcesses'))1429 ln(pjoin(virtual_path,'OLE_order.olc'),pjoin(export_path,'SubProcesses'))
1427 1430
1428 def write_BinothLHA_inc(self, FKSHMultiproc, proc_to_label, SubProcPath):1431 def write_BinothLHA_inc(self, processes, proc_to_label, SubProcPath):
1429 """ Write the file Binoth_proc.inc in each SubProcess directory so as 1432 """ Write the file Binoth_proc.inc in each SubProcess directory so as
1430 to provide the right process_label to use in the OLP call to get the1433 to provide the right process_label to use in the OLP call to get the
1431 loop matrix element evaluation. The proc_to_label is the dictionary of1434 loop matrix element evaluation. The proc_to_label is the dictionary of
1432 the format of the one returned by the function parse_contract_file."""1435 the format of the one returned by the function parse_contract_file."""
1433 1436
1434 for matrix_element in FKSHMultiproc.get('matrix_elements'):1437 for proc in processes:
1435 proc = matrix_element.get('processes')[0]
1436 name = "P%s"%proc.shell_string()1438 name = "P%s"%proc.shell_string()
1437 proc_pdgs=(tuple([leg.get('id') for leg in proc.get('legs') if \1439 proc_pdgs=(tuple([leg.get('id') for leg in proc.get('legs') if \
1438 not leg.get('state')]),1440 not leg.get('state')]),
@@ -1608,26 +1610,19 @@
1608 # write_lh_order1610 # write_lh_order
1609 #===============================================================================1611 #===============================================================================
1610 #test written1612 #test written
1611 def write_lh_order(self, filename, matrix_elements, OLP='MadLoop'):1613 def write_lh_order(self, filename, process_list, OLP='MadLoop'):
1612 """Creates the OLE_order.lh file. This function should be edited according1614 """Creates the OLE_order.lh file. This function should be edited according
1613 to the OLP which is used. For now it is generic."""1615 to the OLP which is used. For now it is generic."""
1614 1616
1615 if isinstance(matrix_elements,fks_helas_objects.FKSHelasProcess):
1616 fksborns=fks_helas_objects.FKSHelasProcessList([matrix_elements])
1617 elif isinstance(matrix_elements,fks_helas_objects.FKSHelasProcessList):
1618 fksborns= matrix_elements
1619 else:
1620 raise fks_common.FKSProcessError('Wrong type of argument for '+\
1621 'matrix_elements in function write_lh_order.')
1622 1617
1623 if len(fksborns)==0:1618 if len(process_list)==0:
1624 raise fks_common.FKSProcessError('No matrix elements provided to '+\1619 raise fks_common.FKSProcessError('No matrix elements provided to '+\
1625 'the function write_lh_order.')1620 'the function write_lh_order.')
1626 return1621 return
1627 1622
1628 # We assume the orders to be common to all Subprocesses1623 # We assume the orders to be common to all Subprocesses
1629 1624
1630 orders = fksborns[0].orders 1625 orders = process_list[0].get('orders')
1631 if 'QED' in orders.keys() and 'QCD' in orders.keys():1626 if 'QED' in orders.keys() and 'QCD' in orders.keys():
1632 QED=orders['QED']1627 QED=orders['QED']
1633 QCD=orders['QCD']1628 QCD=orders['QCD']
@@ -1639,12 +1634,12 @@
1639 QCD=orders['QCD']1634 QCD=orders['QCD']
1640 else:1635 else:
1641 QED, QCD = self.get_qed_qcd_orders_from_weighted(\1636 QED, QCD = self.get_qed_qcd_orders_from_weighted(\
1642 fksborns[0].get_nexternal_ninitial()[0]-1, # -1 is because the function returns nexternal of the real emission1637 len(process_list[0].get('legs')),
1643 orders['WEIGHTED'])1638 orders['WEIGHTED'])
16441639
1645 replace_dict = {}1640 replace_dict = {}
1646 replace_dict['mesq'] = 'CHaveraged'1641 replace_dict['mesq'] = 'CHaveraged'
1647 replace_dict['corr'] = ' '.join(matrix_elements[0].get('processes')[0].\1642 replace_dict['corr'] = ' '.join(process_list[0].\
1648 get('perturbation_couplings'))1643 get('perturbation_couplings'))
1649 replace_dict['irreg'] = 'CDR'1644 replace_dict['irreg'] = 'CDR'
1650 replace_dict['aspow'] = QCD1645 replace_dict['aspow'] = QCD
@@ -1652,8 +1647,10 @@
1652 replace_dict['modelfile'] = './param_card.dat'1647 replace_dict['modelfile'] = './param_card.dat'
1653 replace_dict['params'] = 'alpha_s'1648 replace_dict['params'] = 'alpha_s'
1654 proc_lines=[]1649 proc_lines=[]
1655 for fksborn in fksborns:1650 for proc in process_list:
1656 proc_lines.append(fksborn.get_lh_pdg_string())1651 proc_lines.append('%s -> %s' % \
1652 (' '.join(str(l['id']) for l in proc['legs'] if not l['state']),
1653 ' '.join(str(l['id']) for l in proc['legs'] if l['state'])))
1657 replace_dict['pdgs'] = '\n'.join(proc_lines)1654 replace_dict['pdgs'] = '\n'.join(proc_lines)
1658 replace_dict['symfin'] = 'Yes'1655 replace_dict['symfin'] = 'Yes'
1659 content = \1656 content = \
@@ -3300,7 +3297,7 @@
3300 #===============================================================================3297 #===============================================================================
3301 # write_coef_specs3298 # write_coef_specs
3302 #===============================================================================3299 #===============================================================================
3303 def write_coef_specs_file(self, virt_me_list):3300 def write_coef_specs_file(self, max_loop_vertex_ranks):
3304 """ writes the coef_specs.inc in the DHELAS folder. Should not be called in the 3301 """ writes the coef_specs.inc in the DHELAS folder. Should not be called in the
3305 non-optimized mode"""3302 non-optimized mode"""
3306 filename = os.path.join(self.dir_path, 'Source', 'DHELAS', 'coef_specs.inc')3303 filename = os.path.join(self.dir_path, 'Source', 'DHELAS', 'coef_specs.inc')
@@ -3308,7 +3305,6 @@
3308 replace_dict = {}3305 replace_dict = {}
3309 replace_dict['max_lwf_size'] = 4 3306 replace_dict['max_lwf_size'] = 4
33103307
3311 max_loop_vertex_ranks = [me.get_max_loop_vertex_rank() for me in virt_me_list]
3312 replace_dict['vertex_max_coefs'] = max(\3308 replace_dict['vertex_max_coefs'] = max(\
3313 [q_polynomial.get_number_of_coefs_for_rank(n) 3309 [q_polynomial.get_number_of_coefs_for_rank(n)
3314 for n in max_loop_vertex_ranks])3310 for n in max_loop_vertex_ranks])
33153311
=== modified file 'madgraph/iolibs/template_files/loop_optimized/compute_color_flows.inc'
--- madgraph/iolibs/template_files/loop_optimized/compute_color_flows.inc 2015-07-03 00:47:03 +0000
+++ madgraph/iolibs/template_files/loop_optimized/compute_color_flows.inc 2016-02-24 13:54:50 +0000
@@ -357,7 +357,7 @@
357C 357C
358C LOCAL VARIABLES 358C LOCAL VARIABLES
359C359C
360 INTEGER I, J, M, N, K, SQSOINDEX, DOUBLEFACT360 INTEGER I, J, M, N, K, ISQSO, DOUBLEFACT
361 %(complex_dp_format)s COLOR_COEF361 %(complex_dp_format)s COLOR_COEF
362C 362C
363C FUNCTIONS363C FUNCTIONS
@@ -401,8 +401,8 @@
401 IF (LoopColorFlowMatrix(I,I)%%Denom.LT.0) COLOR_COEF=COLOR_COEF*IMAG1401 IF (LoopColorFlowMatrix(I,I)%%Denom.LT.0) COLOR_COEF=COLOR_COEF*IMAG1
402 DO M=1,NLOOPAMPSO402 DO M=1,NLOOPAMPSO
403 DO N=M,NLOOPAMPSO403 DO N=M,NLOOPAMPSO
404 SQSOINDEX = %(proc_prefix)sML5SQSOINDEX(M,N)404 ISQSO = %(proc_prefix)sML5SQSOINDEX(M,N)
405 IF(CHOSEN_SO_CONFIGS(SQSOINDEX)) THEN405 IF(CHOSEN_SO_CONFIGS(ISQSO)) THEN
406 IF(M.ne.N) THEN406 IF(M.ne.N) THEN
407 DOUBLEFACT=2 407 DOUBLEFACT=2
408 ELSE 408 ELSE
@@ -452,7 +452,7 @@
452C 452C
453C LOCAL VARIABLES 453C LOCAL VARIABLES
454C454C
455 INTEGER I, J, M, N, K, SQSOINDEX, DOUBLEFACT455 INTEGER I, J, M, N, K, ISQSO, DOUBLEFACT
456## if(LoopInduced) {456## if(LoopInduced) {
457 INTEGER LOWERBOUND457 INTEGER LOWERBOUND
458## }458## }
@@ -550,15 +550,15 @@
550 CALL %(proc_prefix)sML5GET_ORDERS_FOR_AMPSOINDEX(M,ORDERS_A)550 CALL %(proc_prefix)sML5GET_ORDERS_FOR_AMPSOINDEX(M,ORDERS_A)
551 CALL %(proc_prefix)sML5GET_ORDERS_FOR_AMPSOINDEX(N,ORDERS_B)551 CALL %(proc_prefix)sML5GET_ORDERS_FOR_AMPSOINDEX(N,ORDERS_B)
552C Now figure out to which SQSOINDEX these orders together correspond to *in the Born subroutine* (Notice the absence of the ML5 prefix in the functions used then)552C Now figure out to which SQSOINDEX these orders together correspond to *in the Born subroutine* (Notice the absence of the ML5 prefix in the functions used then)
553 SQSOINDEX = %(proc_prefix)sSQSOINDEX(%(proc_prefix)sSOINDEX_FOR_AMPORDERS(ORDERS_A),%(proc_prefix)sSOINDEX_FOR_AMPORDERS(ORDERS_B))553 ISQSO = %(proc_prefix)sSQSOINDEX(%(proc_prefix)sSOINDEX_FOR_AMPORDERS(ORDERS_A),%(proc_prefix)sSOINDEX_FOR_AMPORDERS(ORDERS_B))
554 IF(J.ne.I) THEN554 IF(J.ne.I) THEN
555 DOUBLEFACT=2 555 DOUBLEFACT=2
556 ELSE 556 ELSE
557 DOUBLEFACT=1557 DOUBLEFACT=1
558 ENDIF558 ENDIF
559 TEMP(1) = DOUBLEFACT*HEL_MULT*DBLE(COLOR_COEF*JAMPB(I,M)*DCONJG(JAMPB(J,N)))559 TEMP(1) = DOUBLEFACT*HEL_MULT*DBLE(COLOR_COEF*JAMPB(I,M)*DCONJG(JAMPB(J,N)))
560 RES(0,SQSOINDEX) = RES(0,SQSOINDEX) + TEMP(1)560 RES(0,ISQSO) = RES(0,ISQSO) + TEMP(1)
561 IF((.not.FILTER_SO).or.SQSO_TARGET.eq.-1.or.SQSO_TARGET.eq.SQSOINDEX) THEN561 IF((.not.FILTER_SO).or.SQSO_TARGET.eq.-1.or.SQSO_TARGET.eq.ISQSO) THEN
562 RES(0,0) = RES(0,0) + TEMP(1)562 RES(0,0) = RES(0,0) + TEMP(1)
563 ENDIF563 ENDIF
564 ENDDO564 ENDDO
@@ -591,7 +591,7 @@
591 ENDIF591 ENDIF
592 DO N=LOWERBOUND,NLOOPAMPSO592 DO N=LOWERBOUND,NLOOPAMPSO
593## }593## }
594 SQSOINDEX = %(proc_prefix)sML5SQSOINDEX(M,N)594 ISQSO = %(proc_prefix)sML5SQSOINDEX(M,N)
595## if(LoopInduced){595## if(LoopInduced){
596 IF(J.ne.I) THEN596 IF(J.ne.I) THEN
597 DOUBLEFACT=2 597 DOUBLEFACT=2
@@ -626,14 +626,14 @@
626## if(not LoopInduced){ 626## if(not LoopInduced){
627 TEMP(K) = 2.0d0*HEL_MULT*DBLE(COLOR_COEF*JAMPL(K,I,M)*DCONJG(JAMPB(J,N)))627 TEMP(K) = 2.0d0*HEL_MULT*DBLE(COLOR_COEF*JAMPL(K,I,M)*DCONJG(JAMPB(J,N)))
628## }628## }
629 RES(K,SQSOINDEX) = RES(K,SQSOINDEX) + TEMP(K)629 RES(K,ISQSO) = RES(K,ISQSO) + TEMP(K)
630## if(LoopInduced and MadEventOutput){630## if(LoopInduced and MadEventOutput){
631 DO config_i=1,nconfigs631 DO config_i=1,nconfigs
632 AMP2_ALL(K,config_i,SQSOINDEX) = AMP2_ALL(K,config_i,SQSOINDEX) + TEMP_AMP2(K,config_i)632 AMP2_ALL(K,config_i,ISQSO) = AMP2_ALL(K,config_i,ISQSO) + TEMP_AMP2(K,config_i)
633 ENDDO633 ENDDO
634## }634## }
635 ENDDO635 ENDDO
636 IF((.NOT.FILTER_SO).OR.SQSO_TARGET.eq.-1.or.SQSO_TARGET.eq.SQSOINDEX) THEN636 IF((.NOT.FILTER_SO).OR.SQSO_TARGET.eq.-1.or.SQSO_TARGET.eq.ISQSO) THEN
637 DO K=1,3637 DO K=1,3
638 RES(K,0) = RES(K,0) + TEMP(K)638 RES(K,0) = RES(K,0) + TEMP(K)
639## if(LoopInduced and MadEventOutput){639## if(LoopInduced and MadEventOutput){
640640
=== modified file 'madgraph/loop/loop_exporters.py'
--- madgraph/loop/loop_exporters.py 2016-02-23 22:55:27 +0000
+++ madgraph/loop/loop_exporters.py 2016-02-24 13:54:50 +0000
@@ -45,6 +45,7 @@
45import madgraph.various.diagram_symmetry as diagram_symmetry45import madgraph.various.diagram_symmetry as diagram_symmetry
46import madgraph.various.process_checks as process_checks46import madgraph.various.process_checks as process_checks
47import madgraph.various.progressbar as pbar47import madgraph.various.progressbar as pbar
48import madgraph.various.q_polynomial as q_polynomial
48import madgraph.core.color_amp as color_amp49import madgraph.core.color_amp as color_amp
49import madgraph.iolibs.helas_call_writers as helas_call_writers50import madgraph.iolibs.helas_call_writers as helas_call_writers
50import models.check_param_card as check_param_card51import models.check_param_card as check_param_card
@@ -2314,7 +2315,7 @@
23142315
2315 writer.writelines(file,context=self.get_context(matrix_element))2316 writer.writelines(file,context=self.get_context(matrix_element))
2316 2317
2317 def fix_coef_specs(self, overall_max_lwf_size, overall_max_loop_vert_rank):2318 def fix_coef_specs(self, overall_max_lwf_spin, overall_max_loop_vert_rank):
2318 """ If processes with different maximum loop wavefunction size or2319 """ If processes with different maximum loop wavefunction size or
2319 different maximum loop vertex rank have to be output together, then2320 different maximum loop vertex rank have to be output together, then
2320 the file 'coef.inc' in the HELAS Source folder must contain the overall2321 the file 'coef.inc' in the HELAS Source folder must contain the overall
@@ -2325,7 +2326,11 @@
2325 coef_specs_path=os.path.join(self.dir_path,'Source','DHELAS',\2326 coef_specs_path=os.path.join(self.dir_path,'Source','DHELAS',\
2326 'coef_specs.inc')2327 'coef_specs.inc')
2327 os.remove(coef_specs_path)2328 os.remove(coef_specs_path)
2328 2329
2330 spin_to_wf_size = {1:4,2:4,3:4,4:16,5:16}
2331 overall_max_lwf_size = spin_to_wf_size[overall_max_lwf_spin]
2332 overall_max_loop_vert_coefs = q_polynomial.get_number_of_coefs_for_rank(
2333 overall_max_loop_vert_rank)
2329 # Replace it by the appropriate value2334 # Replace it by the appropriate value
2330 IncWriter=writers.FortranWriter(coef_specs_path,'w')2335 IncWriter=writers.FortranWriter(coef_specs_path,'w')
2331 IncWriter.writelines("""INTEGER MAXLWFSIZE2336 IncWriter.writelines("""INTEGER MAXLWFSIZE
@@ -2333,7 +2338,7 @@
2333 INTEGER VERTEXMAXCOEFS2338 INTEGER VERTEXMAXCOEFS
2334 PARAMETER (VERTEXMAXCOEFS=%(vertex_max_coefs)d)"""\2339 PARAMETER (VERTEXMAXCOEFS=%(vertex_max_coefs)d)"""\
2335 %{'max_lwf_size':overall_max_lwf_size,2340 %{'max_lwf_size':overall_max_lwf_size,
2336 'vertex_max_coefs':overall_max_loop_vert_rank})2341 'vertex_max_coefs':overall_max_loop_vert_coefs})
2337 IncWriter.close()2342 IncWriter.close()
23382343
2339 def setup_check_sa_replacement_dictionary(self, matrix_element, \2344 def setup_check_sa_replacement_dictionary(self, matrix_element, \
23402345
=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py 2016-02-18 11:42:56 +0000
+++ madgraph/various/banner.py 2016-02-24 13:54:50 +0000
@@ -1990,7 +1990,7 @@
1990 # check for beam_id1990 # check for beam_id
1991 beam_id = set()1991 beam_id = set()
1992 for proc in proc_def:1992 for proc in proc_def:
1993 for leg in proc[0]['legs']:1993 for leg in proc['legs']:
1994 if not leg['state']:1994 if not leg['state']:
1995 beam_id.add(leg['id'])1995 beam_id.add(leg['id'])
1996 if any(i in beam_id for i in [1,-1,2,-2,3,-3,4,-4,5,-5,21,22]):1996 if any(i in beam_id for i in [1,-1,2,-2,3,-3,4,-4,5,-5,21,22]):
19971997
=== added directory 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall'
=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f 1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f 2016-02-24 13:54:50 +0000
@@ -0,0 +1,278 @@
1 SUBROUTINE CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
2C
3C Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4C By the MadGraph5_aMC@NLO Development Team
5C Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
6C
7C Interface between MG5 and CutTools.
8C
9C Process: d~ u > w+ WEIGHTED=2 QED=1 [ all = QCD ]
10C Process: s~ c > w+ WEIGHTED=2 QED=1 [ all = QCD ]
11C
12C
13C CONSTANTS
14C
15 INTEGER NEXTERNAL
16 PARAMETER (NEXTERNAL=3)
17 LOGICAL CHECKPCONSERVATION
18 PARAMETER (CHECKPCONSERVATION=.TRUE.)
19 REAL*8 NORMALIZATION
20 PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
21 $ *2))
22C
23C ARGUMENTS
24C
25 INTEGER NLOOPLINE, RANK
26 REAL*8 PL(0:3,NLOOPLINE)
27 REAL*8 PCT(0:3,0:NLOOPLINE-1)
28 COMPLEX*16 M2L(NLOOPLINE)
29 COMPLEX*16 M2LCT(0:NLOOPLINE-1)
30 COMPLEX*16 RES(3)
31 LOGICAL STABLE
32C
33C LOCAL VARIABLES
34C
35 COMPLEX*16 R1, ACC
36 INTEGER I, J, K
37 LOGICAL CTINIT, TIRINIT, GOLEMINIT
38 COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT
39C
40C EXTERNAL FUNCTIONS
41C
42 EXTERNAL LOOPNUM
43 EXTERNAL MPLOOPNUM
44C
45C GLOBAL VARIABLES
46C
47 INCLUDE 'coupl.inc'
48 INTEGER CTMODE
49 REAL*8 LSCALE
50 COMMON/CT/LSCALE,CTMODE
51
52 INTEGER ID,SQSOINDEX,R
53 COMMON/LOOP/ID,SQSOINDEX,R
54
55C ----------
56C BEGIN CODE
57C ----------
58
59C INITIALIZE CUTTOOLS IF NEEDED
60 IF (CTINIT) THEN
61 CTINIT=.FALSE.
62 CALL INITCT()
63 ENDIF
64
65C YOU CAN FIND THE DETAILS ABOUT THE DIFFERENT CTMODE AT THE
66C BEGINNING OF THE FILE CTS_CUTS.F90 IN THE CUTTOOLS DISTRIBUTION
67
68C CONVERT THE MASSES TO BE COMPLEX
69 DO I=1,NLOOPLINE
70 M2LCT(I-1)=M2L(I)
71 ENDDO
72
73C CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
74 DO I=0,3
75 DO J=0,(NLOOPLINE-1)
76 PCT(I,J)=0.D0
77 ENDDO
78 ENDDO
79 DO I=0,3
80 DO J=1,NLOOPLINE
81 PCT(I,0)=PCT(I,0)+PL(I,J)
82 ENDDO
83 ENDDO
84 IF (CHECKPCONSERVATION) THEN
85 IF (PCT(0,0).GT.1.D-6) THEN
86 WRITE(*,*) 'energy is not conserved ',PCT(0,0)
87 STOP 'energy is not conserved'
88 ELSEIF (PCT(1,0).GT.1.D-6) THEN
89 WRITE(*,*) 'px is not conserved ',PCT(1,0)
90 STOP 'px is not conserved'
91 ELSEIF (PCT(2,0).GT.1.D-6) THEN
92 WRITE(*,*) 'py is not conserved ',PCT(2,0)
93 STOP 'py is not conserved'
94 ELSEIF (PCT(3,0).GT.1.D-6) THEN
95 WRITE(*,*) 'pz is not conserved ',PCT(3,0)
96 STOP 'pz is not conserved'
97 ENDIF
98 ENDIF
99 DO I=0,3
100 DO J=1,(NLOOPLINE-1)
101 DO K=1,J
102 PCT(I,J)=PCT(I,J)+PL(I,K)
103 ENDDO
104 ENDDO
105 ENDDO
106
107 CALL CTSXCUT(CTMODE,LSCALE,MU_R,NLOOPLINE,LOOPNUM,MPLOOPNUM,RANK
108 $ ,PCT,M2LCT,RES,ACC,R1,STABLE)
109 RES(1)=NORMALIZATION*2.0D0*DBLE(RES(1))
110 RES(2)=NORMALIZATION*2.0D0*DBLE(RES(2))
111 RES(3)=NORMALIZATION*2.0D0*DBLE(RES(3))
112C WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
113 END
114
115 SUBROUTINE INITCT()
116C
117C INITIALISATION OF CUTTOOLS
118C
119C LOCAL VARIABLES
120C
121 REAL*8 THRS
122 LOGICAL EXT_NUM_FOR_R1
123C
124C GLOBAL VARIABLES
125C
126 INCLUDE 'MadLoopParams.inc'
127C ----------
128C BEGIN CODE
129C ----------
130
131C DEFAULT PARAMETERS FOR CUTTOOLS
132C -------------------------------
133C THRS1 IS THE PRECISION LIMIT BELOW WHICH THE MP ROUTINES
134C ACTIVATES
135 THRS=CTSTABTHRES
136C LOOPLIB SET WHAT LIBRARY CT USES
137C 1 -> LOOPTOOLS
138C 2 -> AVH
139C 3 -> QCDLOOP
140 LOOPLIB=CTLOOPLIBRARY
141C MADLOOP'S NUMERATOR IN THE OPEN LOOP IS MUCH FASTER THAN THE
142C RECONSTRUCTED ONE IN CT. SO WE BETTER USE MADLOOP ONE IN THIS
143C CASE.
144 EXT_NUM_FOR_R1=.TRUE.
145C -------------------------------
146
147C The initialization below is for CT v1.8.+
148 CALL CTSINIT(THRS,LOOPLIB,EXT_NUM_FOR_R1)
149C The initialization below is for the older stable CT v1.7, still
150C used for now in the beta release.
151C CALL CTSINIT(THRS,LOOPLIB)
152
153 END
154
155 SUBROUTINE LOOP_3(W1, W2, W3, M1, M2, M3, RANK, SQUAREDSOINDEX
156 $ , LOOPNUM)
157 INTEGER NEXTERNAL
158 PARAMETER (NEXTERNAL=3)
159 INTEGER NLOOPLINE
160 PARAMETER (NLOOPLINE=3)
161 INTEGER NWAVEFUNCS
162 PARAMETER (NWAVEFUNCS=3)
163 INTEGER NLOOPGROUPS
164 PARAMETER (NLOOPGROUPS=1)
165 INTEGER NCOMB
166 PARAMETER (NCOMB=12)
167C These are constants related to the split orders
168 INTEGER NSQUAREDSO
169 PARAMETER (NSQUAREDSO=1)
170C
171C ARGUMENTS
172C
173 INTEGER W1, W2, W3
174 COMPLEX*16 M1, M2, M3
175
176 INTEGER RANK, LSYMFACT
177 INTEGER LOOPNUM, SQUAREDSOINDEX
178C
179C LOCAL VARIABLES
180C
181 REAL*8 PL(0:3,NLOOPLINE)
182 COMPLEX*16 M2L(NLOOPLINE)
183 INTEGER PAIRING(NLOOPLINE),WE(3)
184 INTEGER I, J, K, TEMP,I_LIB
185 LOGICAL COMPLEX_MASS,DOING_QP
186C
187C GLOBAL VARIABLES
188C
189 INCLUDE 'MadLoopParams.inc'
190 INTEGER ID,SQSOINDEX,R
191 COMMON/LOOP/ID,SQSOINDEX,R
192
193 LOGICAL CHECKPHASE, HELDOUBLECHECKED
194 COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED
195
196 INTEGER HELOFFSET
197 INTEGER GOODHEL(NCOMB)
198 LOGICAL GOODAMP(NSQUAREDSO,NLOOPGROUPS)
199 COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET
200
201 COMPLEX*16 LOOPRES(3,NSQUAREDSO,NLOOPGROUPS)
202 LOGICAL S(NSQUAREDSO,NLOOPGROUPS)
203 COMMON/LOOPRES/LOOPRES,S
204
205
206 COMPLEX*16 W(20,NWAVEFUNCS)
207 COMMON/W/W
208 REAL*8 LSCALE
209 INTEGER CTMODE
210 COMMON/CT/LSCALE,CTMODE
211 INTEGER LIBINDEX
212 COMMON/I_LIB/LIBINDEX
213
214C ----------
215C BEGIN CODE
216C ----------
217
218 IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(SQUAREDSOIND
219 $ EX,LOOPNUM)) THEN
220 WE(1)=W1
221 WE(2)=W2
222 WE(3)=W3
223 M2L(1)=M3**2
224 M2L(2)=M1**2
225 M2L(3)=M2**2
226 DO I=1,NLOOPLINE
227 PAIRING(I)=1
228 ENDDO
229
230 R=RANK
231 ID=LOOPNUM
232 SQSOINDEX=SQUAREDSOINDEX
233 DO I=0,3
234 TEMP=1
235 DO J=1,NLOOPLINE
236 PL(I,J)=0.D0
237 DO K=TEMP,(TEMP+PAIRING(J)-1)
238 PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
239 ENDDO
240 TEMP=TEMP+PAIRING(J)
241 ENDDO
242 ENDDO
243C Determine whether the integral is with complex masses or not
244C since some reduction libraries, e.g.PJFry++ and IREGI, are
245C still
246C not able to deal with complex masses
247 COMPLEX_MASS=.FALSE.
248 DO I=1,NLOOPLINE
249 IF(DIMAG(M2L(I)).EQ.0D0)CYCLE
250 IF(ABS(DIMAG(M2L(I)))/MAX(ABS(M2L(I)),1D-2).GT.1D-15)THEN
251 COMPLEX_MASS=.TRUE.
252 EXIT
253 ENDIF
254 ENDDO
255C Determine it uses qp or not
256 DOING_QP=.FALSE.
257 IF(CTMODE.GE.4)DOING_QP=.TRUE.
258C Choose the correct loop library
259 CALL CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
260 $ ,DOING_QP,I_LIB)
261 IF(MLREDUCTIONLIB(I_LIB).EQ.1)THEN
262C USE CUTTOOLS
263 CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,LOOPRES(1,SQUAREDSOINDEX
264 $ ,LOOPNUM),S(SQUAREDSOINDEX,LOOPNUM))
265 ELSE
266C USE TIR
267 CALL TIRLOOP(SQUAREDSOINDEX,LOOPNUM,I_LIB,NLOOPLINE,PL,M2L
268 $ ,RANK,LOOPRES(1,SQUAREDSOINDEX,LOOPNUM),S(SQUAREDSOINDEX
269 $ ,LOOPNUM))
270 ENDIF
271 ELSE
272 LOOPRES(1,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
273 LOOPRES(2,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
274 LOOPRES(3,SQUAREDSOINDEX,LOOPNUM)=(0.0D0,0.0D0)
275 S(SQUAREDSOINDEX,LOOPNUM)=.TRUE.
276 ENDIF
277 END
278
0279
=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f 1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f 2016-02-24 13:54:50 +0000
@@ -0,0 +1,760 @@
1 SUBROUTINE GOLEMLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
2C
3C Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
4C By the MadGraph5_aMC@NLO Development Team
5C Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
6C
7C Interface between MG5 and Golem95.
8C The Golem95 version should be higher than 1.3.0.
9C It supports RANK = NLOOPLINE + 1 tensor integrals when 1 <
10C NLOOPLINE < 6.
11C
12C Process: d~ u > w+ WEIGHTED=2 QED=1 [ all = QCD ]
13C Process: s~ c > w+ WEIGHTED=2 QED=1 [ all = QCD ]
14C
15C
16C MODULES
17C
18 USE MATRICE_S
19 USE FORM_FACTOR_TYPE, ONLY: FORM_FACTOR
20 USE PRECISION_GOLEM, ONLY: KI
21 USE TENS_COMB
22 USE TENS_REC
23 USE FORM_FACTOR_1P, ONLY: A10
24 USE FORM_FACTOR_2P, ONLY: A20
25 USE FORM_FACTOR_3P, ONLY: A30
26 USE FORM_FACTOR_4P, ONLY: A40
27 USE FORM_FACTOR_5P, ONLY: A50
28 USE FORM_FACTOR_6P, ONLY: A60
29 IMPLICIT NONE
30C
31C CONSTANTS
32C
33 INTEGER NEXTERNAL
34 PARAMETER (NEXTERNAL=3)
35 LOGICAL CHECKPCONSERVATION
36 PARAMETER (CHECKPCONSERVATION=.TRUE.)
37 REAL*8 NORMALIZATION
38 PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
39 $ *2))
40 REAL(KI),DIMENSION(0:3),PARAMETER::NULL_VEC = (/0.0_KI,0.0_KI
41 $ ,0.0_KI,0.0_KI/)
42C GOLEM_RUN_MODE = 1: Use directly MadLoop tensorial coefficients
43C GOLEM_RUN_MODE = 2: Reconstruct the tensorial coefficeints
44C directly from
45C numerator using golem internal reconstruction routine
46C GOLEM_RUN_MODE = 3: Cross-checked reconstructed coefficients
47C against
48C MadLoop internal ones.
49 INTEGER GOLEM_RUN_MODE
50 PARAMETER (GOLEM_RUN_MODE=1)
51C The following is the acceptance threshold used for GOLEM_RUN_MODE
52C = 3
53 REAL*8 COEF_CHECK_THRS
54 DATA COEF_CHECK_THRS/1.0D-13/
55 COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
56
57 LOGICAL PASS_COEF_CHECK
58C
59C ARGUMENTS
60C
61 INTEGER NLOOPLINE, RANK
62 REAL*8 PL(0:3,NLOOPLINE)
63 REAL*8 PCT(0:3,0:NLOOPLINE-1)
64 REAL(KI) PGOLEM(NLOOPLINE,0:3)
65 COMPLEX*16 M2L(NLOOPLINE)
66 COMPLEX(KI) M2LGOLEM(NLOOPLINE)
67 COMPLEX*16 RES(3)
68 LOGICAL STABLE
69C
70C LOCAL VARIABLES
71C
72 INTEGER I, J, K
73 TYPE(FORM_FACTOR)::RES_GOLEM
74
75 COMPLEX(KI)::COEFFS0,COEFFS0_REC
76 TYPE(COEFF_TYPE_1)::COEFFS1,COEFFS1_REC
77 TYPE(COEFF_TYPE_2)::COEFFS2,COEFFS2_REC
78 TYPE(COEFF_TYPE_3)::COEFFS3,COEFFS3_REC
79 TYPE(COEFF_TYPE_4)::COEFFS4,COEFFS4_REC
80 TYPE(COEFF_TYPE_5)::COEFFS5,COEFFS5_REC
81 TYPE(COEFF_TYPE_6)::COEFFS6,COEFFS6_REC
82
83C The pinch propagator optimization is not used, so for now it is
84C always 0.
85 INTEGER PINCH
86C
87C EXTERNAL FUNCTIONS
88C
89 COMPLEX(KI) GOLEM_LOOPNUM
90 EXTERNAL GOLEM_LOOPNUM
91 LOGICAL COMPARE_COEFS_0
92 LOGICAL COMPARE_COEFS_1
93 LOGICAL COMPARE_COEFS_2
94 LOGICAL COMPARE_COEFS_3
95 LOGICAL COMPARE_COEFS_4
96 LOGICAL COMPARE_COEFS_5
97 LOGICAL COMPARE_COEFS_6
98C
99C GLOBAL VARIABLES
100C
101 INCLUDE 'coupl.inc'
102 INTEGER CTMODE
103 REAL*8 LSCALE
104 COMMON/CT/LSCALE,CTMODE
105
106 INTEGER ID,SQSOINDEX,R
107 COMMON/LOOP/ID,SQSOINDEX,R
108
109 LOGICAL CTINIT, TIRINIT, GOLEMINIT
110 COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT
111
112 INTEGER NLOOPGROUPS
113 PARAMETER (NLOOPGROUPS=1)
114 INTEGER LOOPMAXCOEFS
115 PARAMETER (LOOPMAXCOEFS=15)
116 INTEGER NSQUAREDSO
117 PARAMETER (NSQUAREDSO=1)
118
119 COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
120 COMMON/LCOEFS/LOOPCOEFS
121C ----------
122C BEGIN CODE
123C ----------
124
125C INITIALIZE CUTTOOLS IF NEEDED
126 IF (GOLEMINIT) THEN
127 GOLEMINIT=.FALSE.
128 CALL INITGOLEM()
129 ENDIF
130
131C No stability test intrisic to Golem95 now
132 STABLE=.TRUE.
133
134C This initialization must be done for each reduction because we
135C have not setup anyoptimization using pinched propagators yet.
136 CALL INITGOLEM95(NLOOPLINE)
137 PINCH = 0
138
139C YOU CAN FIND THE DETAILS ABOUT THE DIFFERENT CTMODE AT THE
140C BEGINNING OF THE FILE CTS_CUTS.F90 IN THE CUTTOOLS DISTRIBUTION
141
142C CONVERT THE MASSES TO BE COMPLEX
143 DO I=1,NLOOPLINE
144 M2LGOLEM(I)=M2L(I)
145 ENDDO
146
147C CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
148 DO I=0,3
149 DO J=0,(NLOOPLINE-1)
150 PCT(I,J)=0.D0
151 ENDDO
152 ENDDO
153 DO I=0,3
154 DO J=1,NLOOPLINE
155 PCT(I,0)=PCT(I,0)+PL(I,J)
156 ENDDO
157 ENDDO
158 IF (CHECKPCONSERVATION) THEN
159 IF (PCT(0,0).GT.1.D-6) THEN
160 WRITE(*,*) 'energy is not conserved ',PCT(0,0)
161 STOP 'energy is not conserved'
162 ELSEIF (PCT(1,0).GT.1.D-6) THEN
163 WRITE(*,*) 'px is not conserved ',PCT(1,0)
164 STOP 'px is not conserved'
165 ELSEIF (PCT(2,0).GT.1.D-6) THEN
166 WRITE(*,*) 'py is not conserved ',PCT(2,0)
167 STOP 'py is not conserved'
168 ELSEIF (PCT(3,0).GT.1.D-6) THEN
169 WRITE(*,*) 'pz is not conserved ',PCT(3,0)
170 STOP 'pz is not conserved'
171 ENDIF
172 ENDIF
173 DO I=0,3
174 DO J=1,(NLOOPLINE-1)
175 DO K=1,J
176 PCT(I,J)=PCT(I,J)+PL(I,K)
177 ENDDO
178 ENDDO
179 ENDDO
180
181C Now convert the loop momenta to Golem95 conventions
182 DO I=0,3
183 PGOLEM(1,I)=0.0E0_KI
184 DO J=2,NLOOPLINE
185 PGOLEM(J,I)=PCT(I,J-1)
186 ENDDO
187 ENDDO
188
189C Fill in the kinematic s-matrix while taking care of on-shell
190C limits.
191 CALL SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2LGOLEM)
192C Construct the golem internal matrices derived from the kinetic
193C one.
194 CALL PREPARESMATRIX()
195
196C Fill in the golem coefficents and compute the loop
197 IF(GOLEM_RUN_MODE.EQ.2)THEN
198 RES_GOLEM = EVALUATE_B(GOLEM_LOOPNUM,PGOLEM,0,RANK)
199 ELSE
200 PASS_COEF_CHECK=.TRUE.
201 SELECT CASE(RANK)
202 CASE(0)
203 CALL FILL_GOLEM_COEFFS_0(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS0)
204 IF(GOLEM_RUN_MODE.EQ.3)THEN
205 COEFFS0_REC = GOLEM_LOOPNUM(NULL_VEC,0.0_KI)
206 PASS_COEF_CHECK=COMPARE_COEFS_0(COEFFS0,COEFFS0_REC)
207 ENDIF
208 CASE(1)
209 CALL FILL_GOLEM_COEFFS_1(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS1)
210 IF(GOLEM_RUN_MODE.EQ.3)THEN
211 CALL RECONSTRUCT1(GOLEM_LOOPNUM,COEFFS1_REC)
212 PASS_COEF_CHECK=COMPARE_COEFS_1(COEFFS1,COEFFS1_REC)
213 ENDIF
214 CASE(2)
215 CALL FILL_GOLEM_COEFFS_2(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS2)
216 IF(GOLEM_RUN_MODE.EQ.3)THEN
217 CALL RECONSTRUCT2(GOLEM_LOOPNUM,COEFFS2_REC)
218 PASS_COEF_CHECK=COMPARE_COEFS_2(COEFFS2,COEFFS2_REC)
219 ENDIF
220 CASE(3)
221 CALL FILL_GOLEM_COEFFS_3(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS3)
222 IF(GOLEM_RUN_MODE.EQ.3)THEN
223 CALL RECONSTRUCT3(GOLEM_LOOPNUM,COEFFS3_REC)
224 PASS_COEF_CHECK=COMPARE_COEFS_3(COEFFS3,COEFFS3_REC)
225 ENDIF
226 CASE(4)
227 CALL FILL_GOLEM_COEFFS_4(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS4)
228 IF(GOLEM_RUN_MODE.EQ.3)THEN
229 CALL RECONSTRUCT4(GOLEM_LOOPNUM,COEFFS4_REC)
230 PASS_COEF_CHECK=COMPARE_COEFS_4(COEFFS4,COEFFS4_REC)
231 ENDIF
232 CASE(5)
233 CALL FILL_GOLEM_COEFFS_5(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS5)
234 IF(GOLEM_RUN_MODE.EQ.3)THEN
235 CALL RECONSTRUCT5(GOLEM_LOOPNUM,COEFFS5_REC)
236 PASS_COEF_CHECK=COMPARE_COEFS_5(COEFFS5,COEFFS5_REC)
237 ENDIF
238 CASE(6)
239 CALL FILL_GOLEM_COEFFS_6(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS6)
240 IF(GOLEM_RUN_MODE.EQ.3)THEN
241 CALL RECONSTRUCT6(GOLEM_LOOPNUM,COEFFS6_REC)
242 PASS_COEF_CHECK=COMPARE_COEFS_6(COEFFS6,COEFFS6_REC)
243 ENDIF
244 CASE DEFAULT
245 WRITE(*,*)'Not yet implemented in Golem95 for rank= ',RANK
246 STOP
247 END SELECT
248
249 IF(.NOT.PASS_COEF_CHECK)THEN
250 WRITE(*,*)'Coefs mismatch for ID ',ID,' and rank ',RANK
251 WRITE(*,*)'Coefs form MadLoop5:'
252 SELECT CASE(RANK)
253 CASE(0)
254 WRITE(*,*)'Constant coef = ',COEFFS0
255 CASE(1)
256 CALL PRINT_COEFFS(COEFFS1)
257 CASE(2)
258 CALL PRINT_COEFFS(COEFFS2)
259 CASE(3)
260 CALL PRINT_COEFFS(COEFFS3)
261 CASE(4)
262 CALL PRINT_COEFFS(COEFFS4)
263 CASE(5)
264 CALL PRINT_COEFFS(COEFFS5)
265 CASE(6)
266 CALL PRINT_COEFFS(COEFFS6)
267 END SELECT
268 WRITE(*,*)'Coefs reconstructed by Golem95:'
269 SELECT CASE(RANK)
270 CASE(0)
271 WRITE(*,*)'Constant coef = ',COEFFS0_REC
272 CASE(1)
273 CALL PRINT_COEFFS(COEFFS1_REC)
274 CASE(2)
275 CALL PRINT_COEFFS(COEFFS2_REC)
276 CASE(3)
277 CALL PRINT_COEFFS(COEFFS3_REC)
278 CASE(4)
279 CALL PRINT_COEFFS(COEFFS4_REC)
280 CASE(5)
281 CALL PRINT_COEFFS(COEFFS5_REC)
282 CASE(6)
283 CALL PRINT_COEFFS(COEFFS6_REC)
284 END SELECT
285 STOP
286 ENDIF
287
288 SELECT CASE(NLOOPLINE)
289 CASE(1)
290 WRITE(*,*)'Golem95 cannot handle with tadpole yet'
291 STOP
292 CASE(2)
293 SELECT CASE(RANK)
294 CASE(0)
295 RES_GOLEM = COEFFS0*A20(PINCH)
296 CASE(1)
297 RES_GOLEM = CONTRACT2_1(COEFFS1,PGOLEM,PINCH)
298 CASE(2)
299 RES_GOLEM = CONTRACT2_2(COEFFS2,PGOLEM,PINCH)
300 CASE(3)
301 RES_GOLEM = CONTRACT2_3(COEFFS3,PGOLEM,PINCH)
302 CASE DEFAULT
303 WRITE(*,*)'Golem95 cannot handle with: N,r = ',2,RANK
304 STOP
305 END SELECT
306 CASE(3)
307 SELECT CASE(RANK)
308 CASE(0)
309 RES_GOLEM = COEFFS0*A30(PINCH)
310 CASE(1)
311 RES_GOLEM = CONTRACT3_1(COEFFS1,PGOLEM,PINCH)
312 CASE(2)
313 RES_GOLEM = CONTRACT3_2(COEFFS2,PGOLEM,PINCH)
314 CASE(3)
315 RES_GOLEM = CONTRACT3_3(COEFFS3,PGOLEM,PINCH)
316 CASE(4)
317 RES_GOLEM = CONTRACT3_4(COEFFS4,PGOLEM,PINCH)
318 CASE DEFAULT
319 WRITE(*,*)'Golem95 cannot handle with: N,r = ',3,RANK
320 STOP
321 END SELECT
322 CASE(4)
323 SELECT CASE(RANK)
324 CASE(0)
325 RES_GOLEM = COEFFS0*A40(PINCH)
326 CASE(1)
327 RES_GOLEM = CONTRACT4_1(COEFFS1,PGOLEM,PINCH)
328 CASE(2)
329 RES_GOLEM = CONTRACT4_2(COEFFS2,PGOLEM,PINCH)
330 CASE(3)
331 RES_GOLEM = CONTRACT4_3(COEFFS3,PGOLEM,PINCH)
332 CASE(4)
333 RES_GOLEM = CONTRACT4_4(COEFFS4,PGOLEM,PINCH)
334 CASE(5)
335 RES_GOLEM = CONTRACT4_5(COEFFS5,PGOLEM,PINCH)
336 CASE DEFAULT
337 WRITE(*,*)'Golem95 cannot handle with: N,r = ',4,RANK
338 STOP
339 END SELECT
340 CASE(5)
341 SELECT CASE(RANK)
342 CASE(0)
343 RES_GOLEM = COEFFS0*A50(PINCH)
344 CASE(1)
345 RES_GOLEM = CONTRACT5_1(COEFFS1,PGOLEM,PINCH)
346 CASE(2)
347 RES_GOLEM = CONTRACT5_2(COEFFS2,PGOLEM,PINCH)
348 CASE(3)
349 RES_GOLEM = CONTRACT5_3(COEFFS3,PGOLEM,PINCH)
350 CASE(4)
351 RES_GOLEM = CONTRACT5_4(COEFFS4,PGOLEM,PINCH)
352 CASE(5)
353 RES_GOLEM = CONTRACT5_5(COEFFS5,PGOLEM,PINCH)
354 CASE(6)
355 RES_GOLEM = CONTRACT5_6(COEFFS6,PGOLEM,PINCH)
356 CASE DEFAULT
357 WRITE(*,*)'Golem95 cannot handle with: N,r = ',5,RANK
358 STOP
359 END SELECT
360 CASE(6)
361 SELECT CASE(RANK)
362 CASE(0)
363 RES_GOLEM = COEFFS0*A60(PINCH)
364 CASE(1)
365 RES_GOLEM = CONTRACT6_1(COEFFS1,PGOLEM,PINCH)
366 CASE(2)
367 RES_GOLEM = CONTRACT6_2(COEFFS2,PGOLEM,PINCH)
368 CASE(3)
369 RES_GOLEM = CONTRACT6_3(COEFFS3,PGOLEM,PINCH)
370 CASE(4)
371 RES_GOLEM = CONTRACT6_4(COEFFS4,PGOLEM,PINCH)
372 CASE(5)
373 RES_GOLEM = CONTRACT6_5(COEFFS5,PGOLEM,PINCH)
374 CASE(6)
375 RES_GOLEM = CONTRACT6_6(COEFFS6,PGOLEM,PINCH)
376 CASE DEFAULT
377 WRITE(*,*)'Golem95 cannot handle with: N,r = ',6,RANK
378 STOP
379 END SELECT
380 CASE DEFAULT
381 WRITE(*,*)'Golem95 cannot handle with: N = ',NLOOPLINE
382 STOP
383 END SELECT
384 ENDIF
385
386 RES(1)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%C+2.0*LOG(MU_R)
387 $ *RES_GOLEM%%B+2.0*LOG(MU_R)**2*RES_GOLEM%%A)
388 RES(2)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%B+2.0*LOG(MU_R)
389 $ *RES_GOLEM%%A)
390 RES(3)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%A)
391C WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
392
393C Finally free golem memory and cache
394 CALL EXITGOLEM95()
395
396 END
397
398 FUNCTION GOLEM_LOOPNUM(Q,MU2)
399 USE PRECISION_GOLEM, ONLY: KI
400 REAL(KI),DIMENSION(0:3),INTENT(IN)::Q
401 REAL(KI),INTENT(IN)::MU2
402 COMPLEX(KI)::GOLEM_LOOPNUM
403
404 COMPLEX*16 QQ(0:3),NUM
405 INTEGER I
406
407 DO I=0,3
408 QQ(I)=CMPLX(Q(I),0.0D0,KIND=16)
409 ENDDO
410
411 CALL LOOPNUM(QQ,NUM)
412 GOLEM_LOOPNUM=NUM
413 RETURN
414 END FUNCTION
415
416 SUBROUTINE INITGOLEM()
417C
418C INITIALISATION OF GOLEM
419C
420C
421C MODULE
422C
423 USE PARAMETRE
424C
425C LOCAL VARIABLES
426C
427 REAL*8 THRS
428 LOGICAL EXT_NUM_FOR_R1
429C
430C GLOBAL VARIABLES
431C
432 INCLUDE 'MadLoopParams.inc'
433C ----------
434C BEGIN CODE
435C ----------
436
437C DEFAULT PARAMETERS FOR GOLEM
438C -------------------------------
439C One can chose here to have either just the rational R1 piece
440C or everything but the R2
441 RAT_OR_TOT_PAR = TOT
442
443 END
444
445 SUBROUTINE SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2L)
446C
447C MODULE
448C
449 USE MATRICE_S
450 USE PRECISION_GOLEM, ONLY: KI
451C
452C ARGUMENTS
453C
454 INTEGER NLOOPLINE
455 REAL(KI) PGOLEM(NLOOPLINE,0:3)
456 COMPLEX(KI) M2L(NLOOPLINE)
457
458C
459C GLOBAL VARIABLES
460C
461 INCLUDE 'MadLoopParams.inc'
462 INTEGER CTMODE
463 REAL*8 LSCALE
464 COMMON/CT/LSCALE,CTMODE
465 REAL*8 REF_NORMALIZATION
466C
467C LOCAL VARIABLES
468C
469 INTEGER I,J,K
470 COMPLEX(KI) DIFFSQ
471C ----------
472C BEGIN CODE
473C ----------
474
475 DO I=1,NLOOPLINE
476 DO J=1,NLOOPLINE
477 IF(I.EQ.J)THEN
478 S_MAT(I,J)=-(M2L(I)+M2L(J))
479 ELSE
480 DIFFSQ = (CMPLX(PGOLEM(I,0),0.0_KI,KIND=KI)-CMPLX(PGOLEM(J
481 $ ,0),0.0_KI,KIND=KI))**2
482 DO K=1,3
483 DIFFSQ = DIFFSQ-(CMPLX(PGOLEM(I,K),0.0_KI,KIND=KI)
484 $ -CMPLX(PGOLEM(J,K),0.0_KI,KIND=KI))**2
485 ENDDO
486 S_MAT(I,J)=DIFFSQ-M2L(I)-M2L(J)
487 IF(M2L(I).NE.0.0D0)THEN
488 IF(ABS((DIFFSQ-M2L(I))/M2L(I)).LT.OSTHRES)THEN
489 S_MAT(I,J)=-M2L(J)
490 ENDIF
491 ENDIF
492 IF(M2L(J).NE.0.0D0)THEN
493 IF(ABS((DIFFSQ-M2L(J))/M2L(J)).LT.OSTHRES)THEN
494 S_MAT(I,J)=-M2L(I)
495 ENDIF
496 ENDIF
497C Chose what seems the most appropriate way to compare
498C massless onshellness. (here we pick the energy component)
499 REF_NORMALIZATION=(PGOLEM(I,0)+PGOLEM(J,0))**2
500 IF(REF_NORMALIZATION.NE.0.0D0)THEN
501 IF(ABS(DIFFSQ/REF_NORMALIZATION).LT.OSTHRES)THEN
502 S_MAT(I,J)=-(M2L(I)+M2L(J))
503 ENDIF
504 ENDIF
505 ENDIF
506 ENDDO
507 ENDDO
508
509 END
510
511 FUNCTION COMPARE_COEFS_0(COEFS_A,COEFS_B)
512
513 USE PRECISION_GOLEM, ONLY: KI
514 COMPLEX(KI) COEFS_A,COEFS_B
515 REAL*8 COEF_CHECK_THRS
516 COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
517 REAL*8 DENOM,NUM
518 LOGICAL COMPARE_COEFS_0
519
520 NUM = ABS(COEFS_A-COEFS_B)
521 DENOM = ABS(COEFS_A+COEFS_B)
522 IF(DENOM.GT.0D0)THEN
523 COMPARE_COEFS_0=((NUM/DENOM).LT.COEF_CHECK_THRS)
524 ELSE
525 COMPARE_COEFS_0=(NUM.LT.COEF_CHECK_THRS)
526 ENDIF
527
528 END
529
530 FUNCTION COMPARE_COEFS_1(COEFS_A,COEFS_B)
531
532 USE TENS_REC
533 TYPE(COEFF_TYPE_1)COEFS_A,COEFS_B
534 REAL*8 COEF_CHECK_THRS
535 COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
536 REAL*8 DENOM,NUM
537 LOGICAL COMPARE_COEFS_1
538
539 NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
540 DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
541 $ ))
542
543 IF(DENOM.GT.0D0)THEN
544 COMPARE_COEFS_1=((NUM/DENOM).LT.COEF_CHECK_THRS)
545 ELSE
546 COMPARE_COEFS_1=(NUM.LT.COEF_CHECK_THRS)
547 ENDIF
548
549 END
550
551 FUNCTION COMPARE_COEFS_2(COEFS_A,COEFS_B)
552
553 USE TENS_REC
554 TYPE(COEFF_TYPE_2) COEFS_A,COEFS_B
555 REAL*8 COEF_CHECK_THRS
556 COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
557 REAL*8 DENOM,NUM
558 LOGICAL COMPARE_COEFS_2
559
560 NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
561 $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))
562 DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
563 $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))
564 IF(DENOM.GT.0D0)THEN
565 COMPARE_COEFS_2=((NUM/DENOM).LT.COEF_CHECK_THRS)
566 ELSE
567 COMPARE_COEFS_2=(NUM.LT.COEF_CHECK_THRS)
568 ENDIF
569
570 END
571
572 FUNCTION COMPARE_COEFS_3(COEFS_A,COEFS_B)
573
574 USE TENS_REC
575 TYPE(COEFF_TYPE_3) COEFS_A, COEFS_B
576 REAL*8 COEF_CHECK_THRS
577 COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
578 REAL*8 DENOM,NUM
579 LOGICAL COMPARE_COEFS_3
580
581 NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
582 $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
583 DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
584 $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3
585 $ ))
586 IF(DENOM.GT.0D0)THEN
587 COMPARE_COEFS_3=((NUM/DENOM).LT.COEF_CHECK_THRS)
588 ELSE
589 COMPARE_COEFS_3=(NUM.LT.COEF_CHECK_THRS)
590 ENDIF
591
592 END
593
594 FUNCTION COMPARE_COEFS_4(COEFS_A,COEFS_B)
595
596 USE TENS_REC
597 TYPE(COEFF_TYPE_4) COEFS_A, COEFS_B
598 REAL*8 COEF_CHECK_THRS
599 COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
600 REAL*8 DENOM,NUM
601 LOGICAL COMPARE_COEFS_4
602
603 NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
604 $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3)
605 $ )+SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
606 DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
607 $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3
608 $ ))+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
609 IF(DENOM.GT.0D0)THEN
610 COMPARE_COEFS_4=((NUM/DENOM).LT.COEF_CHECK_THRS)
611 ELSE
612 COMPARE_COEFS_4=(NUM.LT.COEF_CHECK_THRS)
613 ENDIF
614
615 END
616
617 FUNCTION COMPARE_COEFS_5(COEFS_A,COEFS_B)
618
619 USE TENS_REC
620 TYPE(COEFF_TYPE_5) COEFS_A,COEFS_B
621 REAL*8 COEF_CHECK_THRS
622 COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
623 REAL*8 DENOM,NUM
624 LOGICAL COMPARE_COEFS_5
625
626 NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
627 $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3)
628 $ )+SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
629 DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
630 $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3
631 $ ))+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
632 IF(DENOM.GT.0D0)THEN
633 COMPARE_COEFS_5=((NUM/DENOM).LT.COEF_CHECK_THRS)
634 ELSE
635 COMPARE_COEFS_5=(NUM.LT.COEF_CHECK_THRS)
636 ENDIF
637
638 END
639
640 FUNCTION COMPARE_COEFS_6(COEFS_A,COEFS_B)
641
642 USE TENS_REC
643 TYPE(COEFF_TYPE_6) COEFS_A,COEFS_B
644 REAL*8 COEF_CHECK_THRS
645 COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
646 REAL*8 DENOM,NUM
647 LOGICAL COMPARE_COEFS_6
648
649 NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
650 $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3)
651 $ )+SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
652 DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1
653 $ ))+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3
654 $ ))+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
655 IF(DENOM.GT.0D0)THEN
656 COMPARE_COEFS_6=((NUM/DENOM).LT.COEF_CHECK_THRS)
657 ELSE
658 COMPARE_COEFS_6=(NUM.LT.COEF_CHECK_THRS)
659 ENDIF
660
661 END
662
663
664 SUBROUTINE FILL_GOLEM_COEFFS_0(ML_COEFS,GOLEM_COEFS)
665 USE PRECISION_GOLEM, ONLY: KI
666 INCLUDE 'coef_specs.inc'
667 COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
668 COMPLEX(KI) GOLEM_COEFS
669 GOLEM_COEFS=ML_COEFS(0)
670 END
671
672 SUBROUTINE FILL_GOLEM_COEFFS_1(ML_COEFS,GOLEM_COEFS)
673 USE TENS_REC, ONLY: COEFF_TYPE_1
674 INCLUDE 'coef_specs.inc'
675 COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
676 TYPE(COEFF_TYPE_1) GOLEM_COEFS
677C Constant coefficient
678 GOLEM_COEFS%%C0=ML_COEFS(0)
679C Coefficient q(0)
680 GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
681C Coefficient q(1)
682 GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
683C Coefficient q(2)
684 GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
685C Coefficient q(3)
686 GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
687 END
688
689 SUBROUTINE FILL_GOLEM_COEFFS_2(ML_COEFS,GOLEM_COEFS)
690 USE TENS_REC, ONLY: COEFF_TYPE_2
691 INCLUDE 'coef_specs.inc'
692 COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
693 TYPE(COEFF_TYPE_2) GOLEM_COEFS
694C Constant coefficient
695 GOLEM_COEFS%%C0=ML_COEFS(0)
696C Coefficient q(0)
697 GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
698C Coefficient q(0)^2
699 GOLEM_COEFS%%C1(1,2)= ML_COEFS(5)
700C Coefficient q(1)
701 GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
702C Coefficient q(1)^2
703 GOLEM_COEFS%%C1(2,2)= ML_COEFS(7)
704C Coefficient q(2)
705 GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
706C Coefficient q(2)^2
707 GOLEM_COEFS%%C1(3,2)= ML_COEFS(10)
708C Coefficient q(3)
709 GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
710C Coefficient q(3)^2
711 GOLEM_COEFS%%C1(4,2)= ML_COEFS(14)
712C Coefficient q(0)*q(1)
713 GOLEM_COEFS%%C2(1,1)= ML_COEFS(6)
714C Coefficient q(0)*q(2)
715 GOLEM_COEFS%%C2(2,1)= ML_COEFS(8)
716C Coefficient q(0)*q(3)
717 GOLEM_COEFS%%C2(3,1)= ML_COEFS(11)
718C Coefficient q(1)*q(2)
719 GOLEM_COEFS%%C2(4,1)= ML_COEFS(9)
720C Coefficient q(1)*q(3)
721 GOLEM_COEFS%%C2(5,1)= ML_COEFS(12)
722C Coefficient q(2)*q(3)
723 GOLEM_COEFS%%C2(6,1)= ML_COEFS(13)
724 END
725
726 SUBROUTINE FILL_GOLEM_COEFFS_3(ML_COEFS,GOLEM_COEFS)
727 USE TENS_REC, ONLY: COEFF_TYPE_3
728 INCLUDE 'coef_specs.inc'
729 COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
730 TYPE(COEFF_TYPE_3) GOLEM_COEFS
731C Dummy routine for FILL_GOLEM_COEFS_3
732 STOP 'ERROR: 3 > 2'
733 END
734
735 SUBROUTINE FILL_GOLEM_COEFFS_4(ML_COEFS,GOLEM_COEFS)
736 USE TENS_REC, ONLY: COEFF_TYPE_4
737 INCLUDE 'coef_specs.inc'
738 COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
739 TYPE(COEFF_TYPE_4) GOLEM_COEFS
740C Dummy routine for FILL_GOLEM_COEFS_4
741 STOP 'ERROR: 4 > 2'
742 END
743
744 SUBROUTINE FILL_GOLEM_COEFFS_5(ML_COEFS,GOLEM_COEFS)
745 USE TENS_REC, ONLY: COEFF_TYPE_5
746 INCLUDE 'coef_specs.inc'
747 COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
748 TYPE(COEFF_TYPE_5) GOLEM_COEFS
749C Dummy routine for FILL_GOLEM_COEFS_5
750 STOP 'ERROR: 5 > 2'
751 END
752
753 SUBROUTINE FILL_GOLEM_COEFFS_6(ML_COEFS,GOLEM_COEFS)
754 USE TENS_REC, ONLY: COEFF_TYPE_6
755 INCLUDE 'coef_specs.inc'
756 COMPLEX*16 ML_COEFS(0:LOOP_MAXCOEFS-1)
757 TYPE(COEFF_TYPE_6) GOLEM_COEFS
758C Dummy routine for FILL_GOLEM_COEFS_6
759 STOP 'ERROR: 6 > 2'
760 END
0761
=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f 1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f 2016-02-24 13:54:50 +0000
@@ -0,0 +1,423 @@
1 SUBROUTINE TIRLOOP(I_SQSO,I_LOOPGROUP,I_LIB,NLOOPLINE,PL,M2L
2 $ ,RANK,RES,STABLE)
3C
4C Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
5C By the MadGraph5_aMC@NLO Development Team
6C Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
7C
8C Interface between MG5 and TIR.
9C
10C Process: d~ u > w+ WEIGHTED=2 QED=1 [ all = QCD ]
11C Process: s~ c > w+ WEIGHTED=2 QED=1 [ all = QCD ]
12C
13C
14C CONSTANTS
15C
16 INTEGER NLOOPGROUPS
17 PARAMETER (NLOOPGROUPS=1)
18C These are constants related to the split orders
19 INTEGER NSQUAREDSO
20 PARAMETER (NSQUAREDSO=1)
21 INTEGER LOOPMAXCOEFS
22 PARAMETER (LOOPMAXCOEFS=15)
23 INTEGER NEXTERNAL
24 PARAMETER (NEXTERNAL=3)
25 LOGICAL CHECKPCONSERVATION
26 PARAMETER (CHECKPCONSERVATION=.TRUE.)
27 REAL*8 NORMALIZATION
28 PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
29 $ *2))
30C
31C ARGUMENTS
32C
33 INTEGER I_SQSO,I_LOOPGROUP,I_LIB
34 INTEGER NLOOPLINE, RANK
35 REAL*8 PL(0:3,NLOOPLINE)
36 REAL*8 PCT(0:3,0:NLOOPLINE-1)
37 REAL*8 PDEN(0:3,NLOOPLINE-1)
38 COMPLEX*16 M2L(NLOOPLINE)
39 COMPLEX*16 M2LCT(0:NLOOPLINE-1)
40 COMPLEX*16 RES(3)
41 LOGICAL STABLE
42C
43C LOCAL VARIABLES
44C
45 INTEGER I, J, K
46 INTEGER NLOOPCOEFS
47 LOGICAL CTINIT, TIRINIT, GOLEMINIT
48 COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT
49
50C This variable will be used to detect changes in the TIR library
51C used so as to force the reset of the TIR filter.
52 INTEGER LAST_LIB_USED
53 DATA LAST_LIB_USED/-1/
54
55 COMPLEX*16 TIRCOEFS(0:LOOPMAXCOEFS-1,3)
56 COMPLEX*16 PJCOEFS(0:LOOPMAXCOEFS-1,3)
57C
58C EXTERNAL FUNCTIONS
59C
60C
61C GLOBAL VARIABLES
62C
63 INCLUDE 'MadLoopParams.inc'
64 INCLUDE 'coupl.inc'
65 INTEGER CTMODE
66 REAL*8 LSCALE
67 COMMON/CT/LSCALE,CTMODE
68
69C The argument ILIB is the TIR library to be used for that
70C specific library.
71 INTEGER LIBINDEX
72 COMMON/I_LIB/LIBINDEX
73
74 COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
75 COMMON/LCOEFS/LOOPCOEFS
76C ----------
77C BEGIN CODE
78C ----------
79
80C Initialize for the very first time ML is called LAST_ILIB with
81C the first ILIB used.
82 IF(LAST_LIB_USED.EQ.-1) THEN
83 LAST_LIB_USED = MLREDUCTIONLIB(LIBINDEX)
84 ELSE
85C We changed the TIR library so we must refresh the cache.
86 IF(MLREDUCTIONLIB(LIBINDEX).NE.LAST_LIB_USED) THEN
87 LAST_LIB_USED = MLREDUCTIONLIB(LIBINDEX)
88 CALL CLEAR_TIR_CACHE()
89 ENDIF
90 ENDIF
91
92 IF (MLREDUCTIONLIB(I_LIB).EQ.4) THEN
93C Using Golem95
94C PDEN is dummy for Golem95 so we just initialize it to zero
95C here so as to use it for the function SWITCHORDER
96 DO I=0,3
97 DO J=1,NLOOPLINE-1
98 PDEN(I,J)=0.0D0
99 ENDDO
100 ENDDO
101 CALL SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN,M2L)
102 CALL GOLEMLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
103 RETURN
104 ENDIF
105
106C INITIALIZE TIR IF NEEDED
107 IF (TIRINIT) THEN
108 TIRINIT=.FALSE.
109 CALL INITTIR()
110 ENDIF
111
112C CONVERT THE MASSES TO BE COMPLEX
113 DO I=1,NLOOPLINE
114 M2LCT(I-1)=M2L(I)
115 ENDDO
116
117C CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
118 DO I=0,3
119 DO J=0,(NLOOPLINE-1)
120 PCT(I,J)=0.D0
121 ENDDO
122 ENDDO
123 DO I=0,3
124 DO J=1,NLOOPLINE
125 PCT(I,0)=PCT(I,0)+PL(I,J)
126 ENDDO
127 ENDDO
128 IF (CHECKPCONSERVATION) THEN
129 IF (PCT(0,0).GT.1.D-6) THEN
130 WRITE(*,*) 'energy is not conserved ',PCT(0,0)
131 STOP 'energy is not conserved'
132 ELSEIF (PCT(1,0).GT.1.D-6) THEN
133 WRITE(*,*) 'px is not conserved ',PCT(1,0)
134 STOP 'px is not conserved'
135 ELSEIF (PCT(2,0).GT.1.D-6) THEN
136 WRITE(*,*) 'py is not conserved ',PCT(2,0)
137 STOP 'py is not conserved'
138 ELSEIF (PCT(3,0).GT.1.D-6) THEN
139 WRITE(*,*) 'pz is not conserved ',PCT(3,0)
140 STOP 'pz is not conserved'
141 ENDIF
142 ENDIF
143 DO I=0,3
144 DO J=1,(NLOOPLINE-1)
145 DO K=1,J
146 PCT(I,J)=PCT(I,J)+PL(I,K)
147 ENDDO
148 ENDDO
149 ENDDO
150
151 DO I=0,3
152 DO J=1,(NLOOPLINE-1)
153 PDEN(I,J)=PCT(I,J)
154 ENDDO
155 ENDDO
156C NUMBER OF INDEPEDENT LOOPCOEFS FOR RANK=RANK
157 NLOOPCOEFS=0
158 DO I=0,RANK
159 NLOOPCOEFS=NLOOPCOEFS+(3+I)*(2+I)*(1+I)/6
160 ENDDO
161 SELECT CASE(MLREDUCTIONLIB(I_LIB))
162 CASE(2)
163C PJFry++
164 WRITE(*,*) 'ERROR:: PJFRY++ is not interfaced.'
165 STOP
166 CASE(3)
167C IREGI
168 CALL IMLOOP(CTMODE,IREGIMODE,NLOOPLINE,LOOPMAXCOEFS,RANK,PDEN
169 $ ,M2L,MU_R,PJCOEFS,STABLE)
170C CONVERT TO MADLOOP CONVENTION
171 CALL CONVERT_IREGI_COEFFS(RANK,PJCOEFS,TIRCOEFS)
172 END SELECT
173 DO I=1,3
174 RES(I)=(0.0D0,0.0D0)
175 DO J=0,NLOOPCOEFS-1
176 RES(I)=RES(I)+LOOPCOEFS(J,I_SQSO,I_LOOPGROUP)*TIRCOEFS(J,I)
177 ENDDO
178 ENDDO
179 RES(1)=NORMALIZATION*2.0D0*DBLE(RES(1))
180 RES(2)=NORMALIZATION*2.0D0*DBLE(RES(2))
181 RES(3)=NORMALIZATION*2.0D0*DBLE(RES(3))
182C WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
183 END
184
185 SUBROUTINE SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN,M2L)
186 IMPLICIT NONE
187
188 INTEGER CTMODE,NLOOPLINE
189
190 REAL*8 PL(0:3,NLOOPLINE)
191 REAL*8 PDEN(0:3,NLOOPLINE-1)
192 COMPLEX*16 M2L(NLOOPLINE)
193 REAL*8 NEW_PL(0:3,NLOOPLINE)
194 REAL*8 NEW_PDEN(0:3,NLOOPLINE-1)
195 COMPLEX*16 NEW_M2L(NLOOPLINE)
196
197 INTEGER I,J,K
198
199 IF (CTMODE.NE.2.AND.CTMODE.NE.4) THEN
200 RETURN
201 ENDIF
202
203 IF (NLOOPLINE.LE.2) THEN
204 RETURN
205 ENDIF
206
207 DO I=1,NLOOPLINE-1
208 DO J=0,3
209 NEW_PDEN(J,NLOOPLINE-I) = PDEN(J,I)
210 ENDDO
211 ENDDO
212 DO I=1,NLOOPLINE-1
213 DO J=0,3
214 PDEN(J,I) = NEW_PDEN(J,I)
215 ENDDO
216 ENDDO
217
218 DO I=2,NLOOPLINE
219 NEW_M2L(I)=M2L(NLOOPLINE-I+2)
220 ENDDO
221 DO I=2,NLOOPLINE
222 M2L(I)=NEW_M2L(I)
223 ENDDO
224
225
226 DO I=1,NLOOPLINE
227 DO J=0,3
228 NEW_PL(J,I) = -PL(J,NLOOPLINE+1-I)
229 ENDDO
230 ENDDO
231 DO I=1,NLOOPLINE
232 DO J=0,3
233 PL(J,I) = NEW_PL(J,I)
234 ENDDO
235 ENDDO
236
237 END
238
239 SUBROUTINE INITTIR()
240C
241C INITIALISATION OF TIR
242C
243C LOCAL VARIABLES
244C
245 REAL*8 THRS
246 LOGICAL EXT_NUM_FOR_R1
247C
248C GLOBAL VARIABLES
249C
250 INCLUDE 'MadLoopParams.inc'
251 LOGICAL CTINIT, TIRINIT, GOLEMINIT
252 COMMON/REDUCTIONCODEINIT/CTINIT,TIRINIT,GOLEMINIT
253
254C ----------
255C BEGIN CODE
256C ----------
257
258C DEFAULT PARAMETERS FOR TIR
259C -------------------------------
260C THRS1 IS THE PRECISION LIMIT BELOW WHICH THE MP ROUTINES
261C ACTIVATES
262C USE THE SAME MADLOOP PARAMETER IN CUTTOOLS AND TIR
263C IT IS NECESSARY TO INITIALIZE CT BECAUSE IREGI USES THE VERSION
264C OF ONELOOP
265C FROM CUTTOOLS LIBRARY
266 THRS=CTSTABTHRES
267C LOOPLIB SET WHAT LIBRARY CT USES
268C 1 -> LOOPTOOLS
269C 2 -> AVH
270C 3 -> QCDLOOP
271 LOOPLIB=CTLOOPLIBRARY
272 CALL INITIREGI(IREGIRECY,LOOPLIB,1D-6)
273C The initialization below is for CT v1.9.2+
274 IF (CTINIT) THEN
275 CTINIT=.FALSE.
276 CALL INITCT()
277 ENDIF
278 END
279
280 SUBROUTINE CHOOSE_LOOPLIB(LIBINDEX,NLOOPLINE,RANK,COMPLEX_MASS
281 $ ,DOING_QP,I_LIB)
282C
283C CHOOSE THE CORRECT LOOP LIB
284C Example:
285C MLReductionLib=3|2|1 and LIBINDEX=3
286C IF THE LOOP IS BEYOND THE SCOPE OF LOOP LIB MLReductionLib(3)=1
287C USE LIBINDEX=1, and LIBINDEX=2 ...
288C IF IT IS STILL NOT GOOD,STOP
289C
290 IMPLICIT NONE
291C
292C CONSTANTS
293C
294 INTEGER NLOOPLIB,QP_NLOOPLIB
295 PARAMETER (NLOOPLIB=4,QP_NLOOPLIB=1)
296C
297C ARGUMENTS
298C
299 INTEGER LIBINDEX,NLOOPLINE,RANK,I_LIB
300 LOGICAL COMPLEX_MASS,DOING_QP
301C
302C LOCAL VARIABLES
303C
304 INTEGER I,J_LIB,LIBNUM,SELECT_LIBINDEX
305 LOGICAL LPASS
306C
307C GLOBAL VARIABLES
308C
309 INCLUDE 'MadLoopParams.inc'
310C TILL NOW, ONLY CUTTOOLS PROVIDE QP
311 LOGICAL QP_TOOLS_AVAILABLE
312 INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB)
313 COMMON/LOOP_TOOLS/QP_TOOLS_AVAILABLE,INDEX_QP_TOOLS
314C ----------
315C BEGIN CODE
316C ----------
317
318 IF(DOING_QP)THEN
319C QP EVALUATION, ONLY CUTTOOLS
320 IF(.NOT.QP_TOOLS_AVAILABLE)THEN
321 STOP 'No qp tools available, please make sure MLReductionLi'
322 $ //'b is correct'
323 ENDIF
324 J_LIB=0
325 SELECT_LIBINDEX=LIBINDEX
326 DO WHILE(J_LIB.EQ.0)
327 DO I=1,QP_NLOOPLIB
328 IF(INDEX_QP_TOOLS(I).EQ.SELECT_LIBINDEX)THEN
329 J_LIB=I
330 EXIT
331 ENDIF
332 ENDDO
333 IF(J_LIB.EQ.0)THEN
334 SELECT_LIBINDEX=SELECT_LIBINDEX+1
335 IF(SELECT_LIBINDEX.GT.NLOOPLIB.OR.MLREDUCTIONLIB(SELECT_LIB
336 $ INDEX).EQ.0)SELECT_LIBINDEX=1
337 ENDIF
338 ENDDO
339 I=J_LIB
340 I_LIB=SELECT_LIBINDEX
341 LIBNUM=MLREDUCTIONLIB(I_LIB)
342 DO
343 CALL DETECT_LOOPLIB(LIBNUM,NLOOPLINE,RANK,COMPLEX_MASS,LPASS)
344 IF(LPASS)EXIT
345 I=I+1
346 IF(I.GT.QP_NLOOPLIB.AND.INDEX_QP_TOOLS(I).EQ.0)THEN
347 I=1
348 ENDIF
349 IF(I.EQ.J_LIB)THEN
350 STOP 'No qp loop library can deal with this integral'
351 ENDIF
352 I_LIB=INDEX_QP_TOOLS(I)
353 LIBNUM=MLREDUCTIONLIB(I_LIB)
354 ENDDO
355 ELSE
356C DP EVALUATION
357 I_LIB=LIBINDEX
358 LIBNUM=MLREDUCTIONLIB(I_LIB)
359 DO
360 CALL DETECT_LOOPLIB(LIBNUM,NLOOPLINE,RANK,COMPLEX_MASS,LPASS)
361 IF(LPASS)EXIT
362 I_LIB=I_LIB+1
363 IF(I_LIB.GT.NLOOPLIB.OR.MLREDUCTIONLIB(I_LIB).EQ.0)THEN
364 I_LIB=1
365 ENDIF
366 IF(I_LIB.EQ.LIBINDEX)THEN
367 STOP 'No dp loop library can deal with this integral'
368 ENDIF
369 LIBNUM=MLREDUCTIONLIB(I_LIB)
370 ENDDO
371 ENDIF
372 RETURN
373 END
374
375 SUBROUTINE CLEAR_TIR_CACHE()
376C No TIR caching implemented, this is dummy. (The subroutine is
377C kept as it might be called by the MC).
378 CONTINUE
379 END SUBROUTINE
380
381
382
383 SUBROUTINE CONVERT_IREGI_COEFFS(RANK,IREGICOEFS,TIRCOEFS)
384C GLOABLE VARIABLES
385 INCLUDE 'coef_specs.inc'
386C ARGUMENTS
387 INTEGER RANK
388 COMPLEX*16 IREGICOEFS(0:LOOP_MAXCOEFS-1,3)
389 COMPLEX*16 TIRCOEFS(0:LOOP_MAXCOEFS-1,3)
390C Reduction Coefficient 1
391 TIRCOEFS(0,1:3)=IREGICOEFS(0,1:3)
392 IF(RANK.LE.0)RETURN
393C Reduction Coefficient q(0)
394 TIRCOEFS(1,1:3)=IREGICOEFS(1,1:3)
395C Reduction Coefficient q(1)
396 TIRCOEFS(2,1:3)=IREGICOEFS(2,1:3)
397C Reduction Coefficient q(2)
398 TIRCOEFS(3,1:3)=IREGICOEFS(3,1:3)
399C Reduction Coefficient q(3)
400 TIRCOEFS(4,1:3)=IREGICOEFS(4,1:3)
401 IF(RANK.LE.1)RETURN
402C Reduction Coefficient q(0)^2
403 TIRCOEFS(5,1:3)=IREGICOEFS(5,1:3)
404C Reduction Coefficient q(0)*q(1)
405 TIRCOEFS(6,1:3)=IREGICOEFS(6,1:3)
406C Reduction Coefficient q(1)^2
407 TIRCOEFS(7,1:3)=IREGICOEFS(9,1:3)
408C Reduction Coefficient q(0)*q(2)
409 TIRCOEFS(8,1:3)=IREGICOEFS(7,1:3)
410C Reduction Coefficient q(1)*q(2)
411 TIRCOEFS(9,1:3)=IREGICOEFS(10,1:3)
412C Reduction Coefficient q(2)^2
413 TIRCOEFS(10,1:3)=IREGICOEFS(12,1:3)
414C Reduction Coefficient q(0)*q(3)
415 TIRCOEFS(11,1:3)=IREGICOEFS(8,1:3)
416C Reduction Coefficient q(1)*q(3)
417 TIRCOEFS(12,1:3)=IREGICOEFS(11,1:3)
418C Reduction Coefficient q(2)*q(3)
419 TIRCOEFS(13,1:3)=IREGICOEFS(13,1:3)
420C Reduction Coefficient q(3)^2
421 TIRCOEFS(14,1:3)=IREGICOEFS(14,1:3)
422 IF(RANK.LE.2)RETURN
423 END
0424
=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f 1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f 2016-02-24 13:54:50 +0000
@@ -0,0 +1,540 @@
1 SUBROUTINE SMATRIX(P,ANS_SUMMED)
2C
3C Simple routine wrapper to provide the same interface for
4C backward compatibility for usage without split orders.
5C
6C
7C CONSTANTS
8C
9 INTEGER NEXTERNAL
10 PARAMETER (NEXTERNAL=3)
11 INTEGER NSQAMPSO
12 PARAMETER (NSQAMPSO=1)
13C
14C ARGUMENTS
15C
16 REAL*8 P(0:3,NEXTERNAL), ANS_SUMMED
17C
18C VARIABLES
19C
20 INTEGER I
21 REAL*8 ANS(0:NSQAMPSO)
22C
23C BEGIN CODE
24C
25 CALL SMATRIX_SPLITORDERS(P,ANS)
26 ANS_SUMMED=ANS(0)
27
28 END
29
30 SUBROUTINE SMATRIXHEL(P,HEL,ANS)
31 IMPLICIT NONE
32C
33C CONSTANT
34C
35 INTEGER NEXTERNAL
36 PARAMETER (NEXTERNAL=3)
37 INTEGER NCOMB
38 PARAMETER ( NCOMB=12)
39C
40C ARGUMENTS
41C
42 REAL*8 P(0:3,NEXTERNAL),ANS
43 INTEGER HEL
44C
45C GLOBAL VARIABLES
46C
47 INTEGER USERHEL
48 COMMON/HELUSERCHOICE/USERHEL
49C ----------
50C BEGIN CODE
51C ----------
52 USERHEL=HEL
53 CALL SMATRIX(P,ANS)
54 USERHEL=-1
55
56 END
57
58 SUBROUTINE SMATRIX_SPLITORDERS(P,ANS)
59C
60C Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
61C By the MadGraph5_aMC@NLO Development Team
62C Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
63C
64C MadGraph StandAlone Version
65C
66C Returns amplitude squared summed/avg over colors
67C and helicities
68C for the point in phase space P(0:3,NEXTERNAL)
69C
70C Process: d~ u > w+ WEIGHTED=2 QED=1 [ all = QCD ]
71C Process: s~ c > w+ WEIGHTED=2 QED=1 [ all = QCD ]
72C
73 IMPLICIT NONE
74C
75C CONSTANTS
76C
77 INTEGER NEXTERNAL
78 PARAMETER (NEXTERNAL=3)
79 INTEGER NCOMB
80 PARAMETER ( NCOMB=12)
81 INTEGER NSQAMPSO
82 PARAMETER (NSQAMPSO=1)
83 INTEGER HELAVGFACTOR
84 PARAMETER (HELAVGFACTOR=4)
85 LOGICAL CHOSEN_SO_CONFIGS(NSQAMPSO)
86 DATA CHOSEN_SO_CONFIGS/.TRUE./
87 COMMON/CHOSEN_BORN_SQSO/CHOSEN_SO_CONFIGS
88C
89C ARGUMENTS
90C
91 REAL*8 P(0:3,NEXTERNAL),ANS(0:NSQAMPSO)
92C
93C LOCAL VARIABLES
94C
95 INTEGER NHEL(NEXTERNAL,NCOMB),NTRY
96 REAL*8 T(NSQAMPSO), BUFF
97 INTEGER IHEL,IDEN, I
98 INTEGER JC(NEXTERNAL)
99 LOGICAL GOODHEL(NCOMB)
100 DATA NTRY/0/
101 DATA GOODHEL/NCOMB*.FALSE./
102 DATA (NHEL(I, 1),I=1,3) /-1, 1,-1/
103 DATA (NHEL(I, 2),I=1,3) /-1, 1, 0/
104 DATA (NHEL(I, 3),I=1,3) /-1, 1, 1/
105 DATA (NHEL(I, 4),I=1,3) /-1,-1,-1/
106 DATA (NHEL(I, 5),I=1,3) /-1,-1, 0/
107 DATA (NHEL(I, 6),I=1,3) /-1,-1, 1/
108 DATA (NHEL(I, 7),I=1,3) / 1, 1,-1/
109 DATA (NHEL(I, 8),I=1,3) / 1, 1, 0/
110 DATA (NHEL(I, 9),I=1,3) / 1, 1, 1/
111 DATA (NHEL(I, 10),I=1,3) / 1,-1,-1/
112 DATA (NHEL(I, 11),I=1,3) / 1,-1, 0/
113 DATA (NHEL(I, 12),I=1,3) / 1,-1, 1/
114 DATA IDEN/36/
115C
116C GLOBAL VARIABLES
117C
118 INTEGER USERHEL
119 DATA USERHEL/-1/
120 COMMON/HELUSERCHOICE/USERHEL
121
122C ----------
123C BEGIN CODE
124C ----------
125 NTRY=NTRY+1
126 DO IHEL=1,NEXTERNAL
127 JC(IHEL) = +1
128 ENDDO
129 DO I=1,NSQAMPSO
130 ANS(I) = 0D0
131 ENDDO
132 DO IHEL=1,NCOMB
133 IF (USERHEL.EQ.-1.OR.USERHEL.EQ.IHEL) THEN
134 IF (GOODHEL(IHEL) .OR. NTRY .LT. 2 .OR.USERHEL.NE.-1) THEN
135 CALL MATRIX(P ,NHEL(1,IHEL),JC(1), T)
136 BUFF=0D0
137 DO I=1,NSQAMPSO
138 ANS(I)=ANS(I)+T(I)
139 BUFF=BUFF+T(I)
140 ENDDO
141 IF (BUFF .NE. 0D0 .AND. .NOT. GOODHEL(IHEL)) THEN
142 GOODHEL(IHEL)=.TRUE.
143 ENDIF
144 ENDIF
145 ENDIF
146 ENDDO
147 ANS(0)=0.0D0
148 DO I=1,NSQAMPSO
149 ANS(I)=ANS(I)/DBLE(IDEN)
150 IF (CHOSEN_SO_CONFIGS(I)) THEN
151 ANS(0)=ANS(0)+ANS(I)
152 ENDIF
153 ENDDO
154 IF(USERHEL.NE.-1) THEN
155 ANS(0)=ANS(0)*HELAVGFACTOR
156 DO I=1,NSQAMPSO
157 ANS(I)=ANS(I)*HELAVGFACTOR
158 ENDDO
159 ENDIF
160 END
161
162 SUBROUTINE SMATRIXHEL_SPLITORDERS(P,HEL,ANS)
163 IMPLICIT NONE
164C
165C CONSTANT
166C
167 INTEGER NEXTERNAL
168 PARAMETER (NEXTERNAL=3)
169 INTEGER NCOMB
170 PARAMETER ( NCOMB=12)
171 INTEGER NSQAMPSO
172 PARAMETER (NSQAMPSO=1)
173C
174C ARGUMENTS
175C
176 REAL*8 P(0:3,NEXTERNAL),ANS(0:NSQAMPSO)
177 INTEGER HEL
178C
179C GLOBAL VARIABLES
180C
181 INTEGER USERHEL
182 COMMON/HELUSERCHOICE/USERHEL
183C ----------
184C BEGIN CODE
185C ----------
186 USERHEL=HEL
187 CALL SMATRIX_SPLITORDERS(P,ANS)
188 USERHEL=-1
189
190 END
191
192 SUBROUTINE MATRIX(P,NHEL,IC,RES)
193C
194C Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
195C By the MadGraph5_aMC@NLO Development Team
196C Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
197C
198C Returns amplitude squared summed/avg over colors
199C for the point with external lines W(0:6,NEXTERNAL)
200C
201C Process: d~ u > w+ WEIGHTED=2 QED=1 [ all = QCD ]
202C Process: s~ c > w+ WEIGHTED=2 QED=1 [ all = QCD ]
203C
204 IMPLICIT NONE
205C
206C CONSTANTS
207C
208 INTEGER NGRAPHS
209 PARAMETER (NGRAPHS=1)
210 INTEGER NEXTERNAL
211 PARAMETER (NEXTERNAL=3)
212 INTEGER NWAVEFUNCS, NCOLOR
213 PARAMETER (NWAVEFUNCS=3, NCOLOR=1)
214 INTEGER NAMPSO, NSQAMPSO
215 PARAMETER (NAMPSO=1, NSQAMPSO=1)
216 REAL*8 ZERO
217 PARAMETER (ZERO=0D0)
218 COMPLEX*16 IMAG1
219 PARAMETER (IMAG1=(0D0,1D0))
220C
221C ARGUMENTS
222C
223 REAL*8 P(0:3,NEXTERNAL)
224 INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
225 REAL*8 RES(NSQAMPSO)
226C
227C LOCAL VARIABLES
228C
229 INTEGER I,J,M,N
230 COMPLEX*16 ZTEMP
231 REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
232 COMPLEX*16 AMP(NGRAPHS)
233 COMPLEX*16 JAMP(NCOLOR,NAMPSO)
234 COMPLEX*16 W(20,NWAVEFUNCS)
235 COMPLEX*16 DUM0,DUM1
236 DATA DUM0, DUM1/(0D0, 0D0), (1D0, 0D0)/
237C
238C FUNCTION
239C
240 INTEGER SQSOINDEX
241C
242C GLOBAL VARIABLES
243C
244 INCLUDE 'coupl.inc'
245C
246C COLOR DATA
247C
248 DATA DENOM(1)/1/
249 DATA (CF(I, 1),I= 1, 1) / 3/
250C 1 T(1,2)
251C ----------
252C BEGIN CODE
253C ----------
254 CALL OXXXXX(P(0,1),ZERO,NHEL(1),-1*IC(1),W(1,1))
255 CALL IXXXXX(P(0,2),ZERO,NHEL(2),+1*IC(2),W(1,2))
256 CALL VXXXXX(P(0,3),MDL_MW,NHEL(3),+1*IC(3),W(1,3))
257C Amplitude(s) for diagram number 1
258 CALL FFV2_0(W(1,2),W(1,1),W(1,3),GC_47,AMP(1))
259C JAMPs contributing to orders QCD=0
260 JAMP(1,1)=-AMP(1)
261
262 RES = 0.D0
263 DO M = 1, NAMPSO
264 DO I = 1, NCOLOR
265 ZTEMP = (0.D0,0.D0)
266 DO J = 1, NCOLOR
267 ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M)
268 ENDDO
269 DO N = 1, NAMPSO
270 RES(SQSOINDEX(M,N)) = RES(SQSOINDEX(M,N)) + ZTEMP
271 $ *DCONJG(JAMP(I,N))/DENOM(I)
272 ENDDO
273 ENDDO
274 ENDDO
275
276 END
277
278 SUBROUTINE GET_ME(P, ALPHAS, NHEL ,ANS)
279 IMPLICIT NONE
280C
281C CONSTANT
282C
283 INTEGER NEXTERNAL
284 PARAMETER (NEXTERNAL=3)
285C
286C ARGUMENTS
287C
288 REAL*8 P(0:3,NEXTERNAL),ANS
289 INTEGER NHEL
290 DOUBLE PRECISION ALPHAS
291 REAL*8 PI
292CF2PY INTENT(OUT) :: ANS
293CF2PY INTENT(IN) :: NHEL
294CF2PY INTENT(IN) :: P(0:3,NEXTERNAL)
295CF2PY INTENT(IN) :: ALPHAS
296C ROUTINE FOR F2PY to read the benchmark point.
297C the include file with the values of the parameters and masses
298 INCLUDE 'coupl.inc'
299
300 PI = 3.141592653589793D0
301 G = 2* DSQRT(ALPHAS*PI)
302 CALL UPDATE_AS_PARAM()
303 IF (NHEL.NE.0) THEN
304 CALL SMATRIXHEL(P, NHEL, ANS)
305 ELSE
306 CALL SMATRIX(P, ANS)
307 ENDIF
308 RETURN
309 END
310
311 SUBROUTINE INITIALISE(PATH)
312C ROUTINE FOR F2PY to read the benchmark point.
313 IMPLICIT NONE
314 CHARACTER*180 PATH
315CF2PY INTENT(IN) :: PATH
316 CALL SETPARA(PATH) !first call to setup the paramaters
317 RETURN
318 END
319
320
321C Set of functions to handle the array indices of the split orders
322
323
324 INTEGER FUNCTION SQSOINDEX(ORDERINDEXA, ORDERINDEXB)
325C
326C This functions plays the role of the interference matrix. It can
327C be hardcoded or
328C made more elegant using hashtables if its execution speed ever
329C becomes a relevant
330C factor. From two split order indices, it return the corresponding
331C index in the squared
332C order canonical ordering.
333C
334C CONSTANTS
335C
336
337 INTEGER NSO, NSQUAREDSO, NAMPSO
338 PARAMETER (NSO=1, NSQUAREDSO=1, NAMPSO=1)
339C
340C ARGUMENTS
341C
342 INTEGER ORDERINDEXA, ORDERINDEXB
343C
344C LOCAL VARIABLES
345C
346 INTEGER I, SQORDERS(NSO)
347 INTEGER AMPSPLITORDERS(NAMPSO,NSO)
348 DATA (AMPSPLITORDERS( 1,I),I= 1, 1) / 0/
349 COMMON/AMPSPLITORDERS/AMPSPLITORDERS
350C
351C FUNCTION
352C
353 INTEGER SOINDEX_FOR_SQUARED_ORDERS
354C
355C BEGIN CODE
356C
357 DO I=1,NSO
358 SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I)+AMPSPLITORDERS(ORDERI
359 $ NDEXB,I)
360 ENDDO
361 SQSOINDEX=SOINDEX_FOR_SQUARED_ORDERS(SQORDERS)
362 END
363
364 INTEGER FUNCTION SOINDEX_FOR_SQUARED_ORDERS(ORDERS)
365C
366C This functions returns the integer index identifying the squared
367C split orders list passed in argument which corresponds to the
368C values of the following list of couplings (and in this order).
369C ['QCD']
370C
371C CONSTANTS
372C
373 INTEGER NSO, NSQSO, NAMPSO
374 PARAMETER (NSO=1, NSQSO=1, NAMPSO=1)
375C
376C ARGUMENTS
377C
378 INTEGER ORDERS(NSO)
379C
380C LOCAL VARIABLES
381C
382 INTEGER I,J
383 INTEGER SQSPLITORDERS(NSQSO,NSO)
384 DATA (SQSPLITORDERS( 1,I),I= 1, 1) / 0/
385 COMMON/SQPLITORDERS/SQPLITORDERS
386C
387C BEGIN CODE
388C
389 DO I=1,NSQSO
390 DO J=1,NSO
391 IF (ORDERS(J).NE.SQSPLITORDERS(I,J)) GOTO 1009
392 ENDDO
393 SOINDEX_FOR_SQUARED_ORDERS = I
394 RETURN
395 1009 CONTINUE
396 ENDDO
397
398 WRITE(*,*) 'ERROR:: Stopping in function'
399 WRITE(*,*) 'SOINDEX_FOR_SQUARED_ORDERS'
400 WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
401 STOP
402
403 END
404
405 SUBROUTINE GET_NSQSO_BORN(NSQSO)
406C
407C Simple subroutine returning the number of squared split order
408C contributions returned when calling smatrix_split_orders
409C
410
411 INTEGER NSQUAREDSO
412 PARAMETER (NSQUAREDSO=1)
413
414 INTEGER NSQSO
415
416 NSQSO=NSQUAREDSO
417
418 END
419
420C This is the inverse subroutine of SOINDEX_FOR_SQUARED_ORDERS.
421C Not directly useful, but provided nonetheless.
422 SUBROUTINE GET_SQUARED_ORDERS_FOR_SOINDEX(SOINDEX,ORDERS)
423C
424C This functions returns the orders identified by the squared
425C split order index in argument. Order values correspond to
426C following list of couplings (and in this order):
427C ['QCD']
428C
429C CONSTANTS
430C
431 INTEGER NSO, NSQSO
432 PARAMETER (NSO=1, NSQSO=1)
433C
434C ARGUMENTS
435C
436 INTEGER SOINDEX, ORDERS(NSO)
437C
438C LOCAL VARIABLES
439C
440 INTEGER I
441 INTEGER SQPLITORDERS(NSQSO,NSO)
442 COMMON/SQPLITORDERS/SQPLITORDERS
443C
444C BEGIN CODE
445C
446 IF (SOINDEX.GT.0.AND.SOINDEX.LE.NSQSO) THEN
447 DO I=1,NSO
448 ORDERS(I) = SQPLITORDERS(SOINDEX,I)
449 ENDDO
450 RETURN
451 ENDIF
452
453 WRITE(*,*) 'ERROR:: Stopping function GET_SQUARED_ORDERS_FOR_SOI'
454 $ //'NDEX'
455 WRITE(*,*) 'Could not find squared orders index ',SOINDEX
456 STOP
457
458 END SUBROUTINE
459
460C This is the inverse subroutine of getting amplitude SO orders.
461C Not directly useful, but provided nonetheless.
462 SUBROUTINE GET_ORDERS_FOR_AMPSOINDEX(SOINDEX,ORDERS)
463C
464C This functions returns the orders identified by the split order
465C index in argument. Order values correspond to following list of
466C couplings (and in this order):
467C ['QCD']
468C
469C CONSTANTS
470C
471 INTEGER NSO, NAMPSO
472 PARAMETER (NSO=1, NAMPSO=1)
473C
474C ARGUMENTS
475C
476 INTEGER SOINDEX, ORDERS(NSO)
477C
478C LOCAL VARIABLES
479C
480 INTEGER I
481 INTEGER AMPSPLITORDERS(NAMPSO,NSO)
482 COMMON/AMPSPLITORDERS/AMPSPLITORDERS
483C
484C BEGIN CODE
485C
486 IF (SOINDEX.GT.0.AND.SOINDEX.LE.NAMPSO) THEN
487 DO I=1,NSO
488 ORDERS(I) = AMPSPLITORDERS(SOINDEX,I)
489 ENDDO
490 RETURN
491 ENDIF
492
493 WRITE(*,*) 'ERROR:: Stopping function GET_ORDERS_FOR_AMPSOINDEX'
494 WRITE(*,*) 'Could not find amplitude split orders index ',SOINDEX
495 STOP
496
497 END SUBROUTINE
498
499C This function is not directly useful, but included for completene
500C ss
501 INTEGER FUNCTION SOINDEX_FOR_AMPORDERS(ORDERS)
502C
503C This functions returns the integer index identifying the
504C amplitude split orders passed in argument which correspond to
505C the values of the following list of couplings (and in this
506C order):
507C ['QCD']
508C
509C CONSTANTS
510C
511 INTEGER NSO, NAMPSO
512 PARAMETER (NSO=1, NAMPSO=1)
513C
514C ARGUMENTS
515C
516 INTEGER ORDERS(NSO)
517C
518C LOCAL VARIABLES
519C
520 INTEGER I,J
521 INTEGER AMPSPLITORDERS(NAMPSO,NSO)
522 COMMON/AMPSPLITORDERS/AMPSPLITORDERS
523C
524C BEGIN CODE
525C
526 DO I=1,NAMPSO
527 DO J=1,NSO
528 IF (ORDERS(J).NE.AMPSPLITORDERS(I,J)) GOTO 1009
529 ENDDO
530 SOINDEX_FOR_AMPORDERS = I
531 RETURN
532 1009 CONTINUE
533 ENDDO
534
535 WRITE(*,*) 'ERROR:: Stopping function SOINDEX_FOR_AMPORDERS'
536 WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO)
537 STOP
538
539 END
540
0541
=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa.f 1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa.f 2016-02-24 13:54:50 +0000
@@ -0,0 +1,756 @@
1 PROGRAM DRIVER
2C *****************************************************************
3C ********
4C THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT.
5C IT USES A SIMPLE PHASE SPACE GENERATOR
6C *****************************************************************
7C ********
8 IMPLICIT NONE
9C
10C CONSTANTS
11C
12 REAL*8 ZERO
13 PARAMETER (ZERO=0D0)
14
15 LOGICAL READPS
16 PARAMETER (READPS = .FALSE.)
17
18 INTEGER NPSPOINTS
19 PARAMETER (NPSPOINTS = 4)
20
21C integer nexternal and number particles (incoming+outgoing) in
22C the me
23 INTEGER NEXTERNAL, NINCOMING
24 PARAMETER (NEXTERNAL=3,NINCOMING=2)
25
26 CHARACTER(512) MADLOOPRESOURCEPATH
27
28C
29C INCLUDE FILES
30C
31C the include file with the values of the parameters and masses
32C
33 INCLUDE 'coupl.inc'
34C particle masses
35 REAL*8 PMASS(NEXTERNAL)
36C integer n_max_cg
37 INCLUDE 'ngraphs.inc'
38 INCLUDE 'nsqso_born.inc'
39 INCLUDE 'nsquaredSO.inc'
40
41C
42C LOCAL
43C
44 INTEGER I,J,K
45C four momenta. Energy is the zeroth component.
46 REAL*8 P(0:3,NEXTERNAL)
47 INTEGER MATELEM_ARRAY_DIM
48 REAL*8 , ALLOCATABLE :: MATELEM(:,:)
49 REAL*8 SQRTS,AO2PI,TOTMASS
50C sqrt(s)= center of mass energy
51 REAL*8 PIN(0:3), POUT(0:3)
52 CHARACTER*120 BUFF(NEXTERNAL)
53 INTEGER RETURNCODE, UNITS, TENS, HUNDREDS
54 INTEGER NSQUAREDSO_LOOP
55 REAL*8 , ALLOCATABLE :: PREC_FOUND(:)
56
57C
58C GLOBAL VARIABLES
59C
60C This is from ML code for the list of split orders selected by
61C the process definition
62C
63 INTEGER NLOOPCHOSEN
64 CHARACTER*20 CHOSEN_LOOP_SO_INDICES(NSQUAREDSO)
65 LOGICAL CHOSEN_LOOP_SO_CONFIGS(NSQUAREDSO)
66 COMMON/CHOSEN_LOOP_SQSO/CHOSEN_LOOP_SO_CONFIGS
67 INTEGER NBORNCHOSEN
68 CHARACTER*20 CHOSEN_BORN_SO_INDICES(NSQSO_BORN)
69 LOGICAL CHOSEN_BORN_SO_CONFIGS(NSQSO_BORN)
70 COMMON/CHOSEN_BORN_SQSO/CHOSEN_BORN_SO_CONFIGS
71
72C
73C SAVED VARIABLES
74C
75 LOGICAL INIT
76 DATA INIT/.TRUE./
77 COMMON/INITCHECKSA/INIT
78C
79C EXTERNAL
80C
81 REAL*8 DOT
82 EXTERNAL DOT
83
84C
85C BEGIN CODE
86C
87C
88
89 IF (INIT) THEN
90 INIT=.FALSE.
91 CALL GET_ANSWER_DIMENSION(MATELEM_ARRAY_DIM)
92 ALLOCATE(MATELEM(0:3,0:MATELEM_ARRAY_DIM))
93 CALL GET_NSQSO_LOOP(NSQUAREDSO_LOOP)
94 ALLOCATE(PREC_FOUND(0:NSQUAREDSO_LOOP))
95
96C INITIALIZATION CALLS
97C
98C Call to initialize the values of the couplings, masses and
99C widths
100C used in the evaluation of the matrix element. The primary
101C parameters of the
102C models are read from Cards/param_card.dat. The secondary
103C parameters are calculated
104C in Source/MODEL/couplings.f. The values are stored in common
105C blocks that are listed
106C in coupl.inc .
107C first call to setup the paramaters
108 CALL SETPARA('param_card.dat')
109C set up masses
110 INCLUDE 'pmass.inc'
111
112 ENDIF
113
114
115C Start by initializing what is the squared split orders indices
116C chosen
117 NLOOPCHOSEN=0
118 DO I=1,NSQUAREDSO
119 IF (CHOSEN_LOOP_SO_CONFIGS(I)) THEN
120 NLOOPCHOSEN=NLOOPCHOSEN+1
121 WRITE(CHOSEN_LOOP_SO_INDICES(NLOOPCHOSEN),'(I3,A2)') I,'L)'
122 ENDIF
123 ENDDO
124 NBORNCHOSEN=0
125 DO I=1,NSQSO_BORN
126 IF (CHOSEN_BORN_SO_CONFIGS(I)) THEN
127 NBORNCHOSEN=NBORNCHOSEN+1
128 WRITE(CHOSEN_BORN_SO_INDICES(NBORNCHOSEN),'(I3,A2)') I,'B)'
129 ENDIF
130 ENDDO
131
132 AO2PI=G**2/(8.D0*(3.14159265358979323846D0**2))
133
134 WRITE(*,*) 'AO2PI=',AO2PI
135C Now use a simple multipurpose PS generator (RAMBO) just to get a
136C RANDOM set of four momenta of given masses pmass(i) to be used
137C to evaluate
138C the madgraph matrix-element.
139C Alternatevely, here the user can call or set the four momenta at
140C his will, see below.
141C
142 IF(NINCOMING.EQ.1) THEN
143 SQRTS=PMASS(1)
144 ELSE
145 TOTMASS = 0.0D0
146 DO I=1,NEXTERNAL
147 TOTMASS = TOTMASS + PMASS(I)
148 ENDDO
149C CMS energy in GEV
150 SQRTS=MAX(1000D0,2.0D0*TOTMASS)
151 ENDIF
152
153 CALL PRINTOUT()
154
155
156
157 DO K=1,NPSPOINTS
158
159 IF(READPS) THEN
160 OPEN(967, FILE='PS.input', ERR=976, STATUS='OLD', ACTION='RE'
161 $ //'AD')
162 DO I=1,NEXTERNAL
163 READ(967,*,END=978) P(0,I),P(1,I),P(2,I),P(3,I)
164 ENDDO
165 GOTO 978
166 976 CONTINUE
167 STOP 'Could not read the PS.input phase-space point.'
168 978 CONTINUE
169 CLOSE(967)
170 ELSE
171 IF ((NINCOMING.EQ.2).AND.((NEXTERNAL - NINCOMING .EQ.1))
172 $ ) THEN
173 IF (PMASS(3).EQ.0.0D0) THEN
174 STOP 'Cannot generate 2>1 kin. config. with m3=0.0d0'
175 ELSE
176C deal with the case of only one particle in the final
177C state
178 P(0,1) = PMASS(3)/2D0
179 P(1,1) = 0D0
180 P(2,1) = 0D0
181 P(3,1) = PMASS(3)/2D0
182 IF (PMASS(1).GT.0D0) THEN
183 P(3,1) = DSQRT(PMASS(3)**2/4D0 - PMASS(1)**2)
184 ENDIF
185 P(0,2) = PMASS(3)/2D0
186 P(1,2) = 0D0
187 P(2,2) = 0D0
188 P(3,2) = -PMASS(3)/2D0
189 IF (PMASS(2) > 0D0) THEN
190 P(3,2) = -DSQRT(PMASS(3)**2/4D0 - PMASS(1)**2)
191 ENDIF
192 P(0,3) = PMASS(3)
193 P(1,3) = 0D0
194 P(2,3) = 0D0
195 P(3,3) = 0D0
196 ENDIF
197 ELSE
198 CALL GET_MOMENTA(SQRTS,PMASS,P)
199 ENDIF
200 ENDIF
201
202 DO I=0,3
203 PIN(I)=0.0D0
204 DO J=1,NINCOMING
205 PIN(I)=PIN(I)+P(I,J)
206 ENDDO
207 ENDDO
208
209C In standalone mode, always use sqrt_s as the renormalization
210C scale.
211 SQRTS=DSQRT(DABS(DOT(PIN(0),PIN(0))))
212 MU_R=SQRTS
213
214C Update the couplings with the new MU_R
215 CALL UPDATE_AS_PARAM()
216
217C Optionally the user can set where to find the MadLoop5_resource
218C s folder.
219C Otherwise it will look for it automatically and find it if it
220C has not
221C been moved
222C MadLoopResourcePath = '<MadLoop5_resources_path>'
223C CALL SETMADLOOPPATH(MadLoopResourcePath)
224C To force the stabiliy check to also be performed in the
225C initialization phase
226C CALL FORCE_STABILITY_CHECK(.TRUE.)
227C To chose a particular tartget split order, SOTARGET is an
228C integer labeling
229C the possible squared order couplings contributions (only in
230C optimized mode)
231C CALL SET_COUPLINGORDERS_TARGET(SOTARGET)
232
233
234C
235C Now we can call the matrix element
236C
237 CALL SLOOPMATRIX_THRES(P,MATELEM,-1.0D0,PREC_FOUND,RETURNCODE)
238
239C
240C write the information on the four momenta
241C
242 IF (K.EQ.NPSPOINTS) THEN
243 WRITE (*,*)
244 WRITE (*,*) ' Phase space point:'
245 WRITE (*,*)
246 WRITE (*,*) '---------------------------------'
247 WRITE (*,*) 'n E px py pz m'
248 DO I=1,NEXTERNAL
249 WRITE (*,'(i2,1x,5e15.7)') I, P(0,I),P(1,I),P(2,I),P(3,I)
250 $ ,DSQRT(DABS(DOT(P(0,I),P(0,I))))
251 ENDDO
252 WRITE (*,*) '---------------------------------'
253 WRITE (*,*) 'Detailed result for each coupling order'
254 $ //'s combination.'
255 WRITE(*,*) 'All Born contributions are of split order'
256 $ //'s (QCD=0)'
257 WRITE (*,*) '---------------------------------'
258 WRITE(*,*) 'All loop contributions are of split order'
259 $ //'s (QCD=2)'
260 WRITE (*,*) '---------------------------------'
261 UNITS=MOD(RETURNCODE,10)
262 TENS=(MOD(RETURNCODE,100)-UNITS)/10
263 HUNDREDS=(RETURNCODE-TENS*10-UNITS)/100
264 IF (HUNDREDS.EQ.1) THEN
265 IF (TENS.EQ.3.OR.TENS.EQ.4) THEN
266 WRITE(*,*) 'Unknown numerical stability because MadLoo'
267 $ //'p is in the initialization stage.'
268 ELSE
269 WRITE(*,*) 'Unknown numerical stability, check CTModeRu'
270 $ //'n value in MadLoopParams.dat.'
271 ENDIF
272 ELSEIF (HUNDREDS.EQ.2) THEN
273 WRITE(*,*) 'Stable kinematic configuration (SPS).'
274 ELSEIF (HUNDREDS.EQ.3) THEN
275 WRITE(*,*) 'Unstable kinematic configuration (UPS).'
276 WRITE(*,*) 'Quadruple precision rescue successful.'
277 ELSEIF (HUNDREDS.EQ.4) THEN
278 WRITE(*,*) 'Exceptional kinematic configuration (EPS).'
279 WRITE(*,*) 'Both double an quadruple precision computation'
280 $ //'s, are unstable.'
281 ENDIF
282 IF (TENS.EQ.2.OR.TENS.EQ.4) THEN
283 WRITE(*,*) 'Quadruple precision computation used.'
284 ENDIF
285 IF (HUNDREDS.NE.1) THEN
286 IF (PREC_FOUND(0).GT.0.0D0) THEN
287 WRITE(*,'(1x,a23,1x,1e10.2)') 'Relative accuracy ='
288 $ ,PREC_FOUND(0)
289 ELSEIF (PREC_FOUND(0).EQ.0.0D0) THEN
290 WRITE(*,'(1x,a23,1x,1e10.2,1x,a30)') 'Relative accuracy'
291 $ //' =',PREC_FOUND(0),'(i.e. beyond double precisio'
292 $ //'n)'
293 ELSE
294 WRITE(*,*) 'Estimated accuracy could not be computed fo'
295 $ //'r an unknown reason.'
296 ENDIF
297 ENDIF
298 WRITE (*,*) '---------------------------------'
299 IF (NBORNCHOSEN.EQ.0) THEN
300 WRITE (*,*) 'No Born contribution satisfied the square'
301 $ //'d order constraints.'
302 ELSE IF (NBORNCHOSEN.NE.NSQSO_BORN) THEN
303 WRITE (*,*) 'Selected squared coupling orders combinatio'
304 $ //'n for the Born summed result below:'
305 WRITE (*,*) (CHOSEN_BORN_SO_INDICES(I),I=1,NBORNCHOSEN)
306 ENDIF
307 IF (NLOOPCHOSEN.NE.NSQUAREDSO) THEN
308 WRITE (*,*) 'Selected squared coupling orders combinatio'
309 $ //'n for the loop summed result below:'
310 WRITE (*,*) (CHOSEN_LOOP_SO_INDICES(I),I=1,NLOOPCHOSEN)
311 ENDIF
312 WRITE (*,*) '---------------------------------'
313 WRITE (*,*) 'Matrix element born = ', MATELEM(0,0)
314 $ , ' GeV^',-(2*NEXTERNAL-8)
315 WRITE (*,*) 'Matrix element finite = ', MATELEM(1,0)
316 $ , ' GeV^',-(2*NEXTERNAL-8)
317 WRITE (*,*) 'Matrix element 1eps = ', MATELEM(2,0)
318 $ , ' GeV^',-(2*NEXTERNAL-8)
319 WRITE (*,*) 'Matrix element 2eps = ', MATELEM(3,0)
320 $ , ' GeV^',-(2*NEXTERNAL-8)
321 WRITE (*,*) '---------------------------------'
322 IF (MATELEM(0,0).NE.0.0D0) THEN
323 WRITE (*,*) 'finite / (born*ao2pi) = ', MATELEM(1,0)
324 $ /MATELEM(0,0)/AO2PI
325 WRITE (*,*) '1eps / (born*ao2pi) = ', MATELEM(2,0)
326 $ /MATELEM(0,0)/AO2PI
327 WRITE (*,*) '2eps / (born*ao2pi) = ', MATELEM(3,0)
328 $ /MATELEM(0,0)/AO2PI
329 ELSE
330 WRITE (*,*) 'finite / ao2pi = ', MATELEM(1,0)/AO2PI
331 WRITE (*,*) '1eps / ao2pi = ', MATELEM(2,0)/AO2PI
332 WRITE (*,*) '2eps / ao2pi = ', MATELEM(3,0)/AO2PI
333 ENDIF
334 WRITE (*,*) '---------------------------------'
335
336 OPEN(69, FILE='result.dat', ERR=976, ACTION='WRITE')
337 DO I=1,NEXTERNAL
338 WRITE (69,'(a2,1x,5e25.15)') 'PS',P(0,I),P(1,I),P(2,I),P(3
339 $ ,I)
340 ENDDO
341 WRITE (69,'(a3,1x,i2)') 'EXP',-(2*NEXTERNAL-8)
342 WRITE (69,'(a4,1x,1e25.15)') 'BORN',MATELEM(0,0)
343 IF (MATELEM(0,0).NE.0.0D0) THEN
344 WRITE (69,'(a3,1x,1e25.15)') 'FIN',MATELEM(1,0)/MATELEM(0
345 $ ,0)/AO2PI
346 WRITE (69,'(a4,1x,1e25.15)') '1EPS',MATELEM(2,0)/MATELEM(0
347 $ ,0)/AO2PI
348 WRITE (69,'(a4,1x,1e25.15)') '2EPS',MATELEM(3,0)/MATELEM(0
349 $ ,0)/AO2PI
350 ELSE
351 WRITE (69,'(a3,1x,1e25.15)') 'FIN',MATELEM(1,0)/AO2PI
352 WRITE (69,'(a4,1x,1e25.15)') '1EPS',MATELEM(2,0)/AO2PI
353 WRITE (69,'(a4,1x,1e25.15)') '2EPS',MATELEM(3,0)/AO2PI
354 ENDIF
355 WRITE (69,'(a6,1x,1e25.15)') 'ASO2PI',AO2PI
356 WRITE (69,*) 'Export_Format Default'
357 WRITE (69,'(a7,1x,i3)') 'RETCODE',RETURNCODE
358 WRITE (69,'(a3,1x,1e10.4)') 'ACC',PREC_FOUND(0)
359 WRITE (69,*) 'Born_kept',(CHOSEN_BORN_SO_CONFIGS(I),I=1
360 $ ,NSQSO_BORN)
361 WRITE (69,*) 'Loop_kept',(CHOSEN_LOOP_SO_CONFIGS(I),I=1
362 $ ,NSQUAREDSO)
363 WRITE (69,*) 'Born_SO_Results 0'
364 WRITE (69,*) 'SO_Born BORN ',MATELEM(0,1)
365 WRITE (69,*) 'Split_Orders_Names QCD'
366 WRITE (69,*) 'Loop_SO_Results 2'
367 WRITE (69,*) 'SO_Loop ACC ',PREC_FOUND(1)
368 WRITE (69,*) 'SO_Loop FIN ',MATELEM(1,1)
369 WRITE (69,*) 'SO_Loop 1EPS ',MATELEM(2,1)
370 WRITE (69,*) 'SO_Loop 2EPS ',MATELEM(3,1)
371 CLOSE(69)
372 ELSE
373 WRITE (*,*) 'PS Point #',K,' done.'
374 ENDIF
375 ENDDO
376
377C C
378C C Copy down here (or read in) the four momenta as a string.
379C C
380C C
381C buff(1)=" 1 0.5630480E+04 0.0000000E+00 0.0000000E+00
382C 0.5630480E+04"
383C buff(2)=" 2 0.5630480E+04 0.0000000E+00 0.0000000E+00
384C -0.5630480E+04"
385C buff(3)=" 3 0.5466073E+04 0.4443190E+03 0.2446331E+04
386C -0.4864732E+04"
387C buff(4)=" 4 0.8785819E+03 -0.2533886E+03 0.2741971E+03
388C 0.7759741E+03"
389C buff(5)=" 5 0.4916306E+04 -0.1909305E+03 -0.2720528E+04
390C 0.4088757E+04"
391C C
392C C Here the k,E,px,py,pz are read from the string into the
393C momenta array.
394C C k=1,2 : incoming
395C C k=3,nexternal : outgoing
396C C
397C do i=1,nexternal
398C read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i)
399C enddo
400C
401C C print the momenta out
402C
403C do i=1,nexternal
404C write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i),
405C &dsqrt(dabs(DOT(p(0,i),p(0,i))))
406C enddo
407C
408C CALL SLOOPMATRIX(P,MATELEM)
409C
410C write (*,*) "-------------------------------------------------"
411C write (*,*) "Matrix element = ", MATELEM(1), " GeV^",-(2*nexterna
412C l-8)
413C write (*,*) "-------------------------------------------------"
414
415 DEALLOCATE(MATELEM)
416 DEALLOCATE(PREC_FOUND)
417
418 END
419
420
421
422
423 DOUBLE PRECISION FUNCTION DOT(P1,P2)
424C *************************************************************
425C 4-Vector Dot product
426C *************************************************************
427 IMPLICIT NONE
428 DOUBLE PRECISION P1(0:3),P2(0:3)
429 DOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
430 END
431
432
433 SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P)
434C auxiliary function to change convention between madgraph and
435C rambo
436C four momenta.
437 IMPLICIT NONE
438 INTEGER NEXTERNAL, NINCOMING
439 PARAMETER (NEXTERNAL=3,NINCOMING=2)
440C ARGUMENTS
441 REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT
442C LOCAL
443 INTEGER I
444 REAL*8 ETOT2,MOM,M1,M2,E1,E2
445
446 ETOT2=ENERGY**2
447 M1=PMASS(1)
448 M2=PMASS(2)
449 MOM=(ETOT2**2 - 2*ETOT2*M1**2 + M1**4 - 2*ETOT2*M2**2 - 2*M1**2
450 $ *M2**2 + M2**4)/(4.*ETOT2)
451 MOM=DSQRT(MOM)
452 E1=DSQRT(MOM**2+M1**2)
453 E2=DSQRT(MOM**2+M2**2)
454C write (*,*) e1+e2,mom
455
456 IF(NINCOMING.EQ.2) THEN
457
458 P(0,1)=E1
459 P(1,1)=0D0
460 P(2,1)=0D0
461 P(3,1)=MOM
462
463 P(0,2)=E2
464 P(1,2)=0D0
465 P(2,2)=0D0
466 P(3,2)=-MOM
467
468 CALL RAMBO(NEXTERNAL-2,ENERGY,PMASS(3),PRAMBO,WGT)
469 DO I=3, NEXTERNAL
470 P(0,I)=PRAMBO(4,I-2)
471 P(1,I)=PRAMBO(1,I-2)
472 P(2,I)=PRAMBO(2,I-2)
473 P(3,I)=PRAMBO(3,I-2)
474 ENDDO
475
476 ELSEIF(NINCOMING.EQ.1) THEN
477
478 P(0,1)=ENERGY
479 P(1,1)=0D0
480 P(2,1)=0D0
481 P(3,1)=0D0
482
483 CALL RAMBO(NEXTERNAL-1,ENERGY,PMASS(2),PRAMBO,WGT)
484 DO I=2, NEXTERNAL
485 P(0,I)=PRAMBO(4,I-1)
486 P(1,I)=PRAMBO(1,I-1)
487 P(2,I)=PRAMBO(2,I-1)
488 P(3,I)=PRAMBO(3,I-1)
489 ENDDO
490 ENDIF
491
492 RETURN
493 END
494
495
496 SUBROUTINE RAMBO(N,ET,XM,P,WT)
497C *****************************************************************
498C *****
499C RAMBO *
500C RA(NDOM) M(OMENTA) B(EAUTIFULLY) O(RGANIZED)
501C *
502C *
503C A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR
504C *
505C AUTHORS: S.D. ELLIS, R. KLEISS, W.J. STIRLING
506C *
507C THIS IS VERSION 1.0 - WRITTEN BY R. KLEISS
508C *
509C -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90)
510C *
511C *
512C N = NUMBER OF PARTICLES
513C *
514C ET = TOTAL CENTRE-OF-MASS ENERGY
515C *
516C XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming )
517C *
518C P = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) )
519C *
520C WT = WEIGHT OF THE EVENT
521C *
522C *****************************************************************
523C *****
524 IMPLICIT REAL*8(A-H,O-Z)
525 INTEGER NEXTERNAL, NINCOMING
526 PARAMETER (NEXTERNAL=3,NINCOMING=2)
527 DIMENSION XM(NEXTERNAL-NINCOMING),P(4,NEXTERNAL-NINCOMING)
528 DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4)
529 $ ,B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING)
530 $ ,E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5)
531 SAVE ACC,ITMAX,IBEGIN,IWARN
532 DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
533C
534C INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
535 IF(IBEGIN.NE.0) GOTO 103
536 IBEGIN=1
537 TWOPI=8.*DATAN(1.D0)
538 PO2LOG=LOG(TWOPI/4.)
539 Z(2)=PO2LOG
540 DO 101 K=3,(NEXTERNAL-NINCOMING)
541 101 Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
542 DO 102 K=3,(NEXTERNAL-NINCOMING)
543 102 Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
544C
545C CHECK ON THE NUMBER OF PARTICLES
546 103 IF(N.GT.1.AND.N.LT.101) GOTO 104
547 PRINT 1001,N
548 STOP
549C
550C CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
551 104 XMT=0.
552 NM=0
553 DO 105 I=1,N
554 IF(XM(I).NE.0.D0) NM=NM+1
555 105 XMT=XMT+ABS(XM(I))
556 IF(XMT.LE.ET) GOTO 201
557 PRINT 1002,XMT,ET
558 STOP
559C
560C THE PARAMETER VALUES ARE NOW ACCEPTED
561C
562C GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
563 201 DO 202 I=1,N
564 R1=RN(1)
565 C=2.*R1-1.
566 S=SQRT(1.-C*C)
567 F=TWOPI*RN(2)
568 R1=RN(3)
569 R2=RN(4)
570 Q(4,I)=-LOG(R1*R2)
571 Q(3,I)=Q(4,I)*C
572 Q(2,I)=Q(4,I)*S*COS(F)
573 202 Q(1,I)=Q(4,I)*S*SIN(F)
574C
575C CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
576 DO 203 I=1,4
577 203 R(I)=0.
578 DO 204 I=1,N
579 DO 204 K=1,4
580 204 R(K)=R(K)+Q(K,I)
581 RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
582 DO 205 K=1,3
583 205 B(K)=-R(K)/RMAS
584 G=R(4)/RMAS
585 A=1./(1.+G)
586 X=ET/RMAS
587C
588C TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
589 DO 207 I=1,N
590 BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
591 DO 206 K=1,3
592 206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
593 207 P(4,I)=X*(G*Q(4,I)+BQ)
594C
595C CALCULATE WEIGHT AND POSSIBLE WARNINGS
596 WT=PO2LOG
597 IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
598 IF(WT.GE.-180.D0) GOTO 208
599 IF(IWARN(1).LE.5) PRINT 1004,WT
600 IWARN(1)=IWARN(1)+1
601 208 IF(WT.LE. 174.D0) GOTO 209
602 IF(IWARN(2).LE.5) PRINT 1005,WT
603 IWARN(2)=IWARN(2)+1
604C
605C RETURN FOR WEIGHTED MASSLESS MOMENTA
606 209 IF(NM.NE.0) GOTO 210
607C RETURN LOG OF WEIGHT
608 WT=WT
609 RETURN
610C
611C MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
612 210 XMAX=SQRT(1.-(XMT/ET)**2)
613 DO 301 I=1,N
614 XM2(I)=XM(I)**2
615 301 P2(I)=P(4,I)**2
616 ITER=0
617 X=XMAX
618 ACCU=ET*ACC
619 302 F0=-ET
620 G0=0.
621 X2=X*X
622 DO 303 I=1,N
623 E(I)=SQRT(XM2(I)+X2*P2(I))
624 F0=F0+E(I)
625 303 G0=G0+P2(I)/E(I)
626 IF(ABS(F0).LE.ACCU) GOTO 305
627 ITER=ITER+1
628 IF(ITER.LE.ITMAX) GOTO 304
629 PRINT 1006,ITMAX
630 GOTO 305
631 304 X=X-F0/(X*G0)
632 GOTO 302
633 305 DO 307 I=1,N
634 V(I)=X*P(4,I)
635 DO 306 K=1,3
636 306 P(K,I)=X*P(K,I)
637 307 P(4,I)=E(I)
638C
639C CALCULATE THE MASS-EFFECT WEIGHT FACTOR
640 WT2=1.
641 WT3=0.
642 DO 308 I=1,N
643 WT2=WT2*V(I)/E(I)
644 308 WT3=WT3+V(I)**2/E(I)
645 WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
646C
647C RETURN FOR WEIGHTED MASSIVE MOMENTA
648 WT=WT+WTM
649 IF(WT.GE.-180.D0) GOTO 309
650 IF(IWARN(3).LE.5) PRINT 1004,WT
651 IWARN(3)=IWARN(3)+1
652 309 IF(WT.LE. 174.D0) GOTO 310
653 IF(IWARN(4).LE.5) PRINT 1005,WT
654 IWARN(4)=IWARN(4)+1
655C RETURN LOG OF WEIGHT
656 310 WT=WT
657 RETURN
658C
659 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
660 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',' SMALLE'
661 $ //'R THAN TOTAL ENERGY =',D15.6)
662 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
663 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY OVERFLOW')
664 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE'
665 $ ,' DESIRED ACCURACY =',D15.6)
666 END
667
668 FUNCTION RN(IDUMMY)
669 REAL*8 RN,RAN
670 SAVE INIT
671 DATA INIT /1/
672 IF (INIT.EQ.1) THEN
673 INIT=0
674 CALL RMARIN(1802,9373)
675 END IF
676C
677 10 CALL RANMAR(RAN)
678 IF (RAN.LT.1D-16) GOTO 10
679 RN=RAN
680C
681 END
682
683
684
685 SUBROUTINE RANMAR(RVEC)
686C -----------------
687C Universal random number generator proposed by Marsaglia and
688C Zaman
689C in report FSU-SCRI-87-50
690C In this version RVEC is a double precision variable.
691 IMPLICIT REAL*8(A-H,O-Z)
692 COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
693 COMMON/ RASET2 / IRANMR,JRANMR
694 SAVE /RASET1/,/RASET2/
695 UNI = RANU(IRANMR) - RANU(JRANMR)
696 IF(UNI .LT. 0D0) UNI = UNI + 1D0
697 RANU(IRANMR) = UNI
698 IRANMR = IRANMR - 1
699 JRANMR = JRANMR - 1
700 IF(IRANMR .EQ. 0) IRANMR = 97
701 IF(JRANMR .EQ. 0) JRANMR = 97
702 RANC = RANC - RANCD
703 IF(RANC .LT. 0D0) RANC = RANC + RANCM
704 UNI = UNI - RANC
705 IF(UNI .LT. 0D0) UNI = UNI + 1D0
706 RVEC = UNI
707 END
708
709 SUBROUTINE RMARIN(IJ,KL)
710C -----------------
711C Initializing routine for RANMAR, must be called before
712C generating
713C any pseudorandom numbers with RANMAR. The input values should
714C be in
715C the ranges 0<=ij<=31328 ; 0<=kl<=30081
716 IMPLICIT REAL*8(A-H,O-Z)
717 COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
718 COMMON/ RASET2 / IRANMR,JRANMR
719 SAVE /RASET1/,/RASET2/
720C This shows correspondence between the simplified input seeds
721C IJ, KL
722C and the original Marsaglia-Zaman seeds I,J,K,L.
723C To get the standard values in the Marsaglia-Zaman paper
724C (i=12,j=34
725C k=56,l=78) put ij=1802, kl=9373
726 I = MOD( IJ/177 , 177 ) + 2
727 J = MOD( IJ , 177 ) + 2
728 K = MOD( KL/169 , 178 ) + 1
729 L = MOD( KL , 169 )
730 DO 300 II = 1 , 97
731 S = 0D0
732 T = .5D0
733 DO 200 JJ = 1 , 24
734 M = MOD( MOD(I*J,179)*K , 179 )
735 I = J
736 J = K
737 K = M
738 L = MOD( 53*L+1 , 169 )
739 IF(MOD(L*M,64) .GE. 32) S = S + T
740 T = .5D0*T
741 200 CONTINUE
742 RANU(II) = S
743 300 CONTINUE
744 RANC = 362436D0 / 16777216D0
745 RANCD = 7654321D0 / 16777216D0
746 RANCM = 16777213D0 / 16777216D0
747 IRANMR = 97
748 JRANMR = 33
749 END
750
751
752
753
754
755
756
0757
=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f 1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%check_sa_born_splitOrders.f 2016-02-24 13:54:50 +0000
@@ -0,0 +1,524 @@
1 PROGRAM DRIVER
2C *****************************************************************
3C *********
4C THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT.
5C IT USES A SIMPLE PHASE SPACE GENERATOR
6C Fabio Maltoni - 3rd Febraury 2007
7C *****************************************************************
8C *********
9 IMPLICIT NONE
10C
11C CONSTANTS
12C
13 REAL*8 ZERO
14 PARAMETER (ZERO=0D0)
15 INTEGER NSPLITORDERS
16 PARAMETER (NSPLITORDERS=1)
17C
18C INCLUDE FILES
19C
20C --- the include file with the values of the parameters and
21C masses
22 INCLUDE 'coupl.inc'
23C integer nexternal and number particles (incoming+outgoing) in
24C the me
25 INTEGER NEXTERNAL, NINCOMING
26 PARAMETER (NEXTERNAL=3,NINCOMING=2)
27C --- particle masses
28 REAL*8 PMASS(NEXTERNAL)
29 REAL*8 TOTALMASS
30C --- integer n_max_cg
31 INCLUDE 'ngraphs.inc' !how many diagrams (could be useful to know...)
32
33C
34C LOCAL
35C
36 INTEGER I,J,K
37 REAL*8 P(0:3,NEXTERNAL) ! four momenta. Energy is the zeroth component.
38 REAL*8 SQRTS ! sqrt(s)= center of mass energy
39 REAL*8 MATELEM, MATELEMS(0:NSPLITORDERS)
40 REAL*8 PIN(0:3), POUT(0:3)
41 CHARACTER*120 BUFF(NEXTERNAL)
42
43 INTEGER NCHOSEN
44 CHARACTER*20 CHOSEN_SO_INDICES(NSPLITORDERS)
45
46C
47C EXTERNAL
48C
49 REAL*8 DOT
50 EXTERNAL DOT
51
52 LOGICAL CHOSEN_SO_CONFIGS(NSPLITORDERS)
53 COMMON/CHOSEN_BORN_SQSO/CHOSEN_SO_CONFIGS
54
55C -----
56C BEGIN CODE
57C -----
58C
59C Start by initializing what is the squared split orders indices
60C chosen
61 NCHOSEN=0
62 DO I=1,NSPLITORDERS
63 IF (CHOSEN_SO_CONFIGS(I)) THEN
64 NCHOSEN=NCHOSEN+1
65 WRITE(CHOSEN_SO_INDICES(NCHOSEN),'(I3,A1)') I,')'
66 ENDIF
67 ENDDO
68
69C --- INITIALIZATION CALLS
70C
71C --- Call to initialize the values of the couplings, masses and
72C widths
73C used in the evaluation of the matrix element. The primary
74C parameters of the
75C models are read from Cards/param_card.dat. The secondary
76C parameters are calculated
77C in Source/MODEL/couplings.f. The values are stored in common
78C blocks that are listed
79C in coupl.inc .
80
81 CALL SETPARA('param_card.dat') !first call to setup the paramaters
82 INCLUDE 'pmass.inc' !set up masses
83
84 TOTALMASS = 0.0D0
85 DO I=1,NEXTERNAL
86 TOTALMASS = TOTALMASS + PMASS(I)
87 ENDDO
88
89C --- Now use a simple multipurpose PS generator (RAMBO) just to
90C get a
91C RANDOM set of four momenta of given masses pmass(i) to be used
92C to evaluate
93C the madgraph matrix-element.
94C Alternatevely, here the user can call or set the four momenta at
95C his will, see below.
96C
97 IF(NINCOMING.EQ.1) THEN
98 SQRTS=PMASS(1)
99 ELSE
100 SQRTS=1000D0 !CMS energy in GEV
101 IF (SQRTS.LE.2.0D0*TOTALMASS) THEN
102 SQRTS = 2.0D0*TOTALMASS
103 ENDIF
104 ENDIF
105
106 CALL PRINTOUT()
107
108 CALL GET_MOMENTA(SQRTS,PMASS,P)
109C
110C write the information on the four momenta
111C
112 WRITE (*,*)
113 WRITE (*,*) ' Phase space point:'
114 WRITE (*,*)
115 WRITE (*,*) '-----------------------------'
116 WRITE (*,*) 'n E px py pz m '
117 DO I=1,NEXTERNAL
118 WRITE (*,'(i2,1x,5e15.7)') I, P(0,I),P(1,I),P(2,I),P(3,I)
119 $ ,DSQRT(DABS(DOT(P(0,I),P(0,I))))
120 ENDDO
121 WRITE (*,*) '-----------------------------'
122
123C
124C Now we can call the matrix element!
125C
126 CALL SMATRIX_SPLITORDERS(P,MATELEMS)
127 MATELEM=MATELEMS(0)
128 WRITE(*,*) '1) Matrix element for (QCD=0) = ',MATELEMS(1)
129C
130 IF (NCHOSEN.NE.NSPLITORDERS) THEN
131 WRITE (*,*) 'Selected squared coupling orders combination fo'
132 $ //'r the sum below:'
133 WRITE (*,*) (CHOSEN_SO_INDICES(I),I=1,NCHOSEN)
134 ENDIF
135 WRITE (*,*) 'Total Matrix element = ', MATELEM, ' GeV^',
136 $ -(2*NEXTERNAL-8)
137 WRITE (*,*) '-----------------------------'
138
139C c
140C c Copy down here (or read in) the four momenta as a string.
141C c
142C c
143C buff(1)=" 1 0.5630480E+04 0.0000000E+00 0.0000000E+00
144C 0.5630480E+04"
145C buff(2)=" 2 0.5630480E+04 0.0000000E+00 0.0000000E+00
146C -0.5630480E+04"
147C buff(3)=" 3 0.5466073E+04 0.4443190E+03 0.2446331E+04
148C -0.4864732E+04"
149C buff(4)=" 4 0.8785819E+03 -0.2533886E+03 0.2741971E+03
150C 0.7759741E+03"
151C buff(5)=" 5 0.4916306E+04 -0.1909305E+03 -0.2720528E+04
152C 0.4088757E+04"
153C c
154C c Here the k,E,px,py,pz are read from the string into the
155C momenta array.
156C c k=1,2 : incoming
157C c k=3,nexternal : outgoing
158C c
159C do i=1,nexternal
160C read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i)
161C enddo
162C
163C c- print the momenta out
164C
165C do i=1,nexternal
166C write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i),
167C .dsqrt(dabs(DOT(p(0,i),p(0,i))))
168C enddo
169C
170C CALL SMATRIX(P,MATELEM)
171C
172C write (*,*) "-------------------------------------------------"
173C write (*,*) "Matrix element = ", MATELEM, " GeV^",-(2*nexternal-8
174C )
175C write (*,*) "-------------------------------------------------"
176
177 END
178
179
180
181
182 DOUBLE PRECISION FUNCTION DOT(P1,P2)
183C *****************************************************************
184C ***********
185C 4-Vector Dot product
186C *****************************************************************
187C ***********
188 IMPLICIT NONE
189 DOUBLE PRECISION P1(0:3),P2(0:3)
190 DOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3)
191 END
192
193
194 SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P)
195C ---- auxiliary function to change convention between madgraph
196C and rambo
197C ---- four momenta.
198 IMPLICIT NONE
199C integer nexternal and number particles (incoming+outgoing) in
200C the me
201 INTEGER NEXTERNAL, NINCOMING
202 PARAMETER (NEXTERNAL=3,NINCOMING=2)
203C ARGUMENTS
204 REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT
205C LOCAL
206 INTEGER I
207 REAL*8 ETOT2,MOM,M1,M2,E1,E2
208
209 ETOT2=ENERGY**2
210 M1=PMASS(1)
211 M2=PMASS(2)
212 MOM=(ETOT2**2 - 2*ETOT2*M1**2 + M1**4 - 2*ETOT2*M2**2 - 2*M1**2
213 $ *M2**2 + M2**4)/(4.*ETOT2)
214 MOM=DSQRT(MOM)
215 E1=DSQRT(MOM**2+M1**2)
216 E2=DSQRT(MOM**2+M2**2)
217 WRITE (*,*) E1+E2,MOM
218
219 IF(NINCOMING.EQ.2) THEN
220
221 P(0,1)=E1
222 P(1,1)=0D0
223 P(2,1)=0D0
224 P(3,1)=MOM
225
226 P(0,2)=E2
227 P(1,2)=0D0
228 P(2,2)=0D0
229 P(3,2)=-MOM
230
231 CALL RAMBO(NEXTERNAL-2,ENERGY,PMASS(3),PRAMBO,WGT)
232 DO I=3, NEXTERNAL
233 P(0,I)=PRAMBO(4,I-2)
234 P(1,I)=PRAMBO(1,I-2)
235 P(2,I)=PRAMBO(2,I-2)
236 P(3,I)=PRAMBO(3,I-2)
237 ENDDO
238
239 ELSEIF(NINCOMING.EQ.1) THEN
240
241 P(0,1)=ENERGY
242 P(1,1)=0D0
243 P(2,1)=0D0
244 P(3,1)=0D0
245
246 CALL RAMBO(NEXTERNAL-1,ENERGY,PMASS(2),PRAMBO,WGT)
247 DO I=2, NEXTERNAL
248 P(0,I)=PRAMBO(4,I-1)
249 P(1,I)=PRAMBO(1,I-1)
250 P(2,I)=PRAMBO(2,I-1)
251 P(3,I)=PRAMBO(3,I-1)
252 ENDDO
253 ENDIF
254
255 RETURN
256 END
257
258
259 SUBROUTINE RAMBO(N,ET,XM,P,WT)
260C *****************************************************************
261C ******
262C * RAMBO
263C *
264C * RA(NDOM) M(OMENTA) B(EAUTIFULLY) O(RGANIZED)
265C *
266C *
267C *
268C * A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR
269C *
270C * AUTHORS: S.D. ELLIS, R. KLEISS, W.J. STIRLING
271C *
272C * THIS IS VERSION 1.0 - WRITTEN BY R. KLEISS
273C *
274C * -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC
275C (20-08-90) *
276C *
277C *
278C * N = NUMBER OF PARTICLES
279C *
280C * ET = TOTAL CENTRE-OF-MASS ENERGY
281C *
282C * XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming )
283C *
284C * P = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) )
285C *
286C * WT = WEIGHT OF THE EVENT
287C *
288C *****************************************************************
289C ******
290 IMPLICIT REAL*8(A-H,O-Z)
291C integer nexternal and number particles (incoming+outgoing) in
292C the me
293 INTEGER NEXTERNAL, NINCOMING
294 PARAMETER (NEXTERNAL=3,NINCOMING=2)
295 DIMENSION XM(NEXTERNAL-NINCOMING),P(4,NEXTERNAL-NINCOMING)
296 DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4)
297 $ , B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING)
298 $ , E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5)
299 SAVE ACC,ITMAX,IBEGIN,IWARN
300 DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
301C *
302C * INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
303 IF(IBEGIN.NE.0) GOTO 103
304 IBEGIN=1
305 TWOPI=8.*DATAN(1.D0)
306 PO2LOG=LOG(TWOPI/4.)
307 Z(2)=PO2LOG
308 DO 101 K=3,NEXTERNAL-NINCOMING-1
309 101 Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
310 DO 102 K=3,NEXTERNAL-NINCOMING-1
311 102 Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
312C *
313C * CHECK ON THE NUMBER OF PARTICLES
314 103 IF(N.GT.1.AND.N.LT.101) GOTO 104
315 PRINT 1001,N
316 STOP
317C *
318C * CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
319 104 XMT=0.
320 NM=0
321 DO 105 I=1,N
322 IF(XM(I).NE.0.D0) NM=NM+1
323 105 XMT=XMT+ABS(XM(I))
324 IF(XMT.LE.ET) GOTO 201
325 PRINT 1002,XMT,ET
326 STOP
327C *
328C * THE PARAMETER VALUES ARE NOW ACCEPTED
329C *
330C * GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
331 201 DO 202 I=1,N
332 R1=RN(1)
333 C=2.*R1-1.
334 S=SQRT(1.-C*C)
335 F=TWOPI*RN(2)
336 R1=RN(3)
337 R2=RN(4)
338 Q(4,I)=-LOG(R1*R2)
339 Q(3,I)=Q(4,I)*C
340 Q(2,I)=Q(4,I)*S*COS(F)
341 202 Q(1,I)=Q(4,I)*S*SIN(F)
342C *
343C * CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
344 DO 203 I=1,4
345 203 R(I)=0.
346 DO 204 I=1,N
347 DO 204 K=1,4
348 204 R(K)=R(K)+Q(K,I)
349 RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
350 DO 205 K=1,3
351 205 B(K)=-R(K)/RMAS
352 G=R(4)/RMAS
353 A=1./(1.+G)
354 X=ET/RMAS
355C *
356C * TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
357 DO 207 I=1,N
358 BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
359 DO 206 K=1,3
360 206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
361 207 P(4,I)=X*(G*Q(4,I)+BQ)
362C *
363C * CALCULATE WEIGHT AND POSSIBLE WARNINGS
364 WT=PO2LOG
365 IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
366 IF(WT.GE.-180.D0) GOTO 208
367 IF(IWARN(1).LE.5) PRINT 1004,WT
368 IWARN(1)=IWARN(1)+1
369 208 IF(WT.LE. 174.D0) GOTO 209
370 IF(IWARN(2).LE.5) PRINT 1005,WT
371 IWARN(2)=IWARN(2)+1
372C *
373C * RETURN FOR WEIGHTED MASSLESS MOMENTA
374 209 IF(NM.NE.0) GOTO 210
375C * RETURN LOG OF WEIGHT
376 WT=WT
377 RETURN
378C *
379C * MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
380 210 XMAX=SQRT(1.-(XMT/ET)**2)
381 DO 301 I=1,N
382 XM2(I)=XM(I)**2
383 301 P2(I)=P(4,I)**2
384 ITER=0
385 X=XMAX
386 ACCU=ET*ACC
387 302 F0=-ET
388 G0=0.
389 X2=X*X
390 DO 303 I=1,N
391 E(I)=SQRT(XM2(I)+X2*P2(I))
392 F0=F0+E(I)
393 303 G0=G0+P2(I)/E(I)
394 IF(ABS(F0).LE.ACCU) GOTO 305
395 ITER=ITER+1
396 IF(ITER.LE.ITMAX) GOTO 304
397 PRINT 1006,ITMAX
398 GOTO 305
399 304 X=X-F0/(X*G0)
400 GOTO 302
401 305 DO 307 I=1,N
402 V(I)=X*P(4,I)
403 DO 306 K=1,3
404 306 P(K,I)=X*P(K,I)
405 307 P(4,I)=E(I)
406C *
407C * CALCULATE THE MASS-EFFECT WEIGHT FACTOR
408 WT2=1.
409 WT3=0.
410 DO 308 I=1,N
411 WT2=WT2*V(I)/E(I)
412 308 WT3=WT3+V(I)**2/E(I)
413 WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
414C *
415C * RETURN FOR WEIGHTED MASSIVE MOMENTA
416 WT=WT+WTM
417 IF(WT.GE.-180.D0) GOTO 309
418 IF(IWARN(3).LE.5) PRINT 1004,WT
419 IWARN(3)=IWARN(3)+1
420 309 IF(WT.LE. 174.D0) GOTO 310
421 IF(IWARN(4).LE.5) PRINT 1005,WT
422 IWARN(4)=IWARN(4)+1
423C * RETURN LOG OF WEIGHT
424 310 WT=WT
425 RETURN
426C *
427 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
428 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT SMALLER THA'
429 $ //'N TOTAL ENERGY =',D15.6)
430 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
431 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY OVERFLOW')
432 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE DESIRE'
433 $ //'D ACCURACY =',D15.6)
434 END
435
436 FUNCTION RN(IDUMMY)
437 REAL*8 RN,RAN
438 SAVE INIT
439 DATA INIT /1/
440 IF (INIT.EQ.1) THEN
441 INIT=0
442 CALL RMARIN(1802,9373)
443 END IF
444C *
445 10 CALL RANMAR(RAN)
446 IF (RAN.LT.1D-16) GOTO 10
447 RN=RAN
448C *
449 END
450
451
452
453 SUBROUTINE RANMAR(RVEC)
454C * -----------------
455C * Universal random number generator proposed by Marsaglia and
456C Zaman
457C * in report FSU-SCRI-87-50
458C * In this version RVEC is a double precision variable.
459 IMPLICIT REAL*8(A-H,O-Z)
460 COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
461 COMMON/ RASET2 / IRANMR,JRANMR
462 SAVE /RASET1/,/RASET2/
463 UNI = RANU(IRANMR) - RANU(JRANMR)
464 IF(UNI .LT. 0D0) UNI = UNI + 1D0
465 RANU(IRANMR) = UNI
466 IRANMR = IRANMR - 1
467 JRANMR = JRANMR - 1
468 IF(IRANMR .EQ. 0) IRANMR = 97
469 IF(JRANMR .EQ. 0) JRANMR = 97
470 RANC = RANC - RANCD
471 IF(RANC .LT. 0D0) RANC = RANC + RANCM
472 UNI = UNI - RANC
473 IF(UNI .LT. 0D0) UNI = UNI + 1D0
474 RVEC = UNI
475 END
476
477 SUBROUTINE RMARIN(IJ,KL)
478C * -----------------
479C * Initializing routine for RANMAR, must be called before
480C generating
481C * any pseudorandom numbers with RANMAR. The input values
482C should be in
483C * the ranges 0<=ij<=31328 ; 0<=kl<=30081
484 IMPLICIT REAL*8(A-H,O-Z)
485 COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
486 COMMON/ RASET2 / IRANMR,JRANMR
487 SAVE /RASET1/,/RASET2/
488C * This shows correspondence between the simplified input seeds
489C IJ, KL
490C * and the original Marsaglia-Zaman seeds I,J,K,L.
491C * To get the standard values in the Marsaglia-Zaman paper
492C (i=12,j=34
493C * k=56,l=78) put ij=1802, kl=9373
494 I = MOD( IJ/177 , 177 ) + 2
495 J = MOD( IJ , 177 ) + 2
496 K = MOD( KL/169 , 178 ) + 1
497 L = MOD( KL , 169 )
498 DO 300 II = 1 , 97
499 S = 0D0
500 T = .5D0
501 DO 200 JJ = 1 , 24
502 M = MOD( MOD(I*J,179)*K , 179 )
503 I = J
504 J = K
505 K = M
506 L = MOD( 53*L+1 , 169 )
507 IF(MOD(L*M,64) .GE. 32) S = S + T
508 T = .5D0*T
509 200 CONTINUE
510 RANU(II) = S
511 300 CONTINUE
512 RANC = 362436D0 / 16777216D0
513 RANCD = 7654321D0 / 16777216D0
514 RANCM = 16777213D0 / 16777216D0
515 IRANMR = 97
516 JRANMR = 33
517 END
518
519
520
521
522
523
524
0525
=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%coef_specs.inc'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%coef_specs.inc 1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%coef_specs.inc 2016-02-24 13:54:50 +0000
@@ -0,0 +1,6 @@
1 INTEGER MAXLWFSIZE
2 PARAMETER (MAXLWFSIZE=4)
3 INTEGER LOOP_MAXCOEFS
4 PARAMETER (LOOP_MAXCOEFS=15)
5 INTEGER VERTEXMAXCOEFS
6 PARAMETER (VERTEXMAXCOEFS=5)
07
=== added file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%improve_ps.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%improve_ps.f 1970-01-01 00:00:00 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%improve_ps.f 2016-02-24 13:54:50 +0000
@@ -0,0 +1,981 @@
1 SUBROUTINE IMPROVE_PS_POINT_PRECISION(P)
2 IMPLICIT NONE
3C
4C CONSTANTS
5C
6 INTEGER NEXTERNAL
7 PARAMETER (NEXTERNAL=3)
8C
9C ARGUMENTS
10C
11 DOUBLE PRECISION P(0:3,NEXTERNAL)
12 REAL*16 QP_P(0:3,NEXTERNAL)
13C
14C LOCAL VARIABLES
15C
16 INTEGER I,J
17
18C ----------
19C BEGIN CODE
20C ----------
21
22 DO I=1,NEXTERNAL
23 DO J=0,3
24 QP_P(J,I)=P(J,I)
25 ENDDO
26 ENDDO
27
28 CALL MP_IMPROVE_PS_POINT_PRECISION(QP_P)
29
30 DO I=1,NEXTERNAL
31 DO J=0,3
32 P(J,I)=QP_P(J,I)
33 ENDDO
34 ENDDO
35
36 END
37
38
39 SUBROUTINE MP_IMPROVE_PS_POINT_PRECISION(P)
40 IMPLICIT NONE
41C
42C CONSTANTS
43C
44 INTEGER NEXTERNAL
45 PARAMETER (NEXTERNAL=3)
46C
47C ARGUMENTS
48C
49 REAL*16 P(0:3,NEXTERNAL)
50C
51C LOCAL VARIABLES
52C
53 INTEGER I,J
54 INTEGER ERRCODE,ERRCODETMP
55 REAL*16 NEWP(0:3,NEXTERNAL)
56C
57C FUNCTIONS
58C
59 LOGICAL MP_IS_PHYSICAL
60C
61C SAVED VARIABLES
62C
63 INCLUDE 'MadLoopParams.inc'
64C
65C SAVED VARIABLES
66C
67 INTEGER WARNED
68 DATA WARNED/0/
69
70 LOGICAL TOLD_SUPPRESS
71 DATA TOLD_SUPPRESS/.FALSE./
72C ----------
73C BEGIN CODE
74C ----------
75
76C ERROR CODES CONVENTION
77C
78C 1 :: None physical PS point input
79C 100-1000 :: Error in the origianl method for restoring
80C precision
81C 1000-9999 :: Error when restoring precision ala PSMC
82C
83 ERRCODETMP=0
84 ERRCODE=0
85
86 DO J=1,NEXTERNAL
87 DO I=0,3
88 NEWP(I,J)=P(I,J)
89 ENDDO
90 ENDDO
91
92C Check the sanity of the original PS point
93 IF (.NOT.MP_IS_PHYSICAL(NEWP,WARNED)) THEN
94 ERRCODE = 1
95 WRITE(*,*) 'ERROR:: The input PS point is not precise enough.'
96 GOTO 100
97 ENDIF
98
99C Now restore the precision
100 IF (IMPROVEPSPOINT.EQ.1) THEN
101 CALL MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE,WARNED)
102 ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN
103 CALL MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODE,WARNED)
104 ENDIF
105 IF (ERRCODE.NE.0) THEN
106 IF (WARNED.LT.20) THEN
107 WRITE(*,*) 'INFO:: Attempting to rescue the precisio'
108 $ //'n improvement with an alternative method.'
109 WARNED=WARNED+1
110 ENDIF
111 IF (IMPROVEPSPOINT.EQ.1) THEN
112 CALL MP_ORIG_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODETMP
113 $ ,WARNED)
114 ELSEIF((IMPROVEPSPOINT.EQ.2).OR.(IMPROVEPSPOINT.LE.0)) THEN
115 CALL MP_PSMC_IMPROVE_PS_POINT_PRECISION(NEWP,ERRCODETMP
116 $ ,WARNED)
117 ENDIF
118 IF (ERRCODETMP.NE.0) GOTO 100
119 ENDIF
120
121C Report to the user or update the PS point.
122
123 GOTO 101
124 100 CONTINUE
125 IF (WARNED.LT.20) THEN
126 WRITE(*,*) 'WARNING:: This PS point could not be improved'
127 $ //'. Error code = ',ERRCODE,ERRCODETMP
128 CALL MP_WRITE_MOM(P)
129 WARNED = WARNED +1
130 ENDIF
131 GOTO 102
132 101 CONTINUE
133 DO J=1,NEXTERNAL
134 DO I=0,3
135 P(I,J)=NEWP(I,J)
136 ENDDO
137 ENDDO
138 102 CONTINUE
139
140 IF (WARNED.GE.20.AND..NOT.TOLD_SUPPRESS) THEN
141 WRITE(*,*) 'INFO:: Further warnings from the improve_p'
142 $ //'s routine will now be supressed.'
143 TOLD_SUPPRESS=.TRUE.
144 ENDIF
145
146 END
147
148
149 FUNCTION MP_IS_CLOSE(P,NEWP,WARNED)
150 IMPLICIT NONE
151C
152C CONSTANTS
153C
154 INTEGER NEXTERNAL
155 PARAMETER (NEXTERNAL=3)
156 REAL*16 ZERO
157 PARAMETER (ZERO=0.0E+00_16)
158 REAL*16 THRS_CLOSE
159 PARAMETER (THRS_CLOSE=1.0E-02_16)
160C
161C ARGUMENTS
162C
163 REAL*16 P(0:3,NEXTERNAL), NEWP(0:3,NEXTERNAL)
164 LOGICAL MP_IS_CLOSE
165 INTEGER WARNED
166C
167C LOCAL VARIABLES
168C
169 INTEGER I,J
170 REAL*16 REF,REF2
171 DOUBLE PRECISION BUFFDP
172
The diff has been truncated for viewing.

Subscribers

People subscribed via source and target branches

to all changes: