Merge lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny

Proposed by Rikkert Frederix on 2018-10-04
Status: Merged
Merge reported by: Olivier Mattelaer
Merged at revision: not available
Proposed branch: lp:~maddevelopers/mg5amcnlo/3.0.1
Merge into: lp:~maddevelopers/mg5amcnlo/FKS_EW_granny
Diff against target: 48081 lines (+15230/-26000)
241 files modified
MadSpin/decay.py (+12/-3)
MadSpin/interface_madspin.py (+342/-174)
MadSpin/madspin (+19/-21)
MadSpin/src/driver.f (+5/-5)
Template/LO/Cards/run_card.dat (+1/-1)
Template/LO/SubProcesses/genps.f (+7/-18)
Template/LO/SubProcesses/reweight.f (+102/-10)
Template/LO/SubProcesses/unwgt.f (+11/-0)
Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f (+880/-0)
Template/NLO/MCatNLO/Makefile_MadFKS (+1/-1)
Template/NLO/SubProcesses/MC_integer.f (+3/-6)
Template/NLO/SubProcesses/analysis_lhe.f (+31/-15)
Template/NLO/SubProcesses/check_poles.f (+2/-2)
Template/NLO/SubProcesses/combine_plots_FO.sh (+56/-110)
Template/NLO/SubProcesses/driver_mintFO.f (+53/-23)
Template/NLO/SubProcesses/driver_mintMC.f (+8/-14)
Template/NLO/SubProcesses/fks_singular.f (+19/-9)
Template/NLO/SubProcesses/genps_fks.f (+1/-1)
Template/NLO/SubProcesses/makefile (+0/-6)
Template/NLO/SubProcesses/mint-integrator2.f (+4/-4)
Template/NLO/SubProcesses/montecarlocounter.f (+3259/-3231)
Template/NLO/SubProcesses/read40.for (+1/-1)
Template/NLO/SubProcesses/split40.f (+0/-51)
Template/NLO/SubProcesses/symmetry_fks_v3.f (+52/-9)
UpdateNotes.txt (+55/-3)
VERSION (+2/-3)
aloha/aloha_writers.py (+24/-7)
aloha/create_aloha.py (+3/-4)
madgraph/core/base_objects.py (+8/-2)
madgraph/fks/fks_base.py (+7/-2)
madgraph/fks/fks_helas_objects.py (+1/-1)
madgraph/interface/amcatnlo_interface.py (+23/-10)
madgraph/interface/amcatnlo_run_interface.py (+91/-106)
madgraph/interface/common_run_interface.py (+106/-26)
madgraph/interface/extended_cmd.py (+40/-11)
madgraph/interface/launch_ext_program.py (+15/-10)
madgraph/interface/loop_interface.py (+13/-5)
madgraph/interface/madevent_interface.py (+32/-11)
madgraph/interface/madgraph_interface.py (+57/-16)
madgraph/interface/reweight_interface.py (+13/-5)
madgraph/iolibs/export_fks.py (+16/-6)
madgraph/iolibs/export_v4.py (+43/-34)
madgraph/iolibs/file_writers.py (+3/-3)
madgraph/iolibs/files.py (+3/-0)
madgraph/iolibs/template_files/matrix_standalone_splitOrders_v4.inc (+2/-2)
madgraph/iolibs/template_files/matrix_standalone_v4.inc (+1/-1)
madgraph/iolibs/ufo_expression_parsers.py (+240/-12)
madgraph/loop/loop_exporters.py (+15/-1)
madgraph/loop/loop_helas_objects.py (+3/-3)
madgraph/madevent/gen_crossxhtml.py (+18/-2)
madgraph/madevent/sum_html.py (+7/-1)
madgraph/various/banner.py (+44/-18)
madgraph/various/hepmc_parser.py (+362/-0)
madgraph/various/lhe_parser.py (+61/-7)
madgraph/various/misc.py (+4/-3)
madgraph/various/systematics.py (+1/-0)
mg5decay/decay_objects.py (+1/-1)
models/check_param_card.py (+9/-5)
models/import_ufo.py (+380/-37)
models/model_reader.py (+7/-1)
models/write_param_card.py (+2/-2)
tests/IOTests.py (+11/-5)
tests/acceptance_tests/test_cmd_amcatnlo.py (+13/-6)
tests/acceptance_tests/test_cmd_madevent.py (+36/-9)
tests/acceptance_tests/test_cmd_madloop.py (+2/-2)
tests/acceptance_tests/test_cmd_reweight.py (+2/-2)
tests/acceptance_tests/test_madspin.py (+166/-0)
tests/acceptance_tests/test_madweight.py (+6/-6)
tests/acceptance_tests/test_model_equivalence.py (+8/-6)
tests/input_files/DM_pion/CT_couplings.py (+1731/-0)
tests/input_files/DM_pion/CT_vertices.py (+1291/-0)
tests/input_files/DM_pion/__init__.py (+48/-0)
tests/input_files/DM_pion/coupling_orders.py (+25/-0)
tests/input_files/DM_pion/couplings.py (+563/-0)
tests/input_files/DM_pion/decays.py (+113/-0)
tests/input_files/DM_pion/function_library.py (+71/-0)
tests/input_files/DM_pion/lorentz.py (+182/-0)
tests/input_files/DM_pion/object_library.py (+377/-0)
tests/input_files/DM_pion/param_pion.dat (+183/-0)
tests/input_files/DM_pion/parameters.py (+858/-0)
tests/input_files/DM_pion/particles.py (+489/-0)
tests/input_files/DM_pion/propagators.py (+35/-0)
tests/input_files/DM_pion/py.py (+67/-0)
tests/input_files/DM_pion/vertices.py (+1007/-0)
tests/input_files/DM_pion/write_param_card.py (+207/-0)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%configs_and_props_info.dat (+0/-1358)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_gg_ttx%leshouche_info.dat (+0/-143)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%configs_and_props_info.dat (+0/-504)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uux_ttx%leshouche_info.dat (+0/-123)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%configs_and_props_info.dat (+0/-504)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksreal/%SubProcesses%P0_uxu_ttx%leshouche_info.dat (+0/-123)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ag_ttx%configs_and_props_info.dat (+0/-709)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ag_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ag_ttx%leshouche_info.dat (+0/-64)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%born.f (+1/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%configs_and_props_info.dat (+0/-1668)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%leshouche_info.dat (+0/-147)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%matrix_3.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%matrix_4.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ddx_ttx%matrix_6.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%born.f (+1/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%configs_and_props_info.dat (+0/-1668)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%leshouche_info.dat (+0/-147)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%matrix_3.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%matrix_4.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_dxd_ttx%matrix_6.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ga_ttx%configs_and_props_info.dat (+0/-709)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ga_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_ga_ttx%leshouche_info.dat (+0/-64)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%configs_and_props_info.dat (+0/-2012)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%leshouche_info.dat (+0/-129)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%matrix_1.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%matrix_2.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%matrix_5.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_gg_ttx%matrix_6.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uux_ttx%configs_and_props_info.dat (+0/-1668)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uux_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uux_ttx%leshouche_info.dat (+0/-147)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uxu_ttx%configs_and_props_info.dat (+0/-1668)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uxu_ttx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_pptt_fksrealew/%SubProcesses%P0_uxu_ttx%leshouche_info.dat (+0/-147)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%CT_interface.f (+7/-264)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f (+0/-748)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f (+4/-12)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%born_matrix.f (+2/-2)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%CT_interface.f (+7/-264)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%GOLEM_interface.f (+0/-748)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%TIR_interface.f (+4/-12)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%born_matrix.f (+2/-2)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f (+3/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%configs_and_props_info.dat (+0/-261)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_tdecay_fksreal/%SubProcesses%P0_t_budx%leshouche_info.dat (+0/-33)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%MadLoop5_resources%ColorDenomFactors.dat (+0/-50)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%MadLoop5_resources%ColorNumFactors.dat (+0/-50)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%MadLoop5_resources%HelConfigs.dat (+0/-16)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%CT_interface.f (+19/-270)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%GOLEM_interface.f (+0/-748)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%TIR_interface.f (+4/-12)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%born_matrix.f (+5/-5)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%coef_construction_1.f (+16/-14)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%f2py_wrapper.f (+37/-2)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%helas_calls_ampb_1.f (+7/-8)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%improve_ps.f (+24/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%loop_matrix.f (+21/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%mp_coef_construction_1.f (+16/-16)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%V0_dxu_veep%mp_helas_calls_ampb_1.f (+6/-7)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%configs_and_props_info.dat (+0/-347)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%leshouche_info.dat (+0/-43)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%matrix_1.f (+1/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%matrix_2.f (+1/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_dxu_veep%matrix_3.f (+1/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%MadLoop5_resources%ColorDenomFactors.dat (+0/-50)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%MadLoop5_resources%ColorNumFactors.dat (+0/-50)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%MadLoop5_resources%HelConfigs.dat (+0/-16)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%CT_interface.f (+19/-270)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%GOLEM_interface.f (+0/-748)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%TIR_interface.f (+4/-12)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%born_matrix.f (+5/-5)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%coef_construction_1.f (+16/-14)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%f2py_wrapper.f (+37/-2)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%helas_calls_ampb_1.f (+7/-8)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%improve_ps.f (+24/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%loop_matrix.f (+21/-3)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%mp_coef_construction_1.f (+17/-17)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%V0_udx_veep%mp_helas_calls_ampb_1.f (+6/-7)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%configs_and_props_info.dat (+0/-347)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%iproc.dat (+0/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%leshouche_info.dat (+0/-43)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%matrix_1.f (+1/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%matrix_2.f (+1/-1)
tests/input_files/IOTestsComparison/IOExportFKSTest/test_wprod_fksew/%SubProcesses%P0_udx_veep%matrix_3.f (+1/-1)
tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_standalone/matrix.f (+1/-1)
tests/input_files/IOTestsComparison/MECmdShell/check_html_long_process_strings/info.html (+80/-80)
tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%CT_interface.f (+23/-279)
tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%GOLEM_interface.f (+0/-823)
tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%TIR_interface.f (+4/-12)
tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%born_matrix.f (+2/-2)
tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%loop_matrix.f (+3/-3)
tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_NoSQSO.f (+1/-1)
tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_QCDsq_le_6.f (+2/-2)
tests/input_files/IOTestsComparison/SquaredOrder_IOTest/sqso_uux_uuxuuxx/matrix_ampOrderQED2_eq_2_WGTsq_le_14_QCDsq_gt_4.f (+2/-2)
tests/input_files/IOTestsComparison/TestCmdMatchBox/MatchBoxOutput/%TEST%SubProcesses%P0_wpwm_wpwm%matrix.f (+8/-8)
tests/input_files/IOTestsComparison/TestCmdMatchBox/MatchBoxOutput/%TEST%SubProcesses%P1_uux_uux%CT_interface.f (+15/-273)
tests/input_files/IOTestsComparison/TestCmdMatchBox/MatchBoxOutput/%TEST%SubProcesses%P1_uux_uux%GOLEM_interface.f (+0/-755)
tests/input_files/IOTestsComparison/TestCmdMatchBox/MatchBoxOutput/%TEST%SubProcesses%P1_uux_uux%TIR_interface.f (+4/-12)
tests/input_files/IOTestsComparison/TestCmdMatchBox/MatchBoxOutput/%TEST%SubProcesses%P1_uux_uux%loop_matrix.f (+3/-3)
tests/input_files/IOTestsComparison/long_ML_SMQCD_default/dux_mumvmxg/%..%..%Source%MODEL%model_functions.f (+0/-21)
tests/input_files/IOTestsComparison/long_ML_SMQCD_default/dux_mumvmxg/%..%..%Source%MODEL%model_functions.inc (+0/-2)
tests/input_files/IOTestsComparison/long_ML_SMQCD_default/dux_mumvmxg/born_matrix.f (+1/-1)
tests/input_files/IOTestsComparison/long_ML_SMQCD_default/gg_wmtbx/%..%..%Source%MODEL%model_functions.f (+0/-21)
tests/input_files/IOTestsComparison/long_ML_SMQCD_default/gg_wmtbx/%..%..%Source%MODEL%model_functions.inc (+0/-2)
tests/input_files/IOTestsComparison/long_ML_SMQCD_default/gg_wmtbx/born_matrix.f (+1/-1)
tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/dux_mumvmxg/%..%..%Source%MODEL%model_functions.f (+0/-21)
tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/dux_mumvmxg/%..%..%Source%MODEL%model_functions.inc (+0/-2)
tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/dux_mumvmxg/born_matrix.f (+1/-1)
tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/%..%..%Source%MODEL%model_functions.f (+0/-21)
tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/%..%..%Source%MODEL%model_functions.inc (+0/-2)
tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/born_matrix.f (+1/-1)
tests/input_files/IOTestsComparison/short_ML_SMQCD_LoopInduced/gg_hh/%..%..%Source%MODEL%model_functions.f (+0/-21)
tests/input_files/IOTestsComparison/short_ML_SMQCD_LoopInduced/gg_hh/%..%..%Source%MODEL%model_functions.inc (+0/-2)
tests/input_files/IOTestsComparison/short_ML_SMQCD_default/ddx_ttx/born_matrix.f (+1/-1)
tests/input_files/IOTestsComparison/short_ML_SMQCD_default/gg_ttx/%..%..%Source%MODEL%model_functions.f (+0/-21)
tests/input_files/IOTestsComparison/short_ML_SMQCD_default/gg_ttx/%..%..%Source%MODEL%model_functions.inc (+0/-2)
tests/input_files/IOTestsComparison/short_ML_SMQCD_default/gg_ttx/born_matrix.f (+1/-1)
tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/ddx_ttx/born_matrix.f (+1/-1)
tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/gg_ttx/%..%..%Source%MODEL%model_functions.f (+0/-21)
tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/gg_ttx/%..%..%Source%MODEL%model_functions.inc (+0/-2)
tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/gg_ttx/born_matrix.f (+1/-1)
tests/parallel_tests/madevent_comparator.py (+2/-2)
tests/parallel_tests/me_comparator.py (+3/-5)
tests/parallel_tests/test_ML5.py (+1/-1)
tests/parallel_tests/test_ML5MSSMQCD.py (+1/-1)
tests/parallel_tests/test_aloha.py (+168/-2)
tests/parallel_tests/test_cmd_amcatnlo.py (+5/-3)
tests/parallel_tests/test_madweight.py (+3/-3)
tests/test_manager.py (+1/-1)
tests/time_db (+80/-78)
tests/unit_tests/fks/test_extra_ew.py (+26/-3)
tests/unit_tests/fks/test_fks_base.py (+10/-10)
tests/unit_tests/iolibs/test_export_v4.py (+1/-1)
tests/unit_tests/iolibs/test_helas_call_writers.py (+30/-30)
tests/unit_tests/iolibs/test_link_to_ufo.py (+7/-2)
tests/unit_tests/iolibs/test_ufo_parsers.py (+56/-3)
tests/unit_tests/madspin/test_madspin.py (+7/-0)
tests/unit_tests/various/test_banner.py (+19/-0)
tests/unit_tests/various/test_decay.py (+2/-2)
tests/unit_tests/various/test_import_ufo.py (+138/-8)
tests/unit_tests/various/test_misc.py (+48/-7)
vendor/CutTools/src/makefile (+10/-4)
vendor/CutTools/src/qcdloop/aaxex.f (+1/-1)
vendor/IREGI/src/qcdloop/ff/aaxex.f (+1/-1)
vendor/IREGI/src/qcdloop/ff/makefile (+1/-1)
To merge this branch: bzr merge lp:~maddevelopers/mg5amcnlo/3.0.1
Reviewer Review Type Date Requested Status
Olivier Mattelaer 2018-10-04 Approve on 2018-12-19
marco zaro 2018-10-04 Needs Fixing on 2018-11-21
Hua-Sheng Shao 2018-10-04 Pending
Valentin Hirschi 2018-10-04 Pending
davide.pagani.85 2018-10-04 Pending
Review via email: mp+356133@code.launchpad.net

Commit message

Hi guys,

Before getting super excited and doing lots of test for the 4FS, we should release this code first...

Note: I proposed for merging with the FKS_EW_granny branch -- I don't know if this is exactly the branch relevant to the 3.0.0 release.

Cheers,
Rik

To post a comment you must log in.
marco zaro (marco-zaro) wrote :

have we enabled the NLO+PS for QCD only?

Rikkert Frederix (frederix) wrote :

It's not enabled. We should. Could you please do this?

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-10-04
942. By marco zaro on 2018-10-04

fixes in amcatnlo_run_interface.py to enable NLOPS for QCD-only
not fully working yet

marco zaro (marco-zaro) wrote :

I have pushed a fix, however, it is not complete yet (and the behaviour is quite strange)
for example, I cannot set fixed_order=OFF via the switches.
However, ./bin/generate_events aMC@NLO works as expected (and herwig is compiled and run too), however on the screen the printout one gets is
launch aMC@NLO
The following switches determine which programs are run:
/================== Description ===================|================= values =================|======== other options =========\
| 1. Type of perturbative computation | order = NLO | LO |
| 2. No MC@[N]LO matching / event generation | fixed_order = ON ⇐ ̶O̶F̶F̶ ̶ | ON |
| 3. Shower the generated events | shower = OFF ⇐ ̶H̶E̶R̶W̶I̶G̶6̶ ̶ | OFF|PYTHIA6Q|PYTHIA6PT |
| 4. Decay onshell particles | madspin = OFF | ON|onshell |
| 5. Add weights to events for new hypp. | reweight = OFF | ON|NLO|NLO_TREE|LO |
| 6. Run MadAnalysis5 on the events generated | madanalysis = Not Avail. | Please install module |
\==============================================================================================================================/
Either type the switch number (1 to 6) to change its setting,

which looks weird…
can anyone who knows better the business of these switches have a look?

thanks!
cheers,

Marco

> On 4 Oct 2018, at 15:07, Rikkert Frederix <email address hidden> wrote:
>
> It's not enabled. We should. Could you please do this?
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-10-04
943. By Rikkert Frederix on 2018-10-04

fixed the switches

Rikkert Frederix (frederix) wrote :

Hi,

A small change means that the switches are okay now for the QCD stuff. However, for QED, the formatting looks off. I'm not sure we should be too worried about this, since it works correctly. If there is a simple fix, we might apply it.

Cheers,
Rikkert

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-10-04
944. By marco zaro on 2018-10-04

small fix

Hi Rikkert,

The display is actually perfect for QED (see below):

The following switches determine which programs are run:
/================== Description ==================|=========== values ===========|================ other options ================\
| 1. Type of perturbative computation | order = NLO | LO |
| 2. No MC@[N]LO matching / event generation | fixed_order = ON | No NLO+PS available for EW correction |
\================================================================================================================================/
Either type the switch number (1 to 6) to change its setting,
Set any switch explicitly (e.g. type 'madspin=onshell' at the prompt)
Type 'help' for the list of all valid option
Type '0', 'auto', 'done' or just press enter when you are done.[60s to answer]

Maybe are you referring to the display in debug mode, which include the hidden line (madspin/...)
which are important to see (but only in debug mode). [Also for the hidden option of MadSpin]
In that case, formatting is indeed weird.
I have (slightly) improve it
1) passing the color to green (more consistent with debug coloring in general)
2) fixing the change of length due to the presence of the coloring
This is still not perfect since the length of each column does not take into account those hidden line.
and therefore they can go to overflow.

Cheers,

Olivier

> On 4 Oct 2018, at 17:16, Rikkert Frederix <email address hidden> wrote:
>
> Hi,
>
> A small change means that the switches are okay now for the QCD stuff. However, for QED, the formatting looks off. I'm not sure we should be too worried about this, since it works correctly. If there is a simple fix, we might apply it.
>
> Cheers,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-10-04
945. By olivier-mattelaer on 2018-10-04

fix some question formatting for hidden line [debgug mode only]

Download full text (4.6 KiB)

Hi,

I run (some of) the test this morning (actually 630 of thoses)
and this one sounds quite bad:

ERROR: test_generate_fks_ew (tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW)
check that the generate command works as expected.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/unit_tests/fks/test_extra_ew.py", line 112, in test_generate_fks_ew
    self.assertEqual(set(split), self.interface._fks_multi_proc['splitting_types'])
  File "/Users/omattelaer/Documents/workspace/3.0.1/madgraph/core/base_objects.py", line 69, in __getitem__
    self.is_valid_prop(name) #raise the correct error
  File "/Users/omattelaer/Documents/workspace/3.0.1/madgraph/core/base_objects.py", line 86, in is_valid_prop
    Valid property are %s""" % (name,self.__class__.__name__, self.keys())
PhysicsObjectError: splitting_types is not a valid property for this object: FKSMultiProcess

    Valid property are ['has_isr', 'use_numerical', 'amplitudes', 'born_processes', 'pdgs', 'ignore_six_quark_processes', 'collect_mirror_procs', 'process_definitions', 'OLP', 'real_amplitudes', 'init_lep_split', 'diagram_filter', 'ncores_for_proc_gen', 'has_fsr', 'loop_filter']

some other crashes like this one:

======================================================================
FAIL: testIO_test_wprod_fksew (tests.unit_tests.iolibs.test_export_fks.IOExportFKSTest)
target: SubProcesses/[P0.*\/.+\.(inc|f)]
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 257, in __wrapper
    testKeys=[(testGroup, newTestName)])
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 700, in runIOTests
    self.assertFileContains(open(file_path), goal)
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 389, in assertFileContains
    self.assertEqual(a,b)
AssertionError: ' CALL FFV1_2(W(1,4),W(1,2),GC_3,DCMPLX(ZERO),W(1,8))' != ' CALL FFV1_2(W(1,4),W(1,2),-GC_4,DCMPLX(ZERO),W(1,8))'

Which should be related to model optimization include in 2.6.3. I believe that this is fine (amy does not agree?)

Cheers,

Olivier

On 4 Oct 2018, at 23:24, Olivier Mattelaer <<email address hidden><mailto:<email address hidden>>> wrote:

Hi Rikkert,

The display is actually perfect for QED (see below):

The following switches determine which programs are run:
/================== Description ==================|=========== values ===========|================ other options ================\
| 1. Type of perturbative computation | order = NLO | LO |
| 2. No MC@[N]LO matching / event generation | fixed_order = ON | No NLO+PS available for EW correction |
\================================================================================================================================/
Either type the switch number (1 to 6) to change its setting,
Set any switch explicitly (e.g. type 'madspin=onshell' at the prompt)
T...

Read more...

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-10-05
946. By marco zaro on 2018-10-05

fixed one test

marco zaro (marco-zaro) wrote :
Download full text (5.4 KiB)

Hi Olivier,
the test you mentioned below was simply not updated after a change I did.
Now it is passing.

Cheers,

Marco

> On 5 Oct 2018, at 09:48, Olivier Mattelaer <email address hidden> wrote:
>
> Hi,
>
> I run (some of) the test this morning (actually 630 of thoses)
> and this one sounds quite bad:
>
> ERROR: test_generate_fks_ew (tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW)
> check that the generate command works as expected.
> ----------------------------------------------------------------------
> Traceback (most recent call last):
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/unit_tests/fks/test_extra_ew.py", line 112, in test_generate_fks_ew
> self.assertEqual(set(split), self.interface._fks_multi_proc['splitting_types'])
> File "/Users/omattelaer/Documents/workspace/3.0.1/madgraph/core/base_objects.py", line 69, in __getitem__
> self.is_valid_prop(name) #raise the correct error
> File "/Users/omattelaer/Documents/workspace/3.0.1/madgraph/core/base_objects.py", line 86, in is_valid_prop
> Valid property are %s""" % (name,self.__class__.__name__, self.keys())
> PhysicsObjectError: splitting_types is not a valid property for this object: FKSMultiProcess
>
> Valid property are ['has_isr', 'use_numerical', 'amplitudes', 'born_processes', 'pdgs', 'ignore_six_quark_processes', 'collect_mirror_procs', 'process_definitions', 'OLP', 'real_amplitudes', 'init_lep_split', 'diagram_filter', 'ncores_for_proc_gen', 'has_fsr', 'loop_filter']
>
> some other crashes like this one:
>
> ======================================================================
> FAIL: testIO_test_wprod_fksew (tests.unit_tests.iolibs.test_export_fks.IOExportFKSTest)
> target: SubProcesses/[P0.*\/.+\.(inc|f)]
> ----------------------------------------------------------------------
> Traceback (most recent call last):
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 257, in __wrapper
> testKeys=[(testGroup, newTestName)])
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 700, in runIOTests
> self.assertFileContains(open(file_path), goal)
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/IOTests.py", line 389, in assertFileContains
> self.assertEqual(a,b)
> AssertionError: ' CALL FFV1_2(W(1,4),W(1,2),GC_3,DCMPLX(ZERO),W(1,8))' != ' CALL FFV1_2(W(1,4),W(1,2),-GC_4,DCMPLX(ZERO),W(1,8))'
>
>
> Which should be related to model optimization include in 2.6.3. I believe that this is fine (amy does not agree?)
>
> Cheers,
>
> Olivier
>
> On 4 Oct 2018, at 23:24, Olivier Mattelaer <<email address hidden> <mailto:<email address hidden>><mailto:<email address hidden> <mailto:<email address hidden>>>> wrote:
>
> Hi Rikkert,
>
> The display is actually perfect for QED (see below):
>
>
> The following switches determine which programs are run:
> /================== Description ==================|=========== values ===========|================ other options ================\
> | 1. Type of perturbative computation | order = NLO | LO |
> | 2. No MC@[N]LO matching / e...

Read more...

Rikkert Frederix (frederix) wrote :

Hi,

There is another issue that should be fixed. In particular, when requiring e.g.,

p p > j j QED=0 QCD=4 [QED]

the code also includes QCD corrections (i.e., NLO_1), since at the squared order level, they are allowed. To remove them, one needs to set QCD=2.

In my opinion, the QED=0 and QCD=4 should be used to define which LO contributions should be included, and the ones between the square brackets to determine if one should include QCD and/or QED corrections on top of the Born contributions. Hence,

p p > j j QED=0 QCD=4 [QED]
and
p p > j j QED=0 QCD=2 [QED]

should lead to exactly the same results.

Who knows how to fix this?

Cheers,
Rikkert

On that topic,

Would not be a good idea to put somewhere the plot corresponding to NLO_X LO_Y?
This could be very usefull to know what is computed and what is included by the code.
(what about for 3.0.2?)

Cheers,

Olivier

> On 5 Oct 2018, at 10:48, Rikkert Frederix <email address hidden> wrote:
>
> Hi,
>
> There is another issue that should be fixed. In particular, when requiring e.g.,
>
> p p > j j QED=0 QCD=4 [QED]
>
> the code also includes QCD corrections (i.e., NLO_1), since at the squared order level, they are allowed. To remove them, one needs to set QCD=2.
>
> In my opinion, the QED=0 and QCD=4 should be used to define which LO contributions should be included, and the ones between the square brackets to determine if one should include QCD and/or QED corrections on top of the Born contributions. Hence,
>
> p p > j j QED=0 QCD=4 [QED]
> and
> p p > j j QED=0 QCD=2 [QED]
>
> should lead to exactly the same results.
>
> Who knows how to fix this?
>
> Cheers,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

marco zaro (marco-zaro) wrote :

ciao Rikkert,
I would prefer not to touch the order business now, as I fear to screw up other things with a quick patch.
I agree that it is not optimal, but we may want to upgrade it for the next release with more calm, possibly allowing the user to specify orders at the diagram or amplitude level in a more flexible way.
What is everybody’s idea here?

Cheers,

Marco

> On 5 Oct 2018, at 10:48, Rikkert Frederix <email address hidden> wrote:
>
> Hi,
>
> There is another issue that should be fixed. In particular, when requiring e.g.,
>
> p p > j j QED=0 QCD=4 [QED]
>
> the code also includes QCD corrections (i.e., NLO_1), since at the squared order level, they are allowed. To remove them, one needs to set QCD=2.
>
> In my opinion, the QED=0 and QCD=4 should be used to define which LO contributions should be included, and the ones between the square brackets to determine if one should include QCD and/or QED corrections on top of the Born contributions. Hence,
>
> p p > j j QED=0 QCD=4 [QED]
> and
> p p > j j QED=0 QCD=2 [QED]
>
> should lead to exactly the same results.
>
> Who knows how to fix this?
>
> Cheers,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

Rikkert Frederix (frederix) wrote :

Hi,

Marco: I think you are right. Maybe it's better to wait with this for the next version.

Olivier: I agree with you that for 3.0.2 we should write a better default analysis, which knows about the separate orders, and has a few more plots than just the total cross section.

Cheers,
Rikkert

I was actually not thinking at the analysis level but more at the time of the generation of the diagram
to have the famous bullet plot such that it is simple to see which order are include.

Olivier

> On 5 Oct 2018, at 11:22, Rikkert Frederix <email address hidden> wrote:
>
> Hi,
>
> Marco: I think you are right. Maybe it's better to wait with this for the next version.
>
> Olivier: I agree with you that for 3.0.2 we should write a better default analysis, which knows about the separate orders, and has a few more plots than just the total cross section.
>
> Cheers,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

davide.pagani.85 (davide-pagani) wrote :

Hi,

What Rik suggested is important, the example he is quoting is actually
an error that I did and and made loose quite some time to get what was
going on.

On the other hand, I think is more important to release as soon as
possible this new version and this error is actually due to a wrong
usage of the syntax.

So if it takes long, I would do it for next-to-next release.

-------------------------------------------------------

I think it would be very helpful to add the HwU analysis that Rik wrote
for the results in the paper.

Many times I start from that one or I use it to get fast results.

Do you agree?

Ciao

Davide

On 05/10/18 11:00, marco zaro wrote:
> ciao Rikkert,
> I would prefer not to touch the order business now, as I fear to screw up other things with a quick patch.
> I agree that it is not optimal, but we may want to upgrade it for the next release with more calm, possibly allowing the user to specify orders at the diagram or amplitude level in a more flexible way.
> What is everybody’s idea here?
>
> Cheers,
>
> Marco
>
>
>> On 5 Oct 2018, at 10:48, Rikkert Frederix <email address hidden> wrote:
>>
>> Hi,
>>
>> There is another issue that should be fixed. In particular, when requiring e.g.,
>>
>> p p > j j QED=0 QCD=4 [QED]
>>
>> the code also includes QCD corrections (i.e., NLO_1), since at the squared order level, they are allowed. To remove them, one needs to set QCD=2.
>>
>> In my opinion, the QED=0 and QCD=4 should be used to define which LO contributions should be included, and the ones between the square brackets to determine if one should include QCD and/or QED corrections on top of the Born contributions. Hence,
>>
>> p p > j j QED=0 QCD=4 [QED]
>> and
>> p p > j j QED=0 QCD=2 [QED]
>>
>> should lead to exactly the same results.
>>
>> Who knows how to fix this?
>>
>> Cheers,
>> Rikkert
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>

Rikkert Frederix (frederix) wrote :

Hi guys,

Any updates?

Best,
Rikkert

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-10-26
947. By Pagani <email address hidden> on 2018-10-26

General analysis with options of tagging spefic orders explained

davide.pagani.85 (davide-pagani) wrote :

Hi guys,

sorry for the naive question.
What does it mean merging 3.0.1 into FKS_EW_granny?

Isn't the latter an old version?

Ciao
Davide

On 23/10/18 21:51, Rikkert Frederix wrote:
> Hi guys,
>
> Any updates?
>
> Best,
> Rikkert
>

Rikkert Frederix (frederix) wrote :

FKS_EW_granny is the current public version ---there is no 3.0.0 branch. Hence, this request is to replace the current public version with the 3.0.1 branch.

Best,
Rikkert

davide.pagani.85 (davide-pagani) wrote :

Hi Marco,

I see

942. By marco zaro on 2018-10-04
fixes in amcatnlo_run_interface.py to enable NLOPS for QCD-only
not fully working yet

in which sense is not fully working?

Concerning the Fixed-Order part, the code is ok for me to be replaced.

Cheers
Davide

On 26/10/18 19:57, Rikkert Frederix wrote:
> FKS_EW_granny is the current public version ---there is no 3.0.0 branch. Hence, this request is to replace the current public version with the 3.0.1 branch.
>
> Best,
> Rikkert

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-10-28
948. By olivier-mattelaer on 2018-10-28

merge with 2.6.4

949. By olivier-mattelaer on 2018-10-28

merge with 2.6.4

marco zaro (marco-zaro) wrote :

Ciao Davide,
this was fixed in rev 943 by rikkert
cheers,

Marco

> On 28 Oct 2018, at 09:15, davide.pagani.85 <email address hidden> wrote:
>
> Hi Marco,
>
> I see
>
> 942. By marco zaro on 2018-10-04
> fixes in amcatnlo_run_interface.py to enable NLOPS for QCD-only
> not fully working yet
>
> in which sense is not fully working?
>
> Concerning the Fixed-Order part, the code is ok for me to be replaced.
>
>
> Cheers
> Davide
>
> On 26/10/18 19:57, Rikkert Frederix wrote:
>> FKS_EW_granny is the current public version ---there is no 3.0.0 branch. Hence, this request is to replace the current public version with the 3.0.1 branch.
>>
>> Best,
>> Rikkert
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

davide.pagani.85 (davide-pagani) wrote :

Then I am ok for the merging.

Let me know how I can help you.

Ciao
Davide

On 29/10/18 09:21, marco zaro wrote:
> Ciao Davide,
> this was fixed in rev 943 by rikkert
> cheers,
>
> Marco
>
>> On 28 Oct 2018, at 09:15, davide.pagani.85 <email address hidden> wrote:
>>
>> Hi Marco,
>>
>> I see
>>
>> 942. By marco zaro on 2018-10-04
>> fixes in amcatnlo_run_interface.py to enable NLOPS for QCD-only
>> not fully working yet
>>
>> in which sense is not fully working?
>>
>> Concerning the Fixed-Order part, the code is ok for me to be replaced.
>>
>>
>> Cheers
>> Davide
>>
>> On 26/10/18 19:57, Rikkert Frederix wrote:
>>> FKS_EW_granny is the current public version ---there is no 3.0.0 branch. Hence, this request is to replace the current public version with the 3.0.1 branch.
>>>
>>> Best,
>>> Rikkert
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-11-09
950. By olivier-mattelaer on 2018-11-09

UPdate IOtest

951. By olivier-mattelaer on 2018-11-09

last merge with 2.6.4

Hi Marco,

I realised that you removed that file:
tests/acceptance_tests/test_export_fks.py

Is this expected?

Olivier

Also it seems to have a bunch of (scary) print statement in
test_fks_base.py at line 246
Related to some bypass of tests...

Would be better to clean that right?

Olivier

marco zaro (marco-zaro) wrote :

Hi Olviier,
all IOtests for the NLO output are in
tests/unit_tests/iolibs/test_export_fks.py
feel free to move it to the acceptance tests if it takes too much

cheers,

Marco

> On 9 Nov 2018, at 11:31, Olivier Mattelaer <email address hidden> wrote:
>
> Hi Marco,
>
> I realised that you removed that file:
> tests/acceptance_tests/test_export_fks.py
>
> Is this expected?
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

OK good and concerning the test_fks_base.py?
(print statement that some test are bypassed?)

Cheers,

Olivier

PS: I have found workaround for the lhapdf issue and Mojave issue that I pushed in 2.6.4 (and include those here obviously)

> On 9 Nov 2018, at 12:15, marco zaro <email address hidden> wrote:
>
> Hi Olviier,
> all IOtests for the NLO output are in
> tests/unit_tests/iolibs/test_export_fks.py
> feel free to move it to the acceptance tests if it takes too much
>
> cheers,
>
> Marco
>
>> On 9 Nov 2018, at 11:31, Olivier Mattelaer <email address hidden> wrote:
>>
>> Hi Marco,
>>
>> I realised that you removed that file:
>> tests/acceptance_tests/test_export_fks.py
>>
>> Is this expected?
>>
>> Olivier
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-11-09
952. By marco zaro on 2018-11-09

commented sprint statament in test_fks_base

marco zaro (marco-zaro) wrote :

these are commented now (r952)
which LHAPDF issue?
cheers,

Marco

> On 9 Nov 2018, at 12:30, Olivier Mattelaer <email address hidden> wrote:
>
> OK good and concerning the test_fks_base.py?
> (print statement that some test are bypassed?)
>
>
> Cheers,
>
> Olivier
>
>
> PS: I have found workaround for the lhapdf issue and Mojave issue that I pushed in 2.6.4 (and include those here obviously)
>
>> On 9 Nov 2018, at 12:15, marco zaro <email address hidden> wrote:
>>
>> Hi Olviier,
>> all IOtests for the NLO output are in
>> tests/unit_tests/iolibs/test_export_fks.py
>> feel free to move it to the acceptance tests if it takes too much
>>
>> cheers,
>>
>> Marco
>>
>>> On 9 Nov 2018, at 11:31, Olivier Mattelaer <email address hidden> wrote:
>>>
>>> Hi Marco,
>>>
>>> I realised that you removed that file:
>>> tests/acceptance_tests/test_export_fks.py
>>>
>>> Is this expected?
>>>
>>> Olivier
>>> --
>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>>
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

They were two lhapdf issues:

- it was not compiling on Mojave (the same problem will likely re-appear in a year with NO easy solution this time). The problem is that the yaml parser shipped with lhapdf seems to have issue on mac.

- the download of (any) PDF was crashing. This was due to the change of hepforge.
We were using the ofiicial command (lhapdf install) but that one does not work anymore and after discussing with them they do not seem interested to fix it.
I'm still trying to use that command, but if it fails then I guess the path of the set and do a "simple" wget on that one.

Cheers,

Olivier

PS: The status on the tests are now:
- one unnitest failing
- around 10 acceptance test failing (so far, this is still running)

> On 9 Nov 2018, at 12:42, marco zaro <email address hidden> wrote:
>
> these are commented now (r952)
> which LHAPDF issue?
> cheers,
>
> Marco
>
>> On 9 Nov 2018, at 12:30, Olivier Mattelaer <email address hidden> wrote:
>>
>> OK good and concerning the test_fks_base.py?
>> (print statement that some test are bypassed?)
>>
>>
>> Cheers,
>>
>> Olivier
>>
>>
>> PS: I have found workaround for the lhapdf issue and Mojave issue that I pushed in 2.6.4 (and include those here obviously)
>>
>>> On 9 Nov 2018, at 12:15, marco zaro <email address hidden> wrote:
>>>
>>> Hi Olviier,
>>> all IOtests for the NLO output are in
>>> tests/unit_tests/iolibs/test_export_fks.py
>>> feel free to move it to the acceptance tests if it takes too much
>>>
>>> cheers,
>>>
>>> Marco
>>>
>>>> On 9 Nov 2018, at 11:31, Olivier Mattelaer <email address hidden> wrote:
>>>>
>>>> Hi Marco,
>>>>
>>>> I realised that you removed that file:
>>>> tests/acceptance_tests/test_export_fks.py
>>>>
>>>> Is this expected?
>>>>
>>>> Olivier
>>>> --
>>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>>>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>>>
>>>
>>> --
>>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>>
>>
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

Hi Marco,

The following test is also not passing (going trough ok with 2.6.4):
test_madspin_LOonly

Even though this is a madspin designed tests, it actually crash within the LOonly module.
If you just do
geneate p p > w+ [LOonly]
It will lead to a crash.

I have fixed all the other tests (but a couple of IOtests and i still have to run the short paralel test).

Since this is a beta version, we can also decide to not fix this in this version.

marco zaro (marco-zaro) wrote :

Hi Olivier,
it is strange, it seems that the LOonly mode is broken (it always returns zero)
I will have a look on monday, and let you know.

cheers,

Marco

> On 9 Nov 2018, at 21:10, Olivier Mattelaer <email address hidden> wrote:
>
> Hi Marco,
>
> The following test is also not passing (going trough ok with 2.6.4):
> test_madspin_LOonly
>
> Even though this is a madspin designed tests, it actually crash within the LOonly module.
> If you just do
> geneate p p > w+ [LOonly]
> It will lead to a crash.
>
> I have fixed all the other tests (but a couple of IOtests and i still have to run the short paralel test).
>
> Since this is a beta version, we can also decide to not fix this in this version.
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

marco zaro (marco-zaro) wrote :

Hi Again,
actually, doing p p > e+ ve [LOonly] works ok, at LO, NLO, aMC@LO, aMC@NLO
So it may be something relate to the phase-space of 2->1?
the reason why p p > w+ [LOonly] returns zero is because the function passcuts always returns False, exiting here
c Make sure have reasonable 4-momenta
      if (p(0,1) .le. 0d0) then
         passcuts=.false.
         return
      endif

with P(0,1) = -99

So clearly something goes wrong in the phase space generation. Rikkert, anything that could be related to the granny stuff?

also note that
p p > w+ [QED] has one test failing (soft test for photon emission off the W), where a lot of NaN’s appear
p p > e+ ve [QCD] seems ok though

Cheers,

Marco

> On 9 Nov 2018, at 21:42, marco zaro <email address hidden> wrote:
>
> Hi Olivier,
> it is strange, it seems that the LOonly mode is broken (it always returns zero)
> I will have a look on monday, and let you know.
>
> cheers,
>
> Marco
>
>> On 9 Nov 2018, at 21:10, Olivier Mattelaer <email address hidden> wrote:
>>
>> Hi Marco,
>>
>> The following test is also not passing (going trough ok with 2.6.4):
>> test_madspin_LOonly
>>
>> Even though this is a madspin designed tests, it actually crash within the LOonly module.
>> If you just do
>> geneate p p > w+ [LOonly]
>> It will lead to a crash.
>>
>> I have fixed all the other tests (but a couple of IOtests and i still have to run the short paralel test).
>>
>> Since this is a beta version, we can also decide to not fix this in this version.
>> --
>> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
>> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

Thanks,

Update on my side:
1) the IOTest suite are bugged (not critical for a release) I have ask Valentin to look at it
2) the (fast) parrelel test leads to their own discovery of bug:

2.1) e+ e- > t t~ [real=QCD]
crashes with
FKSProcessError: This process does not have any correction up to NLO in QED,QCD
[ test: test_short_check_generate_events_nlo_py6pt_fsr]
same for e+ e- > p p p [real=QCD] in test: test_short_check_eejjj_lo_lhapdf]

2.2) A division by zero occurs in the jet veto of rik [test_short_jet_veto_xsec]

p p > e+ ve [QCD]
generate_events NLO
                 set ickkw -1
                 set ptj 10

2.3) The value of some squared order amplitude does not match anymore [test_short_sqso ]
Will look at this one (maybe with Valentin).

Cheers,

Olivier

On 9 Nov 2018, at 22:00, marco zaro <<email address hidden><mailto:<email address hidden>>> wrote:

Hi Again,
actually, doing p p > e+ ve [LOonly] works ok, at LO, NLO, aMC@LO, aMC@NLO
So it may be something relate to the phase-space of 2->1?
the reason why p p > w+ [LOonly] returns zero is because the function passcuts always returns False, exiting here
c Make sure have reasonable 4-momenta
     if (p(0,1) .le. 0d0) then
        passcuts=.false.
        return
     endif

with P(0,1) = -99

So clearly something goes wrong in the phase space generation. Rikkert, anything that could be related to the granny stuff?

also note that
p p > w+ [QED] has one test failing (soft test for photon emission off the W), where a lot of NaN’s appear
p p > e+ ve [QCD] seems ok though

Cheers,

Marco

On 9 Nov 2018, at 21:42, marco zaro <<email address hidden><mailto:<email address hidden>>> wrote:

Hi Olivier,
it is strange, it seems that the LOonly mode is broken (it always returns zero)
I will have a look on monday, and let you know.

cheers,

Marco

On 9 Nov 2018, at 21:10, Olivier Mattelaer <<email address hidden><mailto:<email address hidden>>> wrote:

Hi Marco,

The following test is also not passing (going trough ok with 2.6.4):
test_madspin_LOonly

Even though this is a madspin designed tests, it actually crash within the LOonly module.
If you just do
geneate p p > w+ [LOonly]
It will lead to a crash.

I have fixed all the other tests (but a couple of IOtests and i still have to run the short paralel test).

Since this is a beta version, we can also decide to not fix this in this version.
--
https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

--
https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

--
https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

marco zaro (marco-zaro) wrote :
Download full text (4.1 KiB)

Ciao Olivier,
thanks for looking at this.
Let me answer inline

On Fri, Nov 9, 2018 at 10:30 PM Olivier Mattelaer <
<email address hidden>> wrote:

> Thanks,
>
> Update on my side:
> 1) the IOTest suite are bugged (not critical for a release) I have ask
> Valentin to look at it
> 2) the (fast) parrelel test leads to their own discovery of bug:
>
> 2.1) e+ e- > t t~ [real=QCD]
> crashes with
> FKSProcessError: This process does not have any correction up to NLO in
> QED,QCD
> [ test: test_short_check_generate_events_nlo_py6pt_fsr]
> same for e+ e- > p p p [real=QCD] in test:
> test_short_check_eejjj_lo_lhapdf]
>
> for ee collisions, one has to tell the code
set include_lepton_initiated_processes True (there is a message printed on
screen).
The py6isr test should be skipped, as the has_isr flag is not updated
anymore in this branch (there were some print statement about this you
asked me to remove)

  If you want to include it, set the

  'include_lepton_initiated_processes' option to True

> 2.2) A division by zero occurs in the jet veto of rik
> [test_short_jet_veto_xsec]
>
> p p > e+ ve [QCD]
> generate_events NLO
> set ickkw -1
> set ptj 10
>
>
> the jet_veto is likely to be broken in this branch. There are some message
errors in the logs of the fortran executables about this.

> 2.3) The value of some squared order amplitude does not match anymore
> [test_short_sqso ]
> Will look at this one (maybe with Valentin).
>

I will check this right away

Thanks again!

cheers,

Marco

>
> Cheers,
>
> Olivier
>
> On 9 Nov 2018, at 22:00, marco zaro <<email address hidden><mailto:
> <email address hidden>>> wrote:
>
> Hi Again,
> actually, doing p p > e+ ve [LOonly] works ok, at LO, NLO, aMC@LO, aMC@NLO
> So it may be something relate to the phase-space of 2->1?
> the reason why p p > w+ [LOonly] returns zero is because the function
> passcuts always returns False, exiting here
> c Make sure have reasonable 4-momenta
> if (p(0,1) .le. 0d0) then
> passcuts=.false.
> return
> endif
>
> with P(0,1) = -99
>
> So clearly something goes wrong in the phase space generation. Rikkert,
> anything that could be related to the granny stuff?
>
> also note that
> p p > w+ [QED] has one test failing (soft test for photon emission off the
> W), where a lot of NaN’s appear
> p p > e+ ve [QCD] seems ok though
>
> Cheers,
>
> Marco
>
>
>
> On 9 Nov 2018, at 21:42, marco zaro <<email address hidden><mailto:
> <email address hidden>>> wrote:
>
> Hi Olivier,
> it is strange, it seems that the LOonly mode is broken (it always returns
> zero)
> I will have a look on monday, and let you know.
>
> cheers,
>
> Marco
>
> On 9 Nov 2018, at 21:10, Olivier Mattelaer <<email address hidden>
> <mailto:<email address hidden>>> wrote:
>
> Hi Marco,
>
> The following test is also not passing (going trough ok with 2.6.4):
> test_madspin_LOonly
>
> Even though this is a madspin designed tests, it actually crash within the
> LOonly module.
> If you just do
> geneate p p > w+ [LOonly]
> It will lead to a crash.
>
> I have fixed all the other tests (but a couple of IOtests and i still have
> to run the short ...

Read more...

Rikkert Frederix (frederix) wrote :

Hi Olivier, Marco,

I looked into the following two issues:

1. The error in the phase-space for the process 'p p > w+ [LOonly=QCD]'.
This is related to the fact that in the 3.0.1 branch the final state W boson is the fks_j particle, while in the 2.6.4 branch it's one of the initial state quarks. However, the phase-space for 2->1 processes with j-fks the massive final state particle does not work. There are divisions by zero everywhere... I'm not sure if there is an easy fix here, apart from telling the code to use an initial state particle as fake j_fks. This will work for the LOonly 2->1 processes. On the other hand, we need to fix this phase-space at some point, since the same is used (and needed!) for the process 'p p > w+ [QED]'.

2. The error for the VETO X-SEC (ickkw=-1).
The whole split-amplitude business has not been implemented for the veto-xsec stuff. I suggest that we do not try to rush this through now, but leave the 3.0.1 code as beta and fix this at a later stage.

Cheers,
Rik

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-11-16
953. By marco zaro on 2018-11-12

fix for 2->1 processes

954. By marco zaro on 2018-11-12

unit test for previous push and fix for p p > w+ [LOonly=QCD]

955. By marco zaro on 2018-11-16

reverted change from revision 613 (related to split orders and WEIGHTED)
added a debug message

marco zaro (marco-zaro) wrote :

Hi all,
I have just noticed this (fks_singular, line 3818)
      call Qterms_reduced_timelike(j_type,i_type,ch_i,ch_j,t,z,Q)
it should be instead
      call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,t,z,Q)
(switch ch_i <--> ch_j), shouldn't it?

It may be that in the end there is no effect, since, for the relevant cases when either ch_i or ch_j ==0 and the other is not one has always Q=0

c g/a ->qqbar splitting
         Qterms(1) = 4d0 * TR * z*(1d0-z)**2
         Qterms(2) = 4d0 * dble(abs(col1)) * ch1**2 * z*(1d0-z)**2

      elseif ((abs(col1).eq.3 .and. col2.eq.8) .or.
     & (dabs(ch1).gt.0d0 .and. dabs(ch2).eq.0d0)) then
c q->q g/a splitting
         Qterms(1) = 0d0
         Qterms(2) = 0d0

Let me know

Cheers,

Marco

review: Needs Fixing
marco zaro (marco-zaro) wrote :

sorry, the piece of code i wanted to post is

      elseif ((abs(col1).eq.3 .and. col2.eq.8) .or.
     & (dabs(ch1).gt.0d0 .and. dabs(ch2).eq.0d0)) then
c q->q g/a splitting
         Qterms(1) = 0d0
         Qterms(2) = 0d0

      elseif ((col1.eq.8 .and. abs(col2).eq.3) .or.
     & (dabs(ch1).eq.0d0 .and. dabs(ch2).gt.0d0)) then
c q->g/a q splitting
         Qterms(1) = 0d0
         Qterms(2) = 0d0
      else

> On 21 Nov 2018, at 16:16, marco zaro <email address hidden> wrote:
>
> Review: Needs Fixing
>
> Hi all,
> I have just noticed this (fks_singular, line 3818)
> call Qterms_reduced_timelike(j_type,i_type,ch_i,ch_j,t,z,Q)
> it should be instead
> call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,t,z,Q)
> (switch ch_i <--> ch_j), shouldn't it?
>
> It may be that in the end there is no effect, since, for the relevant cases when either ch_i or ch_j ==0 and the other is not one has always Q=0
>
>
> c g/a ->qqbar splitting
> Qterms(1) = 4d0 * TR * z*(1d0-z)**2
> Qterms(2) = 4d0 * dble(abs(col1)) * ch1**2 * z*(1d0-z)**2
>
> elseif ((abs(col1).eq.3 .and. col2.eq.8) .or.
> & (dabs(ch1).gt.0d0 .and. dabs(ch2).eq.0d0)) then
> c q->q g/a splitting
> Qterms(1) = 0d0
> Qterms(2) = 0d0
>
> Let me know
>
> Cheers,
>
> Marco
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

Rikkert Frederix (frederix) wrote :

Hi Marco,

I agree that the current way is wrong, but has no effect on the running of the code whatsoever. Even though, we should replace ch_i and ch_j in the call to the Qterms.
Note that, if there was 'true' bug here, the soft and/or collinear checks and the beginning of any run would not have passed.

Cheers,
Rik

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-11-21
956. By marco zaro on 2018-11-21

fix in a call to
Qterms_reduced_timelike
luckiliy harmless

marco zaro (marco-zaro) wrote :

Hi Rik,
thanks for the cross-check
I have pushed the fix into 3.0.1, rev 956
cheers,

Marco

> On 21 Nov 2018, at 16:24, Rikkert Frederix <email address hidden> wrote:
>
> Hi Marco,
>
> I agree that the current way is wrong, but has no effect on the running of the code whatsoever. Even though, we should replace ch_i and ch_j in the call to the Qterms.
> Note that, if there was 'true' bug here, the soft and/or collinear checks and the beginning of any run would not have passed.
>
> Cheers,
> Rik
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

Hi,

I would like to put this quite high priority for this week (i.e. for me Wednesday-friday).

I have rerun the tests and found that some tests are crashing when run via the full script but are working if run alone. So I have run (the very slow script to determine which of the tests had such side effects and the culprint is:
tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW.test_generate_fks_2to1_no_finalstate_confs
which has a side effect that lead to the following crash:

FAIL: test_generate_fks_ew_dijet_qed2qcd2_qcd (tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW)
check that the generate command works as expected.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/omattelaer/Documents/workspace/3.0.1/tests/unit_tests/fks/test_extra_ew.py", line 869, in test_generate_fks_ew_dijet_qed2qcd2_qcd
    self.assertEqual(real.fks_infos[0]['extra_cnt_index'], -1)
AssertionError: 0 != -1

Might be interesting that you take a look. But I think that we should release this version this week independently.

Any objection? comment?

Olivier

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-12-17
957. By marco zaro on 2018-12-17

tests/unit_tests/fks/test_extra_ew.py has been modified to avoid
side effects

958. By marco zaro on 2018-12-17

fix in amcatnlo_interface to prevent unphysical coupling orders

marco zaro (marco-zaro) wrote :

Hello Olivier,
I agree we should release 3.0.1 before the Christmas break.
concerning the side effect, i reproduced it, and it goes away if I do the following change
 class TestAMCatNLOEW(unittest.TestCase):
     """ a suite of extra tests for the ew stuff """

- interface = mgcmd.MasterCmd()
+ def setUp(self):
+ self.interface = mgcmd.MasterCmd()

I have pushed the fix into rev 957.

The other bogus behavious was related to the order settings in t-channel single top, which we discussed in a separate email.
I have added an error message that should warn the user to specify the coupling orders in that case (rev 958)

Cheers,

Marco

> On 17 Dec 2018, at 09:23, Olivier Mattelaer <email address hidden> wrote:
>
> Hi,
>
> I would like to put this quite high priority for this week (i.e. for me Wednesday-friday).
>
> I have rerun the tests and found that some tests are crashing when run via the full script but are working if run alone. So I have run (the very slow script to determine which of the tests had such side effects and the culprint is:
> tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW.test_generate_fks_2to1_no_finalstate_confs
> which has a side effect that lead to the following crash:
>
> FAIL: test_generate_fks_ew_dijet_qed2qcd2_qcd (tests.unit_tests.fks.test_extra_ew.TestAMCatNLOEW)
> check that the generate command works as expected.
> ----------------------------------------------------------------------
> Traceback (most recent call last):
> File "/Users/omattelaer/Documents/workspace/3.0.1/tests/unit_tests/fks/test_extra_ew.py", line 869, in test_generate_fks_ew_dijet_qed2qcd2_qcd
> self.assertEqual(real.fks_infos[0]['extra_cnt_index'], -1)
> AssertionError: 0 != -1
>
> Might be interesting that you take a look. But I think that we should release this version this week independently.
>
> Any objection? comment?
>
> Olivier
>
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

Rikkert Frederix (frederix) wrote :

Hi Olivier,

Thanks! We should really release it before the break. It will remain a beta version and exist along side the main branch, right? I think that would still be better and removing the beta tag.

best,
Rikkert

Hi,

> Thanks! We should really release it before the break. It will remain a beta version and exist along side the main branch, right?

Correct.

> I think that would still be better and removing the beta tag.

It would not make sense to remove the "beta" tag and keep two official version.
And I do not think that we are ready to keep only this version.

Cheers,

Olivier

> On 18 Dec 2018, at 10:33, Rikkert Frederix <email address hidden> wrote:
>
> Hi Olivier,
>
> Thanks! We should really release it before the break. It will remain a beta version and exist along side the main branch, right? I think that would still be better and removing the beta tag.
>
> best,
> Rikkert
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

Rikkert Frederix (frederix) wrote :

Sorry, I forgot a "NOT" in my suggestion :-D
Hence I agree with not removing the beta tag.

Cheers,
Rik

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-12-19
959. By olivier-mattelaer on 2018-12-19

fixing some additional issue with some test --in particular a real bug for the reweighting--

I have re-run all the tests today (and fix a couple of small stuffs). Except a couple that we have decide to ignore above they all go trough now. So I will make the release tomorrow.

Cheers,

Olivier

review: Approve
Rikkert Frederix (frederix) wrote :

Thanks a lot, Olivier.

Cheers,
Rik

marco zaro (marco-zaro) wrote :

Thanks a lot Olivier!
great that we finally made it!

cheers,

Marco

> On 19 Dec 2018, at 23:42, Rikkert Frederix <email address hidden> wrote:
>
> Thanks a lot, Olivier.
>
> Cheers,
> Rik
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> You are reviewing the proposed merge of lp:~maddevelopers/mg5amcnlo/3.0.1 into lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.

lp:~maddevelopers/mg5amcnlo/3.0.1 updated on 2018-12-20
960. By olivier-mattelaer on 2018-12-20

UpdateNotes/Version update

Done!

davide.pagani.85 (davide-pagani) wrote :

Thanks a lot!

Cheers
Davide

On 20/12/18 09:47, Olivier Mattelaer wrote:
> Done!

By the way,

Do you know if we can close the following bug report:
https://bugs.launchpad.net/mg5amcnlo/+bug/1779831

If not we should at least comment on it (this is one was forbidding ATLAS to use the code)

Cheers,

Olivier

Stefano Frixione (stefano-frixione) wrote :

I have no idea.
Perhaps one should double check with the ATLAS people;
I've never understood how genser works, and their relationship
with experiments.
Cheers, Stefano.

On Fri, 21 Dec 2018, Olivier Mattelaer wrote:

> By the way,
>
> Do you know if we can close the following bug report:
> https://bugs.launchpad.net/mg5amcnlo/+bug/1779831
>
> If not we should at least comment on it (this is one was forbidding
> ATLAS to use the code)
>
>
> Cheers,
>
> Olivier
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/+merge/356133
> Your team MadDevelopers is subscribed to branch lp:~maddevelopers/mg5amcnlo/FKS_EW_granny.
>

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py 2018-04-25 06:47:06 +0000
+++ MadSpin/decay.py 2018-12-20 08:18:55 +0000
@@ -3963,14 +3963,23 @@
3963 elif mode=='full':3963 elif mode=='full':
3964 stdin_text="5 0 0 0 \n" # before closing, write down the seed 3964 stdin_text="5 0 0 0 \n" # before closing, write down the seed
3965 external = self.calculator[('full',path)]3965 external = self.calculator[('full',path)]
3966 external.stdin.write(stdin_text)3966 try:
3967 external.stdin.write(stdin_text)
3968 except Exception:
3969 continue
3967 ranmar_state=external.stdout.readline()3970 ranmar_state=external.stdout.readline()
3968 ranmar_file=pjoin(path,'ranmar_state.dat')3971 ranmar_file=pjoin(path,'ranmar_state.dat')
3969 ranmar=open(ranmar_file, 'w')3972 ranmar=open(ranmar_file, 'w')
3970 ranmar.write(ranmar_state)3973 ranmar.write(ranmar_state)
3971 ranmar.close()3974 ranmar.close()
3972 external.stdin.close()3975 try:
3973 external.stdout.close()3976 external.stdin.close()
3977 except Exception:
3978 continue
3979 try:
3980 external.stdout.close()
3981 except Exception:
3982 continue
3974 external.terminate()3983 external.terminate()
3975 del external3984 del external
3976 else:3985 else:
39773986
=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py 2018-03-11 12:10:19 +0000
+++ MadSpin/interface_madspin.py 2018-12-20 08:18:55 +0000
@@ -35,6 +35,7 @@
35import madgraph.interface.madgraph_interface as mg_interface35import madgraph.interface.madgraph_interface as mg_interface
36import madgraph.interface.master_interface as master_interface36import madgraph.interface.master_interface as master_interface
37import madgraph.interface.madevent_interface as madevent_interface37import madgraph.interface.madevent_interface as madevent_interface
38import madgraph.interface.common_run_interface as common_run_interface
38import madgraph.interface.reweight_interface as rwgt_interface39import madgraph.interface.reweight_interface as rwgt_interface
39import madgraph.various.misc as misc40import madgraph.various.misc as misc
40import madgraph.iolibs.files as files41import madgraph.iolibs.files as files
@@ -53,6 +54,68 @@
53cmd_logger = logging.getLogger('cmdprint2') # -> print54cmd_logger = logging.getLogger('cmdprint2') # -> print
5455
5556
57class MadSpinOptions(banner.ConfigFile):
58
59 def default_setup(self):
60
61 self.add_param("max_weight", -1)
62 self.add_param('curr_dir', os.path.realpath(os.getcwd()))
63 self.add_param('Nevents_for_max_weigth', 0)
64 self.add_param("max_weight_ps_point", 400)
65 self.add_param('BW_cut', -1)
66 self.add_param('nb_sigma', 0.)
67 self.add_param('ms_dir', '')
68 self.add_param('max_running_process', 100)
69 self.add_param('onlyhelicity', False)
70 self.add_param('spinmode', "madspin", allowed=['madspin','none','onshell'])
71 self.add_param('use_old_dir', False, comment='should be use only for faster debugging')
72 self.add_param('run_card', '' , comment='define cut for spinmode==none. Path to run_card to use')
73 self.add_param('fixed_order', False, comment='to activate fixed order handling of counter-event')
74 self.add_param('seed', 0, comment='control the seed of madspin')
75 self.add_param('cross_section', {'__type__':0.}, comment="forcing normalization of cross-section after MS (for none/onshell)" )
76 self.add_param('new_wgt', 'cross-section' ,allowed=['cross-section', 'BR'], comment="if not consistent number of particles, choose what to do for the weight. (BR: means local according to number of part, cross use the force cross-section")
77 self.add_param('input_format', 'auto', allowed=['auto','lhe', 'hepmc', 'lhe_no_banner'])
78
79
80 ############################################################################
81 ## Special post-processing of the options ##
82 ############################################################################
83 def post_set_ms_dir(self, value, change_userdefine, raiseerror):
84 """ special handling for set ms_dir """
85
86 self.__setitem__('curr_dir', value, change_userdefine=change_userdefine)
87
88 ############################################################################
89 def post_set_seed(self, value, change_userdefine, raiseerror):
90 """ special handling for set seed """
91
92 random.seed(value)
93
94 ############################################################################
95 def post_set_run_card(self, value, change_userdefine, raiseerror):
96 """ special handling for set run_card """
97
98 if value == 'default':
99 self.run_card = None
100 elif os.path.isfile(value):
101 self.run_card = banner.RunCard(value)
102
103 args = value.split()
104 if len(args) >2:
105 if not self.options['run_card']:
106 self.run_card = banner.RunCardLO()
107 self.run_card.remove_all_cut()
108 self.run_card[args[0]] = ' '.join(args[1:])
109
110
111 ############################################################################
112 def post_fixed_order(self, value, change_userdefine, raiseerror):
113 """ special handling for set fixed_order """
114
115 if value:
116 logger.warning('Fix order madspin fails to have the correct scale information. This can bias the results!')
117 logger.warning('Not all functionalities of MadSpin handle this mode correctly (only onshell mode so far).')
118
56119
57class MadSpinInterface(extended_cmd.Cmd):120class MadSpinInterface(extended_cmd.Cmd):
58 """Basic interface for madspin"""121 """Basic interface for madspin"""
@@ -77,22 +140,7 @@
77 self.mode = "madspin" # can be flat/bridge change the way the decay is done.140 self.mode = "madspin" # can be flat/bridge change the way the decay is done.
78 # note amc@nlo does not support bridge.141 # note amc@nlo does not support bridge.
79 142
80 self.options = {'max_weight': -1, 143 self.options = MadSpinOptions()
81 'curr_dir': os.path.realpath(os.getcwd()),
82 'Nevents_for_max_weigth': 0,
83 'max_weight_ps_point': 400,
84 'BW_cut':-1,
85 'nb_sigma':0,
86 'ms_dir':None,
87 'max_running_process':100,
88 'onlyhelicity': False,
89 'spinmode': "madspin",
90 'use_old_dir': False, #should be use only for faster debugging
91 'run_card': None, # define cut for spinmode==none.
92 'fixed_order': False # to activate fixed order handling of counter-event
93 }
94
95
96 144
97 self.events_file = None145 self.events_file = None
98 self.decay_processes = {}146 self.decay_processes = {}
@@ -107,7 +155,21 @@
107 logger.info("Extracting the banner ...")155 logger.info("Extracting the banner ...")
108 self.do_import(event_path)156 self.do_import(event_path)
109 157
110 158
159 def setup_for_pure_decay(self):
160 """this is for spinmode=None -> simple decay
161 We go here if they are no banner.
162 -> this requires that a command import model appears in the card!
163 """
164
165 logger.info("Setup the code for pure decay mode")
166 self.proc_option = []
167 self.final_state_full = ''
168 self.final_state_compact = ''
169 self.prod_branches = ''
170 self.final_state = set()
171
172
111 def do_import(self, inputfile):173 def do_import(self, inputfile):
112 """import the event file"""174 """import the event file"""
113 175
@@ -126,6 +188,7 @@
126 if not os.path.exists(inputfile):188 if not os.path.exists(inputfile):
127 if inputfile.endswith('.gz'):189 if inputfile.endswith('.gz'):
128 if not os.path.exists(inputfile[:-3]):190 if not os.path.exists(inputfile[:-3]):
191 misc.sprint(os.getcwd(), os.listdir('.'), inputfile, os.path.exists(inputfile), os.path.exists(inputfile[:-3]))
129 raise self.InvalidCmd('No such file or directory : %s' % inputfile)192 raise self.InvalidCmd('No such file or directory : %s' % inputfile)
130 else: 193 else:
131 inputfile = inputfile[:-3]194 inputfile = inputfile[:-3]
@@ -134,14 +197,22 @@
134 else: 197 else:
135 raise self.InvalidCmd('No such file or directory : %s' % inputfile)198 raise self.InvalidCmd('No such file or directory : %s' % inputfile)
136199
200 self.inputfile = inputfile
201 if self.options['spinmode'] == 'none' and \
202 (self.options['input_format'] not in ['lhe','auto'] or
203 (self.options['input_format'] == 'auto' and '.lhe' not in inputfile[-7:])):
204 self.banner = banner.Banner()
205 self.setup_for_pure_decay()
206 return
207
137 if inputfile.endswith('.gz'):208 if inputfile.endswith('.gz'):
138 misc.gunzip(inputfile)209 misc.gunzip(inputfile)
139 inputfile = inputfile[:-3]210 inputfile = inputfile[:-3]
140
141 # Read the banner of the inputfile211 # Read the banner of the inputfile
142 self.events_file = open(os.path.realpath(inputfile))212 self.events_file = open(os.path.realpath(inputfile))
143 self.banner = banner.Banner(self.events_file)213 self.banner = banner.Banner(self.events_file)
144 214
215
145 # Check the validity of the banner:216 # Check the validity of the banner:
146 if 'slha' not in self.banner:217 if 'slha' not in self.banner:
147 self.events_file = None218 self.events_file = None
@@ -163,7 +234,6 @@
163 self.options['nb_sigma'] = N_sigma234 self.options['nb_sigma'] = N_sigma
164 if self.options['BW_cut'] == -1:235 if self.options['BW_cut'] == -1:
165 self.options['BW_cut'] = float(self.banner.get_detail('run_card', 'bwcutoff'))236 self.options['BW_cut'] = float(self.banner.get_detail('run_card', 'bwcutoff'))
166
167 else:237 else:
168 if not self.options['Nevents_for_max_weigth']:238 if not self.options['Nevents_for_max_weigth']:
169 self.options['Nevents_for_max_weigth'] = 75239 self.options['Nevents_for_max_weigth'] = 75
@@ -273,7 +343,7 @@
273 343
274344
275 #Check the param_card345 #Check the param_card
276 if not bypass_check:346 if not (bypass_check or self.options['input_format'] in ['hepmc', 'lhe_no_banner']):
277 if not hasattr(self.banner, 'param_card'):347 if not hasattr(self.banner, 'param_card'):
278 self.banner.charge_card('slha')348 self.banner.charge_card('slha')
279 param_card = check_param_card.ParamCard(card)349 param_card = check_param_card.ParamCard(card)
@@ -395,42 +465,10 @@
395 465
396 args = self.split_arg(line)466 args = self.split_arg(line)
397 self.check_set(args)467 self.check_set(args)
398 468
399 if args[0] not in ['ms_dir', 'run_card']:469 self.options[args[0]] = ' '.join(args[1:])
400 args[1] = args[1].lower()470
401 471
402 if args[0] in ['max_weight', 'BW_effect','ms_dir', 'spinmode']:
403 self.options[args[0]] = args[1]
404 if args[0] == 'ms_dir':
405 self.options['curr_dir'] = self.options['ms_dir']
406 elif args[0] == 'seed':
407 random.seed(int(args[1]))
408 self.seed = int(args[1])
409 elif args[0] == 'BW_cut':
410 self.options[args[0]] = float(args[1])
411 elif args[0] in ['onlyhelicity', 'use_old_dir']:
412 self.options[args[0]] = banner.ConfigFile.format_variable(args[1], bool, args[0])
413 elif args[0] in ['run_card']:
414 if args[1] == 'default':
415 self.options['run_card'] = None
416 elif os.path.isfile(args[1]):
417 self.options['run_card'] = banner.RunCard(args[1])
418 elif len(args) >2:
419 if not self.options['run_card']:
420 self.options['run_card'] = banner.RunCardLO()
421 self.options['run_card'].remove_all_cut()
422 self.options['run_card'][args[1]] = args[2]
423 elif args[0] in 'fixed_order':
424 if args[1].lower() in ['t', 'true']:
425 self.options['fixed_order'] = True
426 logger.warning('Fix order madspin fails to have the correct scale information. This can bias the results!')
427 logger.warning('Not all functionality of MadSpin handle this mode correctly (only onshell mode so far).')
428 else:
429 self.options['fixed_order'] = False
430 logger.info('fixed_order options set to %s', self.options['fixed_order'])
431 else:
432 self.options[args[0]] = int(args[1])
433
434 def complete_set(self, text, line, begidx, endidx):472 def complete_set(self, text, line, begidx, endidx):
435 473
436474
@@ -559,17 +597,17 @@
559597
560 model_line = self.banner.get('proc_card', 'full_model_line')598 model_line = self.banner.get('proc_card', 'full_model_line')
561599
562 if not self.seed:600 if not self.options['seed']:
563 self.seed = random.randint(0, int(30081*30081))601 self.options['seed'] = random.randint(0, int(30081*30081))
564 self.do_set('seed %s' % self.seed)602 #self.do_set('seed %s' % self.seed)
565 logger.info('Will use seed %s' % self.seed)603 logger.info('Will use seed %s' % self.options['seed'])
566 self.history.insert(0, 'set seed %s' % self.seed)604 self.history.insert(0, 'set seed %s' % self.options['seed'])
567605
568 if self.seed > 30081*30081: # can't use too big random number606 if self.options['seed'] > 30081*30081: # can't use too big random number
569 msg = 'Random seed too large ' + str(self.seed) + ' > 30081*30081'607 msg = 'Random seed too large ' + str(self.options['seed']) + ' > 30081*30081'
570 raise Exception, msg608 raise Exception, msg
571609
572 self.options['seed'] = self.seed610 #self.options['seed'] = self.seed
573 text = '%s\n' % '\n'.join([ line for line in self.history if line])611 text = '%s\n' % '\n'.join([ line for line in self.history if line])
574 self.banner.add_text('madspin' , text)612 self.banner.add_text('madspin' , text)
575 613
@@ -678,9 +716,9 @@
678 716
679 # NOW we have all the information available for RUNNING717 # NOW we have all the information available for RUNNING
680 718
681 if self.seed:719 if self.options['seed']:
682 #seed is specified need to use that one:720 #seed is specified need to use that one:
683 open(pjoin(self.options['ms_dir'],'seeds.dat'),'w').write('%s\n'%self.seed)721 open(pjoin(self.options['ms_dir'],'seeds.dat'),'w').write('%s\n'%self.options['seed'])
684 #remove all ranmar_state722 #remove all ranmar_state
685 for name in misc.glob(pjoin('*', 'SubProcesses','*','ranmar_state.dat'), 723 for name in misc.glob(pjoin('*', 'SubProcesses','*','ranmar_state.dat'),
686 self.options['ms_dir']):724 self.options['ms_dir']):
@@ -741,30 +779,61 @@
741 for name in misc.glob("decay_*_*", self.path_me):779 for name in misc.glob("decay_*_*", self.path_me):
742 shutil.rmtree(name)780 shutil.rmtree(name)
743781
744 self.events_file.close()782 if self.events_file:
745 orig_lhe = lhe_parser.EventFile(self.events_file.name)783 self.events_file.close()
746 784 filename = self.events_file.name
785 else:
786 filename = self.inputfile
787
788 if self.options['input_format'] == 'auto':
789 if '.lhe' in filename :
790 self.options['input_format'] = 'lhe'
791 elif '.hepmc' in filename:
792 self.options['input_format'] = 'hepmc'
793 else:
794 raise Exception, "fail to recognized input format automatically"
795
796 if self.options['input_format'] in ['lhe', 'lhe_no_banner']:
797 orig_lhe = lhe_parser.EventFile(filename)
798 if self.options['input_format'] == 'lhe_no_banner':
799 orig_lhe.allow_empty_event = True
800
801 elif self.options['input_format'] in ['hepmc']:
802 import madgraph.various.hepmc_parser as hepmc_parser
803 orig_lhe = hepmc_parser.HEPMC_EventFile(filename)
804 logger.info("Parsing input event to know how many decay to generate. This can takes few minuts.")
805 else:
806 raise Exception
807
747 to_decay = collections.defaultdict(int)808 to_decay = collections.defaultdict(int)
748 nb_event = 0809 nb_event = 0
810
749 for event in orig_lhe:811 for event in orig_lhe:
750 nb_event +=1812 nb_event +=1
751 for particle in event:813 for particle in event:
752 if particle.status == 1 and particle.pdg in asked_to_decay:814 if particle.status == 1 and particle.pdg in asked_to_decay:
753 # final state and tag as to decay815 # final state and tag as to decay
754 to_decay[particle.pdg] += 1816 to_decay[particle.pdg] += 1
817 if self.options['input_format'] == 'hepmc' and nb_event == 250:
818 currpos = orig_lhe.tell()
819 filesize = orig_lhe.getfilesize()
820 for key in to_decay:
821 to_decay[key] *= 1.05 * filesize/ currpos
822 # 1.05 to avoid accidental coincidence with nevents
823 break
755824
756 # Handle the banner of the output file825 # Handle the banner of the output file
757 if not self.seed:826 if not self.options['seed']:
758 self.seed = random.randint(0, int(30081*30081))827 self.options['seed'] = random.randint(0, int(30081*30081))
759 self.do_set('seed %s' % self.seed)828 #self.do_set('seed %s' % self.seed)
760 logger.info('Will use seed %s' % self.seed)829 logger.info('Will use seed %s' % self.options['seed'])
761 self.history.insert(0, 'set seed %s' % self.seed)830 self.history.insert(0, 'set seed %s' % self.options['seed'])
762831
763 if self.seed > 30081*30081: # can't use too big random number832 if self.options['seed'] > 30081*30081: # can't use too big random number
764 msg = 'Random seed too large ' + str(self.seed) + ' > 30081*30081'833 msg = 'Random seed too large ' + str(self.options['seed']) + ' > 30081*30081'
765 raise Exception, msg834 raise Exception, msg
766835
767 self.options['seed'] = self.seed836 #self.options['seed'] = self.options['seed']
768 837
769 text = '%s\n' % '\n'.join([ line for line in self.history if line])838 text = '%s\n' % '\n'.join([ line for line in self.history if line])
770 self.banner.add_text('madspin' , text)839 self.banner.add_text('madspin' , text)
@@ -780,7 +849,7 @@
780 for pdg, nb_needed in to_decay.items():849 for pdg, nb_needed in to_decay.items():
781 #check if a splitting is needed850 #check if a splitting is needed
782 if nb_needed == nb_event:851 if nb_needed == nb_event:
783 evt_decayfile[pdg] = self.generate_events(pdg, nb_needed, mg5)852 evt_decayfile[pdg] = self.generate_events(pdg, min(nb_needed,100000), mg5)
784 elif nb_needed % nb_event == 0:853 elif nb_needed % nb_event == 0:
785 nb_mult = nb_needed // nb_event854 nb_mult = nb_needed // nb_event
786 part = self.model.get_particle(pdg)855 part = self.model.get_particle(pdg)
@@ -788,64 +857,107 @@
788 if name not in self.list_branches:857 if name not in self.list_branches:
789 continue858 continue
790 elif len(self.list_branches[name]) == nb_mult:859 elif len(self.list_branches[name]) == nb_mult:
791 evt_decayfile[pdg] = self.generate_events(pdg, nb_event, mg5)860 evt_decayfile[pdg] = self.generate_events(pdg, min(nb_event,100000), mg5)
792 else:861 else:
793 evt_decayfile[pdg] = self.generate_events(pdg, nb_needed, mg5, cumul=True)862 evt_decayfile[pdg] = self.generate_events(pdg, min(nb_needed,100000), mg5, cumul=True)
863 elif self.options['cross_section']:
864 #cross-section hard-coded -> allow
865 part = self.model.get_particle(pdg)
866 name = part.get_name()
867
868 if name not in self.list_branches:
869 continue
870 else:
871 try:
872 evt_decayfile[pdg] = self.generate_events(pdg, min(nb_needed,100000), mg5, cumul=True)
873 except common_run_interface.ZeroResult:
874 logger.warning("Branching ratio is zero for this particle. Not decaying it")
875 del to_decay[pdg]
794 else:876 else:
795 part = self.model.get_particle(pdg)877 part = self.model.get_particle(pdg)
796 name = part.get_name()878 name = part.get_name()
797 if name not in self.list_branches or len(self.list_branches[name]) == 0:879 if name not in self.list_branches or len(self.list_branches[name]) == 0:
798 continue880 continue
799 raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")881 raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed. One workaround is to force the final cross-section manually.")
800 882
801 883
802
803
804 # Compute the branching ratio.884 # Compute the branching ratio.
805 br = 1885 if not self.options['cross_section']:
806 for (pdg, event_files) in evt_decayfile.items():886 br = 1
807 if not event_files:887 for (pdg, event_files) in evt_decayfile.items():
808 continue888 if not event_files:
809 totwidth = float(self.banner.get('param', 'decay', abs(pdg)).value)889 continue
810 if to_decay[pdg] == nb_event:890 totwidth = float(self.banner.get('param', 'decay', abs(pdg)).value)
811 # Exactly one particle of this type to decay by event891 if to_decay[pdg] == nb_event:
812 pwidth = sum([event_files[k].cross for k in event_files])892 # Exactly one particle of this type to decay by event
813 if pwidth > 1.01 * totwidth:893 pwidth = sum([event_files[k].cross for k in event_files])
814 logger.critical("Branching ratio larger than one for %s " % pdg)
815 br *= pwidth / totwidth
816 elif to_decay[pdg] % nb_event == 0:
817 # More than one particle of this type to decay by event
818 # Need to check the number of event file to check if we have to
819 # make separate type of decay or not.
820 nb_mult = to_decay[pdg] // nb_event
821 if nb_mult == len(event_files):
822 for k in event_files:
823 pwidth = event_files[k].cross
824 if pwidth > 1.01 * totwidth:
825 logger.critical("Branching ratio larger than one for %s " % pdg)
826 br *= pwidth / totwidth
827 br *= math.factorial(nb_mult)
828 else:
829 pwidth = sum(event_files[k].cross for k in event_files)
830 if pwidth > 1.01 * totwidth:894 if pwidth > 1.01 * totwidth:
831 logger.critical("Branching ratio larger than one for %s " % pdg) 895 logger.critical("Branching ratio larger than one for %s " % pdg)
832 br *= (pwidth / totwidth)**nb_mult896 br *= pwidth / totwidth
833 else:897 elif to_decay[pdg] % nb_event == 0:
834 raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")898 # More than one particle of this type to decay by event
899 # Need to check the number of event file to check if we have to
900 # make separate type of decay or not.
901 nb_mult = to_decay[pdg] // nb_event
902 if nb_mult == len(event_files):
903 for k in event_files:
904 pwidth = event_files[k].cross
905 if pwidth > 1.01 * totwidth:
906 logger.critical("Branching ratio larger than one for %s " % pdg)
907 br *= pwidth / totwidth
908 br *= math.factorial(nb_mult)
909 else:
910 pwidth = sum(event_files[k].cross for k in event_files)
911 if pwidth > 1.01 * totwidth:
912 logger.critical("Branching ratio larger than one for %s " % pdg)
913 br *= (pwidth / totwidth)**nb_mult
914 else:
915 raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
916 else:
917 br = 1
835 self.branching_ratio = br918 self.branching_ratio = br
836 self.efficiency = 1919 self.efficiency = 1
837 self.cross, self.error = self.banner.get_cross(witherror=True)920 try:
921 self.cross, self.error = self.banner.get_cross(witherror=True)
922 except:
923 if self.options['input_format'] != 'lhe':
924 self.cross, self.error = 0, 0
838 self.cross *= br925 self.cross *= br
839 self.error *= br926 self.error *= br
840 927
841
842 # modify the cross-section in the init block of the banner928 # modify the cross-section in the init block of the banner
843 self.banner.scale_init_cross(self.branching_ratio)929 if not self.options['cross_section']:
844 930 self.banner.scale_init_cross(self.branching_ratio)
845 931 else:
932
933 if self.options['input_format'] in ['lhe_no_banner','hepmc'] and 'init' not in self.banner:
934 self.cross = sum(self.options['cross_section'].values())
935 self.error = 0
936 self.branching_ratio = 1
937 else:
938 self.banner.modify_init_cross(self.options['cross_section'])
939 new_cross, new_error = self.banner.get_cross(witherror=True)
940 self.branching_ratio = new_cross / self.cross
941 self.cross = new_cross
942 self.error = new_error
943
846 # 3. Merge the various file together.944 # 3. Merge the various file together.
847 output_lhe = lhe_parser.EventFile(orig_lhe.name.replace('.lhe', '_decayed.lhe.gz'), 'w')945 if self.options['input_format'] == 'hepmc':
848 self.banner.write(output_lhe, close_tag=False)946 name = orig_lhe.name.replace('.hepmc', '_decayed.lhe')
947 if not name.endswith('.gz'):
948 name = '%s.gz' % name
949
950 output_lhe = lhe_parser.EventFile(name, 'w')
951 else:
952 name = orig_lhe.name.replace('.lhe', '_decayed.lhe')
953 if not name.endswith('.gz'):
954 name = '%s.gz' % name
955 output_lhe = lhe_parser.EventFile(name, 'w')
956 try:
957 self.banner.write(output_lhe, close_tag=False)
958 except Exception:
959 if self.options['input_format'] == 'lhe':
960 raise
849 961
850 # initialise object which store not use event due to wrong helicity962 # initialise object which store not use event due to wrong helicity
851 bufferedEvents_decay = {}963 bufferedEvents_decay = {}
@@ -856,8 +968,9 @@
856 start = time.time()968 start = time.time()
857 counter = 0969 counter = 0
858 orig_lhe.seek(0)970 orig_lhe.seek(0)
971
859 for event in orig_lhe:972 for event in orig_lhe:
860 if counter and counter % 1000 == 0 and float(str(counter)[1:]) ==0:973 if counter and counter % 100 == 0 and float(str(counter)[1:]) ==0:
861 print "decaying event number %s [%s s]" % (counter, time.time()-start)974 print "decaying event number %s [%s s]" % (counter, time.time()-start)
862 counter +=1975 counter +=1
863 976
@@ -866,10 +979,17 @@
866 particles = [p for p in event if int(p.status) == 1.0]979 particles = [p for p in event if int(p.status) == 1.0]
867 random.shuffle(particles)980 random.shuffle(particles)
868 ids = [particle.pid for particle in particles]981 ids = [particle.pid for particle in particles]
982 br = 1 #br for that particular events (for special/weighted case)
983 hepmc_output = lhe_parser.Event() #for hepmc case: collect the decay particle
869 for i,particle in enumerate(particles):984 for i,particle in enumerate(particles):
985 #misc.sprint(i, particle.pdg, particle.pid)
986 #misc.sprint(self.final_state, evt_decayfile)
870 # check if we need to decay the particle 987 # check if we need to decay the particle
871 if particle.pdg not in self.final_state or particle.pdg not in evt_decayfile:988 if self.final_state and particle.pdg not in self.final_state:
872 continue # nothing to do for this particle989 continue # nothing to do for this particle
990 if particle.pdg not in evt_decayfile:
991 continue # nothing to do for this particle
992
873 # check how the decay need to be done993 # check how the decay need to be done
874 nb_decay = len(evt_decayfile[particle.pdg])994 nb_decay = len(evt_decayfile[particle.pdg])
875 if nb_decay == 0:995 if nb_decay == 0:
@@ -888,14 +1008,28 @@
888 cumul = 01008 cumul = 0
889 for j,events in evt_decayfile[particle.pdg].items():1009 for j,events in evt_decayfile[particle.pdg].items():
890 cumul += events.cross1010 cumul += events.cross
891 if r < cumul:1011 if r <= cumul:
892 decay_file = events1012 decay_file = events
893 decay_file_nb = j1013 decay_file_nb = j
894 else:
895 break1014 break
896 1015 else:
1016 # security for numerical accuracy issue... (unlikely but better safe)
1017 if (cumul-tot)/tot < 1e-5:
1018 decay_file = events
1019 decay_file_nb = j
1020 else:
1021 misc.sprint(j,cumul, events.cross, tot, (tot-cumul)/tot)
1022 raise Exception
1023
1024 if self.options['new_wgt'] == 'BR':
1025 tot_width = float(self.banner.get('param', 'decay', abs(pdg)).value)
1026 if tot_width:
1027 br = decay_file.cross / tot_width
897 # ok start the procedure1028 # ok start the procedure
898 helicity = particle.helicity1029 if hasattr(particle,'helicity'):
1030 helicity = particle.helicity
1031 else:
1032 helicity = 9
899 bufferedEvents = bufferedEvents_decay[particle.pdg][decay_file_nb]1033 bufferedEvents = bufferedEvents_decay[particle.pdg][decay_file_nb]
900 1034
901 # now that we have the file to read. find the associate event1035 # now that we have the file to read. find the associate event
@@ -910,8 +1044,8 @@
910 except StopIteration:1044 except StopIteration:
911 # check how far we are1045 # check how far we are
912 ratio = counter / nb_event 1046 ratio = counter / nb_event
913 needed = 1.05 * to_decay[particle.pdg]/ratio - counter1047 needed = 1.05 * to_decay[particle.pdg] - counter
914 needed = min(1000, max(needed, 1000))1048 needed = min(100000, max(needed, 6000))
915 with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):1049 with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
916 new_file = self.generate_events(particle.pdg, needed, mg5, [decay_file_nb])1050 new_file = self.generate_events(particle.pdg, needed, mg5, [decay_file_nb])
917 evt_decayfile[particle.pdg].update(new_file)1051 evt_decayfile[particle.pdg].update(new_file)
@@ -927,14 +1061,36 @@
927 # only add to the buffering if the buffer is not too large1061 # only add to the buffering if the buffer is not too large
928 bufferedEvents[helicity].append(decay)1062 bufferedEvents[helicity].append(decay)
929 # now that we have the event make the merge1063 # now that we have the event make the merge
930 particle.add_decay(decay)1064 if self.options['input_format'] != 'hepmc':
1065 particle.add_decay(decay)
1066 else:
1067 if len(hepmc_output) == 0:
1068 hepmc_output.append(lhe_parser.Particle(event=hepmc_output))
1069 hepmc_output[0].color2 = 0
1070 decayed_particle = lhe_parser.Particle(particle, hepmc_output)
1071 hepmc_output.append(decayed_particle)
1072 decayed_particle.add_decay(decay)
931 # change the weight associate to the event1073 # change the weight associate to the event
932 event.wgt *= self.branching_ratio1074 if self.options['new_wgt'] == 'cross-section':
933 wgts = event.parse_reweight()1075 event.wgt *= self.branching_ratio
934 for key in wgts:1076 br = self.branching_ratio
935 wgts[key] *= self.branching_ratio1077 else:
936 # all particle have been decay if needed1078 event.wgt *= br
937 output_lhe.write(str(event))1079
1080 if self.options['input_format'] != 'hepmc':
1081 wgts = event.parse_reweight()
1082 for key in wgts:
1083 wgts[key] *= br
1084 # all particle have been decay if needed
1085 output_lhe.write(str(event))
1086 else:
1087 hepmc_output.wgt = event.wgt
1088 hepmc_output.nexternal = len(hepmc_output) # the append does not update nexternal
1089 hepmc_output.assign_mother()
1090 output_lhe.write(str(hepmc_output))
1091 else:
1092 if counter==0:
1093 raise Exception
938 output_lhe.write('</LesHouchesEvents>\n') 1094 output_lhe.write('</LesHouchesEvents>\n')
939 1095
940 1096
@@ -971,7 +1127,12 @@
971 cumul allow to merge all the definition in one run (add process)1127 cumul allow to merge all the definition in one run (add process)
972 to generate events according to cross-section1128 to generate events according to cross-section
973 """1129 """
974 1130 if not hasattr(self, 'me_int'):
1131 self.me_int = {}
1132
1133
1134
1135 nb_event = int(nb_event) # in case of hepmc request the nb_event is not an integer
975 if cumul:1136 if cumul:
976 width = 0.1137 width = 0.
977 else: 1138 else:
@@ -1003,31 +1164,36 @@
1003 options = dict(mg5.options)1164 options = dict(mg5.options)
1004 if self.options['ms_dir']:1165 if self.options['ms_dir']:
1005 # we are in gridpack mode -> create it1166 # we are in gridpack mode -> create it
1006 me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\1167 if decay_dir in self.me_int:
1007 decay_dir), options=options)1168 me5_cmd = self.me_int[decay_dir]
1008 me5_cmd.options["automatic_html_opening"] = False1169 else:
1009 me5_cmd.options["madanalysis5_path"] = None1170 me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
1010 me5_cmd.options["madanalysis_path"] = None1171 decay_dir), options=options)
1011 me5_cmd.allow_notification_center = False1172 me5_cmd.options["automatic_html_opening"] = False
1012 try:1173 me5_cmd.options["madanalysis5_path"] = None
1013 os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))1174 me5_cmd.options["madanalysis_path"] = None
1014 os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))1175 me5_cmd.allow_notification_center = False
1015 except Exception,error:1176 try:
1016 logger.debug(error)1177 os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
1017 pass 1178 os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
1179 except Exception,error:
1180 logger.debug(error)
1181 pass
1182 self.me_int[decay_dir] = me5_cmd
1183
1018 if self.options["run_card"]:1184 if self.options["run_card"]:
1019 run_card = self.options["run_card"]1185 run_card = self.run_card
1020 else:1186 else:
1021 run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat")) 1187 run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))
1022 1188
1023 run_card["iseed"] = self.seed1189 run_card["iseed"] = self.options['seed']
1024 run_card['gridpack'] = True1190 run_card['gridpack'] = True
1025 run_card['systematics_program'] = 'False'1191 run_card['systematics_program'] = 'False'
1026 run_card['use_syst'] = False1192 run_card['use_syst'] = False
1027 run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))1193 run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
1028 param_card = self.banner['slha']1194 param_card = self.banner['slha']
1029 open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)1195 open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
1030 self.seed += 11196 self.options['seed'] += 1
1031 # actually creation1197 # actually creation
1032 me5_cmd.exec_cmd("generate_events run_01 -f")1198 me5_cmd.exec_cmd("generate_events run_01 -f")
1033 if output_width:1199 if output_width:
@@ -1041,26 +1207,35 @@
1041 misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir)1207 misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir)
1042 1208
1043 # Now generate the events1209 # Now generate the events
1044
1045 if not self.options['ms_dir']:1210 if not self.options['ms_dir']:
1046 me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\1211 if decay_dir in self.me_int:
1047 decay_dir), options=mg5.options)1212 me5_cmd = self.me_int[decay_dir]
1048 me5_cmd.options["automatic_html_opening"] = False1213 else:
1049 me5_cmd.options["automatic_html_opening"] = False1214 me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
1050 me5_cmd.options["madanalysis5_path"] = None1215 decay_dir), options=mg5.options)
1051 me5_cmd.options["madanalysis_path"] = None1216 me5_cmd.options["automatic_html_opening"] = False
1052 me5_cmd.allow_notification_center = False1217 me5_cmd.options["automatic_html_opening"] = False
1053 try:1218 me5_cmd.options["madanalysis5_path"] = None
1054 os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))1219 me5_cmd.options["madanalysis_path"] = None
1055 os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))1220 me5_cmd.allow_notification_center = False
1056 except Exception,error:1221 try:
1057 logger.debug(error)1222 os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
1058 pass 1223 os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
1224 except Exception,error:
1225 logger.debug(error)
1226 pass
1227 self.me_int[decay_dir] = me5_cmd
1059 if self.options["run_card"]:1228 if self.options["run_card"]:
1060 run_card = self.options["run_card"]1229 run_card = self.run_card
1061 else:1230 else:
1062 run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))1231 run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))
1063 run_card["nevents"] = int(1.2*nb_event)1232 run_card["nevents"] = int(1.2*nb_event)
1233 # Handle the banner of the output file
1234 if not self.seed:
1235 self.seed = random.randint(0, int(30081*30081))
1236 self.do_set('seed %s' % self.seed)
1237 logger.info('Will use seed %s' % self.seed)
1238 self.history.insert(0, 'set seed %s' % self.seed)
1064 run_card["iseed"] = self.seed1239 run_card["iseed"] = self.seed
1065 run_card["systematics_program"] = 'None'1240 run_card["systematics_program"] = 'None'
1066 run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))1241 run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
@@ -1140,8 +1315,7 @@
1140 if particle.status == 1 and particle.pdg in asked_to_decay:1315 if particle.status == 1 and particle.pdg in asked_to_decay:
1141 # final state and tag as to decay1316 # final state and tag as to decay
1142 to_decay[particle.pdg] += 11317 to_decay[particle.pdg] += 1
1143 #misc.sprint(to_decay)1318
1144 #misc.sprint("import the mode -> temporary with logging")
1145 with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):1319 with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
1146 mg5 = self.mg5cmd1320 mg5 = self.mg5cmd
1147 if not self.model:1321 if not self.model:
@@ -1436,10 +1610,6 @@
1436 full_event = lhe_parser.Event(str(production))1610 full_event = lhe_parser.Event(str(production))
1437 full_event = full_event.add_decays(decays)1611 full_event = full_event.add_decays(decays)
1438 full_me = self.calculate_matrix_element(full_event)1612 full_me = self.calculate_matrix_element(full_event)
1439 #misc.sprint(full_event)
1440 #misc.sprint([p.pdg for p in production])
1441 #misc.sprint([p.pdg for p in full_event])
1442 #misc.sprint(full_me, production_me, decay_me)
1443 return full_event, full_me/(production_me*decay_me)1613 return full_event, full_me/(production_me*decay_me)
1444 1614
1445 1615
@@ -1457,7 +1627,6 @@
1457 init = (-init[0],)1627 init = (-init[0],)
1458 final = tuple(-i for i in final)1628 final = tuple(-i for i in final)
1459 tag = (init, final)1629 tag = (init, final)
1460 misc.sprint([k for k in self.all_me.keys() if len(k[0])==1])
1461 orig_order = self.all_me[tag]['order']1630 orig_order = self.all_me[tag]['order']
1462 1631
1463 pdir = self.all_me[tag]['pdir']1632 pdir = self.all_me[tag]['pdir']
@@ -1520,7 +1689,6 @@
1520 1689
1521 rwgt_interface.ReweightInterface.get_LO_definition_from_NLO() 1690 rwgt_interface.ReweightInterface.get_LO_definition_from_NLO()
1522 1691
1523
15241692
15251693
1526if __name__ == '__main__':1694if __name__ == '__main__':
15271695
=== modified file 'MadSpin/madspin'
--- MadSpin/madspin 2013-11-12 23:52:05 +0000
+++ MadSpin/madspin 2018-12-20 08:18:55 +0000
@@ -49,6 +49,7 @@
49if len(args) == 0:49if len(args) == 0:
50 args = ''50 args = ''
5151
52
52import subprocess53import subprocess
5354
54# Check if optimize mode is (and should be) activated55# Check if optimize mode is (and should be) activated
@@ -114,23 +115,6 @@
114 pass115 pass
115import MadSpin.interface_madspin as interface116import MadSpin.interface_madspin as interface
116117
117try:
118 cmd_line = interface.MadSpinInterface()
119 cmd_line.use_rawinput = True
120 cmd_line.cmdloop()
121except:
122 pass
123try:
124 cmd_line.exec_cmd('quit all', printcmd=False)
125 readline.set_history_length(100)
126 if not os.path.exists(os.path.join(os.environ['HOME'], '.mg5')):
127 os.mkdir(os.path.join(os.environ['HOME'], '.mg5'))
128 readline.write_history_file(history_file)
129except Exception, error:
130 pass
131
132sys.exit(0)
133
134# Call the cmd interface main loop118# Call the cmd interface main loop
135try:119try:
136 if options.file or args:120 if options.file or args:
@@ -143,12 +127,14 @@
143 cmd_line = interface.MasterCmdWeb()127 cmd_line = interface.MasterCmdWeb()
144 cmd_line.debug_output = os.path.join(os.path.dirname(input_file),'generation.log')128 cmd_line.debug_output = os.path.join(os.path.dirname(input_file),'generation.log')
145 cmd_line.use_rawinput = False129 cmd_line.use_rawinput = False
146 cmd_line.run_cmd('import ' + input_file)130 cmd_line.import_command_file(input_file)
147 cmd_line.run_cmd('quit')131 cmd_line.run_cmd('quit')
148 else:132 else:
133 print "using input+file", input_file
149 cmd_line = interface.MadSpinInterface()134 cmd_line = interface.MadSpinInterface()
150 cmd_line.use_rawinput = False135 cmd_line.use_rawinput = False
151 cmd_line.run_cmd('import ' + input_file)136 cmd_line.haspiping = False
137 cmd_line.import_command_file(input_file)
152 cmd_line.run_cmd('quit')138 cmd_line.run_cmd('quit')
153 else:139 else:
154 # Interactive mode140 # Interactive mode
@@ -160,8 +146,20 @@
160 cmd_line = interface.MasterCmdWeb()146 cmd_line = interface.MasterCmdWeb()
161 cmd_line.cmdloop()147 cmd_line.cmdloop()
162 else:148 else:
163 cmd_line = interface.MasterCmd(mgme_dir = options.mgme_dir)149 try:
164 cmd_line.cmdloop()150 cmd_line = interface.MadSpinInterface()
151 cmd_line.use_rawinput = True
152 cmd_line.cmdloop()
153 except:
154 pass
155 try:
156 cmd_line.exec_cmd('quit all', printcmd=False)
157 readline.set_history_length(100)
158 if not os.path.exists(os.path.join(os.environ['HOME'], '.mg5')):
159 os.mkdir(os.path.join(os.environ['HOME'], '.mg5'))
160 readline.write_history_file(history_file)
161 except Exception, error:
162 pass
165except KeyboardInterrupt:163except KeyboardInterrupt:
166 print 'writting history and quit on KeyboardInterrupt' 164 print 'writting history and quit on KeyboardInterrupt'
167 pass165 pass
168166
=== modified file 'MadSpin/src/driver.f'
--- MadSpin/src/driver.f 2017-08-03 18:49:25 +0000
+++ MadSpin/src/driver.f 2018-12-20 08:18:55 +0000
@@ -340,10 +340,10 @@
340 if (jac.lt.0d0) then340 if (jac.lt.0d0) then
341 counter2=counter2+1 341 counter2=counter2+1
342 counter3=counter3+1 342 counter3=counter3+1
343c if (counter3.gt.500) then343 if (counter3.gt.500*8*100) then
344c write(*,*) "500_pts_failed_stop_executation"344 write(*,*) "500_pts_failed_stop_executation"
345c stop345 stop
346c endif346 endif
347 if (counter2.ge.8) then ! use another topology to generate PS points347 if (counter2.ge.8) then ! use another topology to generate PS points
348 do k=1,n_max_cg348 do k=1,n_max_cg
349 amp2(k)=0d0349 amp2(k)=0d0
@@ -960,7 +960,7 @@
960 include 'nexternal_prod.inc'960 include 'nexternal_prod.inc'
961 !include 'run.inc'961 !include 'run.inc'
962c arguments962c arguments
963 double precision jac,x(36),p(0:3,nexternal), p_prod(0:3,nexternal)963 double precision jac,x(36),p(0:3,nexternal), p_prod(0:3,nexternal_prod)
964 integer itree(2,-nexternal:-1)964 integer itree(2,-nexternal:-1)
965 double precision qmass(-nexternal:0),qwidth(-nexternal:0)965 double precision qmass(-nexternal:0),qwidth(-nexternal:0)
966c common966c common
967967
=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat 2018-04-25 21:20:23 +0000
+++ Template/LO/Cards/run_card.dat 2018-12-20 08:18:55 +0000
@@ -12,7 +12,7 @@
12# mind the format: value = variable ! comment *12# mind the format: value = variable ! comment *
13# *13# *
14# To display more options, you can type the command: *14# To display more options, you can type the command: *
15* update full_run_card *15# update full_run_card *
16#*********************************************************************16#*********************************************************************
17#17#
18#******************* 18#*******************
1919
=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f 2018-04-20 08:38:20 +0000
+++ Template/LO/SubProcesses/genps.f 2018-12-20 08:18:55 +0000
@@ -41,7 +41,7 @@
41 integer mincfig,maxcfig !Range of configurations41 integer mincfig,maxcfig !Range of configurations
42 integer invar42 integer invar
43 double precision wgt !(input and output)43 double precision wgt !(input and output)
44 double precision x(maxdim),p(maxdim) !x,p (output) [p(0:3,nexternal)]44 double precision x(*),p(*) !x,p (output) [p(0:3,nexternal)]
45c45c
46c Local46c Local
47c47c
@@ -104,7 +104,7 @@
104c104c
105 integer iconfig,mincfig,maxcfig,invar105 integer iconfig,mincfig,maxcfig,invar
106 double precision p1(0:3,nexternal+1)106 double precision p1(0:3,nexternal+1)
107 double precision x(maxinvar)107 double precision x(*)
108 double precision wgt108 double precision wgt
109c109c
110c Local110c Local
@@ -191,8 +191,6 @@
191c191c
192 include 'configs.inc'192 include 'configs.inc'
193 data firsttime/.true./193 data firsttime/.true./
194 integer isym(0:100)
195c data isym /2,1,5,27,42,47,0,0,0,0,0/
196 data jfig/1/194 data jfig/1/
197c-----195c-----
198c Begin Code196c Begin Code
@@ -215,11 +213,6 @@
215 write(*,'(a,12e10.3)') ' Masses:',(m(i),i=1,nparticles)213 write(*,'(a,12e10.3)') ' Masses:',(m(i),i=1,nparticles)
216 endif !First_time214 endif !First_time
217215
218 if (.false.) then
219 iconfig = isym(jfig)
220 jfig = jfig+1
221 if (jfig .gt. isym(0)) jfig=1
222 endif
223 this_config = iconfig !Pass iconfig to amplitude routine216 this_config = iconfig !Pass iconfig to amplitude routine
224C217C
225C Get fraction of beam energy if pdf's are used218C Get fraction of beam energy if pdf's are used
@@ -279,19 +272,18 @@
279 cm_rap=.5d0*dlog((p0+p3)/(p0-p3))272 cm_rap=.5d0*dlog((p0+p3)/(p0-p3))
280 set_cm_rap=.true.273 set_cm_rap=.true.
281c Set shat274c Set shat
282 s(-nbranch) = m2**2+2*xbk(1)*ebeam(1) *275 s(-nbranch) = x(ndim)*stot
283 $ (ebeam(2)+sqrt(ebeam(2)**2-m2**2))
284 elseif (abs(lpp(2)) .ge. 1) then276 elseif (abs(lpp(2)) .ge. 1) then
285 call sample_get_x(sjac,x(ndim),ndim,mincfig,0d0,1d0)277 call sample_get_x(sjac,x(ndim),ndim,mincfig,0d0,1d0)
286 xbk(2) = x(ndim)278 xbk(2) = x(ndim)
279
287c Set CM rapidity for use in the rap() function280c Set CM rapidity for use in the rap() function
288 p0=ebeam(1)+xbk(2)*ebeam(2)281 p0=ebeam(1)+xbk(2)*ebeam(2)
289 p3=sqrt(ebeam(1)**2-m1**2)-xbk(2)*ebeam(2)282 p3=sqrt(ebeam(1)**2-m1**2)-xbk(2)*ebeam(2)
290 cm_rap=.5d0*dlog((p0+p3)/(p0-p3))283 cm_rap=.5d0*dlog((p0+p3)/(p0-p3))
291 set_cm_rap=.true.284 set_cm_rap=.true.
292c Set shat285c Set shat
293 s(-nbranch) = m1**2+2*(ebeam(1)+sqrt(ebeam(1)**2-m1**2))286 s(-nbranch) = x(ndim)*stot
294 $ * xbk(2)*ebeam(2)
295 else287 else
296c Set CM rapidity for use in the rap() function288c Set CM rapidity for use in the rap() function
297 p0=ebeam(1) + ebeam(2)289 p0=ebeam(1) + ebeam(2)
@@ -576,7 +568,7 @@
576 if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0568 if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
577 if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0569 if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
578 if (mass_ion(1).ge.0d0) m1 = mass_ion(1)570 if (mass_ion(1).ge.0d0) m1 = mass_ion(1)
579 if (mass_ion(2).ge.0d0) m1 = mass_ion(2)571 if (mass_ion(2).ge.0d0) m2 = mass_ion(2)
580 if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1572 if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1
581 if(ebeam(2).lt.m2.and.lpp(2).ne.9) ebeam(2)=m2573 if(ebeam(2).lt.m2.and.lpp(2).ne.9) ebeam(2)=m2
582 pi1(0)=ebeam(1)574 pi1(0)=ebeam(1)
@@ -598,7 +590,6 @@
598 maxcfig=iconfig590 maxcfig=iconfig
599 call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)591 call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
600 maxwgt=0d0592 maxwgt=0d0
601c write(*,'(a,12i4)') 'Summing configs',(isym(i),i=1,isym(0))
602 nparticles = nexternal593 nparticles = nexternal
603 nfinal = nparticles-nincoming594 nfinal = nparticles-nincoming
604 nbranch = nparticles-2595 nbranch = nparticles-2
@@ -617,8 +608,6 @@
617c608c
618c do i=1,mapconfig(0)609c do i=1,mapconfig(0)
619 do i=mincfig,maxcfig610 do i=mincfig,maxcfig
620c do k=1,isym(0)
621c i = isym(k)
622 write(*,'(15i4)') i,(minvar(j,i),j=1,ndim)611 write(*,'(15i4)') i,(minvar(j,i),j=1,ndim)
623 do j=1,ndim612 do j=1,ndim
624 ipole = minvar(j,i)613 ipole = minvar(j,i)
@@ -663,7 +652,7 @@
663c double precision spole(-max_branch:0),swidth(-max_branch:0)652c double precision spole(-max_branch:0),swidth(-max_branch:0)
664 double precision jac,pswgt653 double precision jac,pswgt
665 integer nbranch654 integer nbranch
666 double precision x(40) ! ja 3/2/11 21->40 after strange segfault655 double precision x(*) ! ja 3/2/11 21->40 after strange segfault
667c656c
668c Local657c Local
669c658c
670659
=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f 2018-03-11 12:29:51 +0000
+++ Template/LO/SubProcesses/reweight.f 2018-12-20 08:18:55 +0000
@@ -1,3 +1,21 @@
1c for cross-checking change in this file.
2c here is a minimal list of process that we have to
3c test
4c
5c SM
6c ----
7c p p > t t~ (up to 2jet)
8c p p > w+ (up to 3 jet)
9c p p > j j w+ (ordering seems important)
10c p p > z t t~ j j (no MLM needed)
11c
12c
13c HEFT
14c ----
15c p p > h j b b~
16c q q > a a g q q
17c g g > h g q q
18
1 double precision function gamma(q0)19 double precision function gamma(q0)
2c**************************************************20c**************************************************
3c calculates the branching probability21c calculates the branching probability
@@ -226,13 +244,15 @@
226 integer iddgluon, iddother, idgluon, idother244 integer iddgluon, iddother, idgluon, idother
227 logical isqcd245 logical isqcd
228 external isqcd246 external isqcd
247 integer get_color
248 external get_color
229249
230 idmo=ipdg(imo)250 idmo=ipdg(imo)
231 idda1=ipdg(ida1)251 idda1=ipdg(ida1)
232 idda2=ipdg(ida2)252 idda2=ipdg(ida2)
233253
234 if (btest(mlevel,4)) then254 if (btest(mlevel,4)) then
235 write(*,*) ' updating ipart for: ',ida1,ida2,' -> ',imo255 write(*,*) 'updating ipart for: ',ida1,ida2,' -> ',imo
236 endif256 endif
237257
238 if (btest(mlevel,4)) then258 if (btest(mlevel,4)) then
@@ -274,7 +294,6 @@
274 $ ' (',ipdg(imo),')'294 $ ' (',ipdg(imo),')'
275 return295 return
276 endif 296 endif
277
278c FS clustering297c FS clustering
279c Transmit parton PDG code for parton vertex298c Transmit parton PDG code for parton vertex
280 if(isjet(idmo)) then299 if(isjet(idmo)) then
@@ -353,14 +372,26 @@
353c quark -> gluon-quark or Z-quark or h-quark or W-quark372c quark -> gluon-quark or Z-quark or h-quark or W-quark
354 ipart(1,imo)=ipart(1,ida2)373 ipart(1,imo)=ipart(1,ida2)
355 ipart(2,imo)=0374 ipart(2,imo)=0
356 else375 else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda1)).eq.3.and.get_color(idda2).eq.1) then
376c exotic q > q' Scalar
377 ipart(1,imo)=ipart(1,ida1)
378 ipart(2,imo)=0
379 else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda2)).eq.3.and.get_color(idda1).eq.1) then
380c exotic q > Scalar q'
381 ipart(1,imo)=ipart(1,ida2)
382 ipart(2,imo)=0
383 else if (get_color(idmo).eq.1) then
357c Color singlet384c Color singlet
358 ipart(1,imo)=ipart(1,ida1)385 ipart(1,imo)=ipart(1,ida1)
359 ipart(2,imo)=ipart(1,ida2)386 ipart(2,imo)=ipart(1,ida2)
387 else
388 write(*,*) idmo,'>', idda1, idda2, 'color', get_color(idmo),'>', get_color(idda1), get_color(idda2)
389 write(*,*) "failed for ipartupdate. Please retry without MLM/default dynamical scale"
390 stop 3
360 endif391 endif
361 392
362 if (btest(mlevel,4)) then393 if (btest(mlevel,4)) then
363 write(*,*) ' -> ',(ipart(i,imo),i=1,2),' (',ipdg(imo),')'394 write(*,*) 'XY -> ',(ipart(i,imo),i=1,2),' (',ipdg(imo),')'
364 endif395 endif
365396
366 return397 return
@@ -522,7 +553,8 @@
522 common/to_stot/stot,m1,m2553 common/to_stot/stot,m1,m2
523554
524C local variables555C local variables
525 integer i, j, idi, idj, k556 integer i, j, idi, idj, k,m
557 integer get_color
526 real*8 PI558 real*8 PI
527 parameter( PI = 3.14159265358979323846d0 )559 parameter( PI = 3.14159265358979323846d0 )
528 integer iforest(2,-max_branch:-1,lmaxconfigs)560 integer iforest(2,-max_branch:-1,lmaxconfigs)
@@ -545,6 +577,7 @@
545 logical chclusold,fail,increasecode577 logical chclusold,fail,increasecode
546 save chclusold578 save chclusold
547 integer tmpindex579 integer tmpindex
580 integer pdgm, pdgid1, pdgid2
548581
549 logical isqcd,isjet,isparton,cluster,isjetvx,is_octet582 logical isqcd,isjet,isparton,cluster,isjetvx,is_octet
550 integer ifsno583 integer ifsno
@@ -676,6 +709,7 @@
676c increasecode gives whether we should increase jcode at next vertex709c increasecode gives whether we should increase jcode at next vertex
677 increasecode=.false.710 increasecode=.false.
678 do n=1,nexternal-2711 do n=1,nexternal-2
712c write(*,*) 'QCD jet status (before n= ',n,'):',(iqjets(i),i=3,nexternal)
679 do i=1,2 ! index of the child in the interaction713 do i=1,2 ! index of the child in the interaction
680 do j=1,2 ! j index of the beam714 do j=1,2 ! j index of the beam
681 if(idacl(n,i).eq.ibeam(j))then715 if(idacl(n,i).eq.ibeam(j))then
@@ -746,8 +780,16 @@
746c The ishft gives the FS particle corresponding to imocl780c The ishft gives the FS particle corresponding to imocl
747 if(.not.is_octet(ipdgcl(ishft(1,ipart(1,imocl(n))-1),igraphs(1),iproc)))then781 if(.not.is_octet(ipdgcl(ishft(1,ipart(1,imocl(n))-1),igraphs(1),iproc)))then
748 ! split case for q a > q and for g > g h (with the gluon splitting into quark)782 ! split case for q a > q and for g > g h (with the gluon splitting into quark)
749 if (ipart(2,imocl(n)).eq.0) then ! q a > q case783 ! also check for case of three scalar interaction (then do nothing)
750 iqjets(ipart(1,imocl(n)))=0784 pdgm = ipdgcl(imocl(n),igraphs(1),iproc)
785 pdgid1 = ipdgcl(idacl(n,1),igraphs(1),iproc)
786 pdgid2 = ipdgcl(idacl(n,2),igraphs(1),iproc)
787 if (.not.isqcd(pdgm).and..not.isqcd(pdgid1).and..not.isqcd(pdgid2)) then
788 ! this is to avoid to do weird stuff for w+ w- z (or h h h)
789 ! this fix an issue for qq_zttxqq G1594.08
790 continue
791 elseif (ipart(2,imocl(n)).eq.0) then ! q a > q case
792 iqjets(ipart(1,imocl(n)))=0
751 else ! octet. want to be sure that both are tagged as jet before removing one793 else ! octet. want to be sure that both are tagged as jet before removing one
752 ! this prevent that both are removed in case of g > g h , g > q1 q2, q1 > a q1.794 ! this prevent that both are removed in case of g > g h , g > q1 q2, q1 > a q1.
753 ! at least one of the two should be kept as jet795 ! at least one of the two should be kept as jet
@@ -791,7 +833,57 @@
791 goodjet(imocl(n))=833 goodjet(imocl(n))=
792 $ (isjet(ipdgcl(imocl(n),igraphs(1),iproc)).and.834 $ (isjet(ipdgcl(imocl(n),igraphs(1),iproc)).and.
793 $ goodjet(idacl(n,1)).and.goodjet(idacl(n,2)))835 $ goodjet(idacl(n,1)).and.goodjet(idacl(n,2)))
794 endif836
837c check case with g > g g
838c where the hardest gluon is not goodjet but the other is.
839c in that case change ipart(1,) of the mother gluon
840c pure QCD jet
841c need to take care of the following case:
842c tttttttttt
843c gggggggggggggggg
844c ggggggggg tttttttttt
845c gggggggg
846c
847c in that case the up gluon can be tag as the hardest one
848c but this one is also lead to no QCD one.
849c so in that case we have to change ipart(1) to the sofest gluon
850 pdgm = ipdgcl(imocl(n),igraphs(1),iproc)
851 pdgid1 = ipdgcl(idacl(n,1),igraphs(1),iproc)
852 pdgid2 = ipdgcl(idacl(n,2),igraphs(1),iproc)
853 if (is_octet(pdgm).and.is_octet(pdgid1).and.is_octet(pdgid2))then
854c write(*,*) 'pure QCD vertex (2)'
855c write(*,*) pdgm , '>', pdgid1,' ', pdgid2
856c write(*,*) 'ipart', ipart(1,imocl(n)), ipart(2,imocl(n))
857c write(*,*) 'id', imocl(n), idacl(n,1),idacl(n,2)
858c write(*,*) 'ipart',ipart(1,imocl(n)),'/', ipart(2,imocl(n)), ipart(1,idacl(n,1)),'/', ipart(2,idacl(n,1)),
859c $ ipart(1,idacl(n,2)),'/', ipart(2,idacl(n,2))
860c write(*,*) 'googjet', goodjet(idacl(n,1)),goodjet(idacl(n,2))
861 if (ipart(1,imocl(n)).eq.ipart(1, idacl(n,1))) then
862 if (.not.goodjet(idacl(n,1)).and.goodjet(idacl(n,2))) then
863c write(*,*) 'ggg with hard jet set a QED the second jet lead to', ipart(1,idacl(n,2)), ipart(2,idacl(n,2))
864 do m =n_max_cl,n,-1
865 if(ipart(1,m).eq.ipart(1,imocl(n)).and.ipart(2,m).eq.ipart(2,imocl(n)))then
866 ipart(1,m) = ipart(1,idacl(n,2))
867 ipart(2,m) = ipart(2,idacl(n,2))
868 endif
869 enddo
870 endif
871 else
872 if (.not.goodjet(idacl(n,2)).and.goodjet(idacl(n,1))) then
873c write(*,*) 'ggg with hard jet set a QED the second jet lead to', ipart(1,idacl(n,1)), ipart(2,idacl(n,1))
874 do m =n_max_cl,n,-1
875 if(ipart(1,m).eq.ipart(1,imocl(n)).and.ipart(2,m).eq.ipart(2,imocl(n)))then
876 ipart(1,m) = ipart(1,idacl(n,1))
877 ipart(2,m) = ipart(2,idacl(n,1))
878 endif
879 enddo
880 endif
881 endif
882 endif
883
884
885
886 endif
795 enddo887 enddo
796888
797 if (btest(mlevel,4))then889 if (btest(mlevel,4))then
@@ -836,8 +928,8 @@
836 endif928 endif
837 enddo929 enddo
838 njetstore(iconfig)=njets930 njetstore(iconfig)=njets
839 if (btest(mlevel,4))931 if (btest(mlevel,4))
840 $ write(*,*) 'Storing jets: ',(iqjetstore(i,iconfig),i=1,njets)932 $ write(*,*) 'Storing jets: ',(iqjetstore(i,iconfig),i=1,njets)
841c Recluster without requiring chcluster933c Recluster without requiring chcluster
842 goto 100934 goto 100
843 else935 else
844936
=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f 2017-06-22 08:03:29 +0000
+++ Template/LO/SubProcesses/unwgt.f 2018-12-20 08:18:55 +0000
@@ -513,6 +513,10 @@
513513
514 double precision stot,m1,m2514 double precision stot,m1,m2
515 common/to_stot/stot,m1,m2515 common/to_stot/stot,m1,m2
516
517 INTEGER IMIRROR
518 INTEGER IPROC
519 COMMON/TO_MIRROR/IMIRROR, IPROC
516c520c
517c Data521c Data
518c522c
@@ -641,6 +645,13 @@
641 pb(4,isym(j,jsym))=pmass(j)645 pb(4,isym(j,jsym))=pmass(j)
642 enddo646 enddo
643 endif647 endif
648
649 if (IMIRROR.eq.2.and.pmass(1).ne.pmass(2)) then
650c Note that in this context isym(1,jsym) should never be "2" since the mass differ
651 pb(4,isym(1,jsym))=pmass(2)
652 pb(4,isym(2,jsym))=pmass(1)
653 endif
654
644c655c
645c Add info on resonant mothers656c Add info on resonant mothers
646c657c
647658
=== added file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f'
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f 1970-01-01 00:00:00 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f 2018-12-20 08:18:55 +0000
@@ -0,0 +1,880 @@
1cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2 subroutine analysis_begin(nwgt,weights_info)
3cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4 implicit none
5 integer nwgt
6 character*(*) weights_info(*)
7 integer i,l
8 character*6 cc(2)
9 data cc/'|T@NLO','|T@LO '/
10
11c Also specific perturbative orders can be directly plotted, adding for examples
12c the following further entries in the variable data
13c $ ' |T@QCD4QED0 ',
14c $ ,' |T@QCD2QED2 ',' |T@QCD0QED4 ',' |T@QCD6QED0 '
15c $ ,' |T@QCD4QED2 ',' |T@QCD2QED4 ',' |T@QCD0QED6 '
16c $ ,' |T@QCD8QED0 ',' |T@QCD6QED2 ',' |T@QCD4QED4 '
17c $ ,' |T@QCD2QED6 ',' |T@QCD0QED8 '
18c
19c See also line 376 in this file
20 call HwU_inithist(nwgt,weights_info)
21 do i=1,2
22 l=(i-1)*55
23c transverse momenta
24 call HwU_book(l+ 1,'total rate '//cc(i),1,0.5d0,1.5d0)
25 call HwU_book(l+ 2,'1st charged lepton log pT '//cc(i),25,-0.2d0,3.8d0)
26 call HwU_book(l+ 3,'2nd charged lepton log pT '//cc(i),25,-0.2d0,3.8d0)
27 call HwU_book(l+ 4,'3rd charged lepton log pT '//cc(i),25,-0.2d0,3.8d0)
28 call HwU_book(l+ 5,'4th charged lepton log pT '//cc(i),25,-0.2d0,3.8d0)
29 call HwU_book(l+ 6,'electron log pT '//cc(i),25,-0.2d0,3.8d0)
30 call HwU_book(l+ 7,'positron log pT '//cc(i),25,-0.2d0,3.8d0)
31 call HwU_book(l+ 8,'mu-plus log pT '//cc(i),25,-0.2d0,3.8d0)
32 call HwU_book(l+ 9,'mu-minus log pT '//cc(i),25,-0.2d0,3.8d0)
33 call HwU_book(l+14,'epem log pT '//cc(i),25,-0.2d0,3.8d0)
34 call HwU_book(l+15,'mupmum log pT '//cc(i),25,-0.2d0,3.8d0)
35 call HwU_book(l+16,'1st OSSF lep pair log pT '//cc(i),25,-0.2d0,3.8d0)
36 call HwU_book(l+17,'2nd OSSF lep pair log pT '//cc(i),25,-0.2d0,3.8d0)
37 call HwU_book(l+18,'epve log pT '//cc(i),25,-0.2d0,3.8d0)
38 call HwU_book(l+19,'mumvm log pT '//cc(i),25,-0.2d0,3.8d0)
39 call HwU_book(l+20,'1st SF lep-neu pair log pT '//cc(i),25,-0.2d0,3.8d0)
40 call HwU_book(l+21,'2nd SF lep-neu pair log pT '//cc(i),25,-0.2d0,3.8d0)
41 call HwU_book(l+22,'epemmupmum log pT '//cc(i),25,-0.2d0,3.8d0)
42 call HwU_book(l+23,'epvemumvm log pT '//cc(i),25,-0.2d0,3.8d0)
43 call HwU_book(l+24,'top quark log pT '//cc(i),25,-0.2d0,3.8d0)
44 call HwU_book(l+25,'anti-top quark log pT '//cc(i),25,-0.2d0,3.8d0)
45 call HwU_book(l+26,'ttbar log pT '//cc(i),25,-0.2d0,3.8d0)
46 call HwU_book(l+27,'log missing Et (neutrinos) '//cc(i),25,-0.2d0,3.8d0)
47 call HwU_book(l+28,'1st higgs log pT '//cc(i),25,-0.2d0,3.8d0)
48 call HwU_book(l+29,'2nd higgs log pT '//cc(i),25,-0.2d0,3.8d0)
49 call HwU_book(l+30,'higgs-pair log pT '//cc(i),25,-0.2d0,3.8d0)
50 call HwU_book(l+31,'1st V-boson log pT '//cc(i),25,-0.2d0,3.8d0)
51 call HwU_book(l+32,'2nd V-boson log pT '//cc(i),25,-0.2d0,3.8d0)
52 call HwU_book(l+33,'3rd V-boson log pT '//cc(i),25,-0.2d0,3.8d0)
53 call HwU_book(l+34,'1st jet log pT '//cc(i),25,-0.2d0,3.8d0)
54 call HwU_book(l+35,'2nd jet log pT '//cc(i),25,-0.2d0,3.8d0)
55 call HwU_book(l+36,'3rd jet log pT '//cc(i),25,-0.2d0,3.8d0)
56c
57c invariant masses
58 call HwU_book(l+37,'epem log M '//cc(i),25,-0.2d0,3.8d0)
59 call HwU_book(l+38,'mupmum log M '//cc(i),25,-0.2d0,3.8d0)
60 call HwU_book(l+39,'1st OSSF lep pair log M '//cc(i),25,-0.2d0,3.8d0)
61 call HwU_book(l+40,'2nd OSSF lep pair log M '//cc(i),25,-0.2d0,3.8d0)
62 call HwU_book(l+41,'epve log M '//cc(i),25,-0.2d0,3.8d0)
63 call HwU_book(l+42,'mumvm log M '//cc(i),25,-0.2d0,3.8d0)
64 call HwU_book(l+43,'1st SF lep-neu pair log M '//cc(i),25,-0.2d0,3.8d0)
65 call HwU_book(l+44,'2nd SF lep-neu pair log M '//cc(i),25,-0.2d0,3.8d0)
66 call HwU_book(l+45,'epemmupmum log M '//cc(i),25,-0.2d0,3.8d0)
67 call HwU_book(l+46,'epvemumvm log M '//cc(i),25,-0.2d0,3.8d0)
68 call HwU_book(l+47,'ttbar log M '//cc(i),25,-0.2d0,3.8d0)
69 call HwU_book(l+48,'higgs-pair log M '//cc(i),25,-0.2d0,3.8d0)
70 call HwU_book(l+49,'j1j2 log M '//cc(i),25,-0.2d0,3.8d0)
71 call HwU_book(l+50,'j2j3 log M '//cc(i),25,-0.2d0,3.8d0)
72 call HwU_book(l+51,'j1j3 log M '//cc(i),25,-0.2d0,3.8d0)
73 call HwU_book(l+52,'j1j2j3 log M '//cc(i),25,-0.2d0,3.8d0)
74 call HwU_book(l+53,'3w log M '//cc(i),25,-0.2d0,3.8d0)
75c
76c HT
77 call HwU_book(l+54,'log HT (partons) '//cc(i),25,-0.2d0,3.8d0)
78 call HwU_book(l+55,'log HT (reconstructed) '//cc(i),25,-0.2d0,3.8d0)
79 enddo
80 return
81 end
82
83
84
85
86cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
87 subroutine analysis_end(dummy)
88cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
89 implicit none
90 double precision dummy
91 call HwU_write_file
92 return
93 end
94
95
96cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
97 subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
98cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
99 implicit none
100 include 'nexternal.inc'
101 include 'cuts.inc'
102 integer istatus(nexternal)
103 integer iPDG(nexternal)
104 double precision p(0:4,nexternal)
105 double precision wgts(*)
106 integer ibody
107 integer i,j,k,kk,l
108 double precision www,pQCD(0:3,nexternal),palg,rfj,sycut,yjmax
109 $ ,pjet(0:3,nexternal),tmp,ptlep(4),ptem,ptep,ptmm,ptmp,ptepem
110 $ ,ptmpmm,ptepve,ptmmvm,ptepemmpmm,ptepvemmvm,ptt,ptat,pttt
111 $ ,etmiss,pth(2),pthh,ptv(3),ptjet(3),Mepem,Mmpmm,Mepve,Mmmvm
112 $ ,Mepemmpmm,Mepvemmvm,Mtt,Mhh,Mj1j2,Mj1j3,Mj2j3,Mj1j2j3,Mvvv
113 $ ,HTparton,HTreco,p_reco(0:4,nexternal)
114 integer nQCD,jet(nexternal),njet,itop,iatop,iem,iep,imp,imm,ive
115 $ ,ivm,iv1,iv2,iv3,ih1,ih2,il,ipdg_reco(nexternal)
116 double precision getptv4,getptv4_2,getptv4_4,getinvm4_2
117 $ ,getinvm4_3,getinvm4_4,l10
118 external getptv4,getptv4_2,getptv4_4,getinvm4_2,getinvm4_3
119 $ ,getinvm4_4,l10
120 integer orders_tag_plot
121 common /corderstagplot/ orders_tag_plot
122c First, try to recombine photons with leptons
123 if (.not.quarkphreco) then
124 write (*,*) 'quark-photon recombination is turned off. '/
125 $ /'Do need it'
126 stop
127 endif
128 if (.not. lepphreco) then
129 write (*,*) 'lepton-photon recombination is turned off. '
130 $ //'Do need it.'
131 stop
132 endif
133 call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
134 $ p, iPDG, p_reco, iPDG_reco)
135
136c Put all (light) QCD partons(+photon) in momentum array for jet clustering.
137 nQCD=0
138 do j=nincoming+1,nexternal
139 if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 .or.
140 $ ipdg_reco(j).eq.22)then
141 nQCD=nQCD+1
142 do i=0,3
143 pQCD(i,nQCD)=p_reco(i,j)
144 enddo
145 endif
146 enddo
147
148C---CLUSTER THE EVENT
149 palg = jetalgo
150 rfj = jetradius
151 sycut = ptj
152 yjmax = etaj
153
154c******************************************************************************
155c call FASTJET to get all the jets
156c
157c INPUT:
158c input momenta: pQCD(0:3,nexternal), energy is 0th component
159c number of input momenta: nQCD
160c radius parameter: rfj
161c minumum jet pt: sycut
162c jet algorithm: palg, 1.0=kt, 0.0=C/A, -1.0 = anti-kt
163c
164c OUTPUT:
165c jet momenta: pjet(0:3,nexternal), E is 0th cmpnt
166c the number of jets (with pt > SYCUT): njet
167c the jet for a given particle 'i': jet(i), note that this is
168c the particle in pQCD, which doesn't necessarily correspond to the particle
169c label in the process
170c
171 call amcatnlo_fastjetppgenkt_etamax(pQCD,nQCD,rfj,sycut,yjmax,palg
172 $ ,pjet,njet,jet)
173c
174c******************************************************************************
175
176c Look for the other physics objects
177 itop=0
178 iatop=0
179 iem=0
180 iep=0
181 iem=0
182 iep=0
183 imp=0
184 imm=0
185 ive=0
186 ivm=0
187 iv1=0
188 iv2=0
189 iv3=0
190 ih1=0
191 ih2=0
192 do i=1,nexternal
193 if (ipdg_reco(i).eq.6) then
194 itop=i
195 elseif(ipdg_reco(i).eq.-6) then
196 iatop=i
197 elseif(ipdg_reco(i).eq.11) then
198 iem=i
199 elseif(ipdg_reco(i).eq.13) then
200 imm=i
201 elseif(ipdg_reco(i).eq.-11) then
202 iep=i
203 elseif(ipdg_reco(i).eq.-13) then
204 imp=i
205 elseif(abs(ipdg_reco(i)).eq.12) then
206 ive=i
207 elseif(abs(ipdg_reco(i)).eq.14) then
208 ivm=i
209 elseif(abs(ipdg_reco(i)).eq.24 .or. ipdg_reco(i).eq.23) then
210 if (iv1.eq.0) then
211 iv1=i
212 else
213 if (iv2.eq.0) then
214 iv2=i
215 else
216 if(iv3.eq.0) then
217 iv3=i
218 else
219 write (*,*) 'too many vector bosons'
220 stop
221 endif
222 endif
223 endif
224 elseif(abs(ipdg_reco(i)).eq.25) then
225 if (ih1.eq.0) then
226 ih1=i
227 else
228 if (ih2.eq.0) then
229 ih2=i
230 else
231 write (*,*) 'too many higgs bosons'
232 stop
233 endif
234 endif
235 endif
236 enddo
237 if (itop.ne.0) ptt=getptv4(p_reco(0,itop))
238 if (iatop.ne.0) ptat=getptv4(p_reco(0,iatop))
239 if (itop.ne.0 .and. iatop.ne.0) then
240 pttt=getptv4_2(p_reco(0,itop),p_reco(0,iatop))
241 Mtt=getinvm4_2(p_reco(0,itop),p_reco(0,iatop))
242 endif
243 if (iem.ne.0) ptem=getptv4(p_reco(0,iem))
244 if (iep.ne.0) ptep=getptv4(p_reco(0,iep))
245 if (imm.ne.0) ptmm=getptv4(p_reco(0,imm))
246 if (imp.ne.0) ptmp=getptv4(p_reco(0,imp))
247 if (iem.ne.0 .and. iep.ne.0) then
248 ptepem=getptv4_2(p_reco(0,iem),p_reco(0,iep))
249 Mepem=getinvm4_2(p_reco(0,iem),p_reco(0,iep))
250 endif
251 if (imm.ne.0 .and. imp.ne.0) then
252 ptmpmm=getptv4_2(p_reco(0,imm),p_reco(0,imp))
253 Mmpmm=getinvm4_2(p_reco(0,imm),p_reco(0,imp))
254 endif
255 if (iep.ne.0 .and. ive.ne.0) then
256 ptepve=getptv4_2(p_reco(0,iep),p_reco(0,ive))
257 Mepve=getinvm4_2(p_reco(0,iep),p_reco(0,ive))
258 endif
259 if (imm.ne.0 .and. ivm.ne.0) then
260 ptmmvm=getptv4_2(p_reco(0,imm),p_reco(0,ivm))
261 Mmmvm=getinvm4_2(p_reco(0,imm),p_reco(0,ivm))
262 endif
263 if (iem.ne.0 .and. iep.ne.0 .and. imm.ne.0 .and. imp.ne.0) then
264 ptepemmpmm=getptv4_4(p_reco(0,iem),p_reco(0,iep),p_reco(0,imm),p_reco(0,imp))
265 Mepemmpmm=getinvm4_4(p_reco(0,iem),p_reco(0,iep),p_reco(0,imm),p_reco(0,imp))
266 endif
267 if (ive.ne.0 .and. iep.ne.0 .and. imm.ne.0 .and. ivm.ne.0) then
268 ptepvemmvm=getptv4_4(p_reco(0,iep),p_reco(0,ive),p_reco(0,imm),p_reco(0,ivm))
269 Mepvemmvm=getinvm4_4(p_reco(0,iep),p_reco(0,ive),p_reco(0,imm),p_reco(0,ivm))
270 endif
271
272 il=0
273 if (iem.ne.0) then
274 il=il+1
275 ptlep(il)=ptem
276 endif
277 if (iep.ne.0) then
278 il=il+1
279 ptlep(il)=ptep
280 endif
281 if (imm.ne.0) then
282 il=il+1
283 ptlep(il)=ptmm
284 endif
285 if (imp.ne.0) then
286 il=il+1
287 ptlep(il)=ptmp
288 endif
289 if (il.gt.1) then
290 do i=1,il-1
291 do j=1,il-i
292 if (ptlep(j).lt.ptlep(j+1)) then
293 tmp=ptlep(j)
294 ptlep(j)=ptlep(j+1)
295 ptlep(j+1)=tmp
296 endif
297 enddo
298 enddo
299 endif
300
301
302c missing Et
303 if (ive.ne.0 .and. ivm.ne.0) then
304 etmiss=getptv4_2(p_reco(0,ivm),p_reco(0,ive))
305 elseif (ive.ne.0 .or. ivm.ne.0) then
306 etmiss=getptv4(p_reco(0,ive+ivm))
307 endif
308
309 if (ih1.ne.0)pth(1)=getptv4(p_reco(0,ih1))
310 if (ih2.ne.0)pth(2)=getptv4(p_reco(0,ih2))
311 if (ih1.ne.0 .and. ih2.ne.0) then
312 pthh=getptv4_2(p_reco(0,ih1),p_reco(0,ih2))
313 Mhh=getinvm4_2(p_reco(0,ih1),p_reco(0,ih2))
314c order the higgs bosons (if there are 2)
315 if (pth(1).lt.pth(2)) then
316 tmp=pth(1)
317 pth(1)=pth(2)
318 pth(2)=tmp
319 endif
320 endif
321 if (iv1.ne.0) ptv(1)=getptv4(p_reco(0,iv1))
322 if (iv2.ne.0) ptv(2)=getptv4(p_reco(0,iv2))
323 if (iv3.ne.0) ptv(3)=getptv4(p_reco(0,iv3))
324 if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) then
325 Mvvv=getinvm4_3(p_reco(0,iv1),p_reco(0,iv2),p_reco(0,iv3))
326c order the vector bosons (if there are 3)
327 do i=1,2
328 do j=1,3-i
329 if (ptv(j).lt.ptv(j+1)) then
330 tmp=ptv(j)
331 ptv(j)=ptv(j+1)
332 ptv(j+1)=tmp
333 endif
334 enddo
335 enddo
336 elseif (iv1.ne.0 .and. iv2.ne.0) then
337c order the vector bosons (if there are 2)
338 if (ptv(1).lt.ptv(2)) then
339 tmp=ptv(1)
340 ptv(1)=ptv(2)
341 ptv(2)=tmp
342 endif
343 endif
344 do i=1,njet
345 ptjet(i)=getptv4(pjet(0,i))
346 enddo
347 if (njet.ge.2) then
348 Mj1j2=getinvm4_2(pjet(0,1),pjet(0,2))
349 endif
350 if (njet.ge.3) then
351 Mj2j3=getinvm4_2(pjet(0,2),pjet(0,3))
352 Mj1j3=getinvm4_2(pjet(0,1),pjet(0,3))
353 Mj1j2j3=getinvm4_3(pjet(0,1),pjet(0,2),pjet(0,3))
354 endif
355c
356 HTparton=0d0
357 HTreco=0d0
358 do i=1,nexternal
359 HTparton=HTparton+getptv4(p(0,i))
360 if (abs(ipdg_reco(i)).gt.5 .and. ipdg_reco(i).ne.21 .and.
361 $ ipdg_reco(i).ne.22 .and. abs(ipdg_reco(i)).ne.12 .and.
362 $ abs(ipdg_reco(i)).ne.14 .and. abs(ipdg_reco(i)).ne.16)
363 $ then
364 HTreco=HTreco+getptv4(p_reco(0,i))
365 endif
366 enddo
367 do i=1,njet
368 HTreco=HTreco+getptv4(pjet(0,i))
369 enddo
370 if (ive.ne.0 .or. ivm.ne.0) HTreco=HTreco+etmiss
371
372 do i=1,2
373 l=(i-1)*55
374 if (ibody.ne.3 .and.i.eq.2) cycle
375
376c How to tag orders (QCD+QED*100)
377c
378c if (i.eq. 3.and.orders_tag_plot.ne.4) cycle
379c if (i.eq. 4.and.orders_tag_plot.ne.202) cycle
380c if (i.eq. 5.and.orders_tag_plot.ne.400) cycle
381c if (i.eq. 6.and.orders_tag_plot.ne.6) cycle
382c if (i.eq. 7.and.orders_tag_plot.ne.204) cycle
383c if (i.eq. 8.and.orders_tag_plot.ne.402) cycle
384c if (i.eq. 9.and.orders_tag_plot.ne.600) cycle
385c if (i.eq.10.and.orders_tag_plot.ne.8) cycle
386c if (i.eq.11.and.orders_tag_plot.ne.206) cycle
387c if (i.eq.12.and.orders_tag_plot.ne.404) cycle
388c if (i.eq.13.and.orders_tag_plot.ne.602) cycle
389c if (i.eq.14.and.orders_tag_plot.ne.800) cycle
390
391
392
393c total rate
394 call HwU_fill(l+ 1,1d0,wgts)
395c transverse momenta
396 if (il .ge.1) call HwU_fill(l+ 2,l10(ptlep(1)),wgts)
397 if (il .ge.2) call HwU_fill(l+ 3,l10(ptlep(2)),wgts)
398 if (il .ge.3) call HwU_fill(l+ 4,l10(ptlep(3)),wgts)
399 if (il .ge.4) call HwU_fill(l+ 5,l10(ptlep(4)),wgts)
400 if (iem.ne.0) call HwU_fill(l+ 6,l10(ptem),wgts)
401 if (iep.ne.0) call HwU_fill(l+ 7,l10(ptep),wgts)
402 if (imp.ne.0) call HwU_fill(l+ 8,l10(ptmp),wgts)
403 if (imm.ne.0) call HwU_fill(l+ 9,l10(ptmm),wgts)
404 if (iep.ne.0 .and. iem.ne.0) call HwU_fill(l+14,l10(ptepem),wgts)
405 if (imp.ne.0 .and. imm.ne.0) call HwU_fill(l+15,l10(ptmpmm),wgts)
406 if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0) then
407 if (ptepem.gt.ptmpmm) then
408 call HwU_fill(l+16,l10(ptepem),wgts)
409 call HwU_fill(l+17,l10(ptmpmm),wgts)
410 else
411 call HwU_fill(l+17,l10(ptepem),wgts)
412 call HwU_fill(l+16,l10(ptmpmm),wgts)
413 endif
414 elseif (iep.ne.0 .and. iem.ne.0) then
415 call HwU_fill(l+17,l10(ptepem),wgts)
416 elseif (imp.ne.0 .and. imm.ne.0) then
417 call HwU_fill(l+17,l10(ptmpmm),wgts)
418 endif
419 if (iep.ne.0 .and. ive.ne.0) call HwU_fill(l+18,l10(ptepve),wgts)
420 if (imm.ne.0 .and. ivm.ne.0) call HwU_fill(l+19,l10(ptmmvm),wgts)
421 if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0) then
422 if (ptepve.gt.ptmmvm) then
423 call HwU_fill(l+20,l10(ptepve),wgts)
424 call HwU_fill(l+21,l10(ptmmvm),wgts)
425 else
426 call HwU_fill(l+21,l10(ptepve),wgts)
427 call HwU_fill(l+20,l10(ptmmvm),wgts)
428 endif
429 elseif (iep.ne.0 .and. ive.ne.0) then
430 call HwU_fill(l+20,l10(ptepve),wgts)
431 elseif (imm.ne.0 .and. ivm.ne.0) then
432 call HwU_fill(l+20,l10(ptmmvm),wgts)
433 endif
434 if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0)
435 $ call HwU_fill(l+22,l10(ptepemmpmm),wgts)
436 if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0)
437 $ call HwU_fill(l+23,l10(ptepvemmvm),wgts)
438 if (itop.ne.0) call HwU_fill(l+24,l10(ptt),wgts)
439 if (iatop.ne.0) call HwU_fill(l+25,l10(ptat),wgts)
440 if (itop.ne.0 .and. iatop.ne.0) call HwU_fill(l+26,l10(pttt),wgts)
441 if (ive.ne.0 .or. ivm .ne.0) call HwU_fill(l+27,l10(etmiss),wgts)
442 if (ih1.ne.0) call HwU_fill(l+28,l10(pth(1)),wgts)
443 if (ih2.ne.0) call HwU_fill(l+29,l10(pth(2)),wgts)
444 if (ih1.ne.0 .and. ih2.ne.0) call HwU_fill(l+30,l10(pthh),wgts)
445 if (iv1.ne.0) call HwU_fill(l+31,l10(ptv(1)),wgts)
446 if (iv2.ne.0) call HwU_fill(l+32,l10(ptv(2)),wgts)
447 if (iv3.ne.0) call HwU_fill(l+33,l10(ptv(3)),wgts)
448 if (njet.ge.1) call HwU_fill(l+34,l10(ptjet(1)),wgts)
449 if (njet.ge.2) call HwU_fill(l+35,l10(ptjet(2)),wgts)
450 if (njet.ge.3) call HwU_fill(l+36,l10(ptjet(3)),wgts)
451c invariant masses
452 if (iep.ne.0 .and. iem.ne.0) call HwU_fill(l+37,l10(Mepem),wgts)
453 if (imp.ne.0 .and. imm.ne.0) call HwU_fill(l+38,l10(Mmpmm),wgts)
454 if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0) then
455 if (Mepem.gt.Mmpmm) then
456 call HwU_fill(l+39,l10(Mepem),wgts)
457 call HwU_fill(l+40,l10(Mmpmm),wgts)
458 else
459 call HwU_fill(l+40,l10(Mepem),wgts)
460 call HwU_fill(l+39,l10(Mmpmm),wgts)
461 endif
462 elseif (iep.ne.0 .and. iem.ne.0) then
463 call HwU_fill(l+39,l10(Mepem),wgts)
464 elseif (imp.ne.0 .and. imm.ne.0) then
465 call HwU_fill(l+39,l10(Mmpmm),wgts)
466 endif
467 if (iep.ne.0 .and. ive.ne.0) call HwU_fill(l+41,l10(Mepve),wgts)
468 if (imm.ne.0 .and. ivm.ne.0) call HwU_fill(l+42,l10(Mmmvm),wgts)
469 if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0) then
470 if (Mepve.gt.Mmmvm) then
471 call HwU_fill(l+43,l10(Mepve),wgts)
472 call HwU_fill(l+44,l10(Mmmvm),wgts)
473 else
474 call HwU_fill(l+44,l10(Mepve),wgts)
475 call HwU_fill(l+43,l10(Mmmvm),wgts)
476 endif
477 elseif (iep.ne.0 .and. ive.ne.0) then
478 call HwU_fill(l+43,l10(Mepve),wgts)
479 elseif (imm.ne.0 .and. ivm.ne.0) then
480 call HwU_fill(l+43,l10(Mmmvm),wgts)
481 endif
482 if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0)
483 $ call HwU_fill(l+45,l10(Mepemmpmm),wgts)
484 if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0)
485 $ call HwU_fill(l+46,l10(Mepvemmvm),wgts)
486 if (itop.ne.0 .and. iatop.ne.0) call HwU_fill(l+47,l10(Mtt),wgts)
487 if (ih1.ne.0 .and. ih2.ne.0) call HwU_fill(l+48,l10(Mhh),wgts)
488 if (njet.ge.2) call HwU_fill(l+49,l10(Mj1j2),wgts)
489 if (njet.ge.3) call HwU_fill(l+50,l10(Mj2j3),wgts)
490 if (njet.ge.3) call HwU_fill(l+51,l10(Mj1j3),wgts)
491 if (njet.ge.3) call HwU_fill(l+52,l10(Mj1j2j3),wgts)
492 if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) call HwU_fill(l+53,l10(Mvvv),wgts)
493c HT
494 call HwU_fill(l+54,l10(HTparton),wgts)
495 call HwU_fill(l+55,l10(HTreco),wgts)
496 enddo
497
498 999 return
499 end
500
501
502 double precision function l10(var)
503 implicit none
504 double precision var
505 if (var.gt.0) then
506 l10=log10(var)
507 else
508 l10=-99d99
509 endif
510 return
511 end
512
513 function getrapidity(en,pl)
514 implicit none
515 real*8 getrapidity,en,pl,tiny,xplus,xminus,y
516 parameter (tiny=1.d-8)
517c
518 xplus=en+pl
519 xminus=en-pl
520 if(xplus.gt.tiny.and.xminus.gt.tiny)then
521 if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny )then
522 y=0.5d0*log( xplus/xminus )
523 else
524 y=sign(1.d0,pl)*1.d8
525 endif
526 else
527 y=sign(1.d0,pl)*1.d8
528 endif
529 getrapidity=y
530 return
531 end
532
533
534 function getpseudorap(en,ptx,pty,pl)
535 implicit none
536 real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
537 parameter (tiny=1.d-5)
538c
539 pt=sqrt(ptx**2+pty**2)
540 if(pt.lt.tiny.and.abs(pl).lt.tiny)then
541 eta=sign(1.d0,pl)*1.d8
542 else
543 th=atan2(pt,pl)
544 eta=-log(tan(th/2.d0))
545 endif
546 getpseudorap=eta
547 return
548 end
549
550
551 function getinvm(en,ptx,pty,pl)
552 implicit none
553 real*8 getinvm,en,ptx,pty,pl,tiny,tmp
554 parameter (tiny=1.d-5)
555c
556 tmp=en**2-ptx**2-pty**2-pl**2
557 if(tmp.gt.0.d0)then
558 tmp=sqrt(tmp)
559 elseif(tmp.gt.-tiny)then
560 tmp=0.d0
561 else
562 write(*,*)'Attempt to compute a negative mass'
563 stop
564 endif
565 getinvm=tmp
566 return
567 end
568
569
570 function getdelphi(ptx1,pty1,ptx2,pty2)
571 implicit none
572 real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
573 parameter (tiny=1.d-5)
574c
575 pt1=sqrt(ptx1**2+pty1**2)
576 pt2=sqrt(ptx2**2+pty2**2)
577 if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
578 tmp=ptx1*ptx2+pty1*pty2
579 tmp=tmp/(pt1*pt2)
580 if(abs(tmp).gt.1.d0+tiny)then
581 write(*,*)'Cosine larger than 1'
582 stop
583 elseif(abs(tmp).ge.1.d0)then
584 tmp=sign(1.d0,tmp)
585 endif
586 tmp=acos(tmp)
587 else
588 tmp=1.d8
589 endif
590 getdelphi=tmp
591 return
592 end
593
594
595 function getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
596 implicit none
597 real*8 getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
598 # getpseudorap,getdelphi
599c
600 deta=getpseudorap(en1,ptx1,pty1,pl1)-
601 # getpseudorap(en2,ptx2,pty2,pl2)
602 dphi=getdelphi(ptx1,pty1,ptx2,pty2)
603 getdr=sqrt(dphi**2+deta**2)
604 return
605 end
606
607
608 function getdry(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
609 implicit none
610 real*8 getdry,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
611 # getrapidity,getdelphi
612c
613 deta=getrapidity(en1,pl1)-
614 # getrapidity(en2,pl2)
615 dphi=getdelphi(ptx1,pty1,ptx2,pty2)
616 getdry=sqrt(dphi**2+deta**2)
617 return
618 end
619
620
621 function getptv(p)
622 implicit none
623 real*8 getptv,p(5)
624c
625 getptv=sqrt(p(1)**2+p(2)**2)
626 return
627 end
628
629
630 function getpseudorapv(p)
631 implicit none
632 real*8 getpseudorapv,p(5)
633 real*8 getpseudorap
634c
635 getpseudorapv=getpseudorap(p(4),p(1),p(2),p(3))
636 return
637 end
638
639
640 function getrapidityv(p)
641 implicit none
642 real*8 getrapidityv,p(5)
643 real*8 getrapidity
644c
645 getrapidityv=getrapidity(p(4),p(3))
646 return
647 end
648
649
650 function getdrv(p1,p2)
651 implicit none
652 real*8 getdrv,p1(5),p2(5)
653 real*8 getdr
654c
655 getdrv=getdr(p1(4),p1(1),p1(2),p1(3),
656 # p2(4),p2(1),p2(2),p2(3))
657 return
658 end
659
660
661 function getinvmv(p)
662 implicit none
663 real*8 getinvmv,p(5)
664 real*8 getinvm
665c
666 getinvmv=getinvm(p(4),p(1),p(2),p(3))
667 return
668 end
669
670
671 function getdelphiv(p1,p2)
672 implicit none
673 real*8 getdelphiv,p1(5),p2(5)
674 real*8 getdelphi
675c
676 getdelphiv=getdelphi(p1(1),p1(2),
677 # p2(1),p2(2))
678 return
679 end
680
681
682 function getptv4(p)
683 implicit none
684 real*8 getptv4,p(0:3)
685c
686 getptv4=sqrt(max(p(1)**2+p(2)**2,0d0))
687 return
688 end
689
690 function getptv4_2(p1,p2)
691 implicit none
692 real*8 getptv4_2,p1(0:3),p2(0:3),psum(0:3)
693 integer i
694 do i=0,3
695 psum(i)=p1(i)+p2(i)
696 enddo
697 getptv4_2=sqrt(max(psum(1)**2+psum(2)**2,0d0))
698 return
699 end
700
701 function getptv4_4(p1,p2,p3,p4)
702 implicit none
703 real*8 getptv4_4,p1(0:3),p2(0:3),p3(0:3),p4(0:3),psum(0:3)
704 integer i
705 do i=0,3
706 psum(i)=p1(i)+p2(i)+p3(i)+p4(i)
707 enddo
708 getptv4_4=sqrt(max(psum(1)**2+psum(2)**2,0d0))
709 return
710 end
711
712
713 function getpseudorapv4(p)
714 implicit none
715 real*8 getpseudorapv4,p(0:3)
716 real*8 getpseudorap
717c
718 getpseudorapv4=getpseudorap(p(0),p(1),p(2),p(3))
719 return
720 end
721
722
723 function getrapidityv4(p)
724 implicit none
725 real*8 getrapidityv4,p(0:3)
726 real*8 getrapidity
727c
728 getrapidityv4=getrapidity(p(0),p(3))
729 return
730 end
731
732
733 function getdrv4(p1,p2)
734 implicit none
735 real*8 getdrv4,p1(0:3),p2(0:3)
736 real*8 getdr
737c
738 getdrv4=getdr(p1(0),p1(1),p1(2),p1(3),
739 # p2(0),p2(1),p2(2),p2(3))
740 return
741 end
742
743
744 function getinvm4(p)
745 implicit none
746 real*8 getinvm4,p(0:3)
747 real*8 getinvm
748c
749 getinvm4=getinvm(p(0),p(1),p(2),p(3))
750 return
751 end
752
753 function getinvm4_2(p1,p2)
754 implicit none
755 real*8 getinvm4_2,p1(0:3),p2(0:3),p(0:3)
756 real*8 getinvm
757 integer i
758 do i=0,3
759 p(i)=p1(i)+p2(i)
760 enddo
761 getinvm4_2=getinvm(p(0),p(1),p(2),p(3))
762 return
763 end
764
765 function getinvm4_3(p1,p2,p3)
766 implicit none
767 real*8 getinvm4_3,p1(0:3),p2(0:3),p3(0:3),p(0:3)
768 real*8 getinvm
769 integer i
770 do i=0,3
771 p(i)=p1(i)+p2(i)+p3(i)
772 enddo
773 getinvm4_3=getinvm(p(0),p(1),p(2),p(3))
774 return
775 end
776
777 function getinvm4_4(p1,p2,p3,p4)
778 implicit none
779 real*8 getinvm4_4,p1(0:3),p2(0:3),p3(0:3),p4(0:3),p(0:3)
780 real*8 getinvm
781 integer i
782 do i=0,3
783 p(i)=p1(i)+p2(i)+p3(i)+p4(i)
784 enddo
785 getinvm4_4=getinvm(p(0),p(1),p(2),p(3))
786 return
787 end
788
789
790 function getdelphiv4(p1,p2)
791 implicit none
792 real*8 getdelphiv4,p1(0:3),p2(0:3)
793 real*8 getdelphi
794c
795 getdelphiv4=getdelphi(p1(1),p1(2),
796 # p2(1),p2(2))
797 return
798 end
799
800
801 function getcosv4(q1,q2)
802 implicit none
803 real*8 getcosv4,q1(0:3),q2(0:3)
804 real*8 xnorm1,xnorm2,tmp
805c
806 if(q1(0).lt.0.d0.or.q2(0).lt.0.d0)then
807 getcosv4=-1.d10
808 return
809 endif
810 xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
811 xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
812 if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
813 tmp=-1.d10
814 else
815 tmp=q1(1)*q2(1)+q1(2)*q2(2)+q1(3)*q2(3)
816 tmp=tmp/(xnorm1*xnorm2)
817 if(abs(tmp).gt.1.d0.and.abs(tmp).le.1.001d0)then
818 tmp=sign(1.d0,tmp)
819 elseif(abs(tmp).gt.1.001d0)then
820 write(*,*)'Error in getcosv4',tmp
821 stop
822 endif
823 endif
824 getcosv4=tmp
825 return
826 end
827
828
829
830 function getmod(p)
831 implicit none
832 double precision p(0:3),getmod
833
834 getmod=sqrt(p(1)**2+p(2)**2+p(3)**2)
835
836 return
837 end
838
839
840
841 subroutine getperpenv4(q1,q2,qperp)
842c Normal to the plane defined by \vec{q1},\vec{q2}
843 implicit none
844 real*8 q1(0:3),q2(0:3),qperp(0:3)
845 real*8 xnorm1,xnorm2
846 integer i
847c
848 xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
849 xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
850 if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
851 do i=1,4
852 qperp(i)=-1.d10
853 enddo
854 else
855 qperp(1)=q1(2)*q2(3)-q1(3)*q2(2)
856 qperp(2)=q1(3)*q2(1)-q1(1)*q2(3)
857 qperp(3)=q1(1)*q2(2)-q1(2)*q2(1)
858 do i=1,3
859 qperp(i)=qperp(i)/(xnorm1*xnorm2)
860 enddo
861 qperp(0)=1.d0
862 endif
863 return
864 end
865
866
867
868
869 subroutine getwedge(p1,p2,pout)
870 implicit none
871 real*8 p1(0:3),p2(0:3),pout(0:3)
872
873 pout(1)=p1(2)*p2(3)-p1(3)*p2(2)
874 pout(2)=p1(3)*p2(1)-p1(1)*p2(3)
875 pout(3)=p1(1)*p2(2)-p1(2)*p2(1)
876 pout(0)=0d0
877
878 return
879 end
880
0881
=== modified file 'Template/NLO/MCatNLO/Makefile_MadFKS'
--- Template/NLO/MCatNLO/Makefile_MadFKS 2016-02-24 00:19:50 +0000
+++ Template/NLO/MCatNLO/Makefile_MadFKS 2018-12-20 08:18:55 +0000
@@ -1,7 +1,7 @@
1-include ../../Source/make_opts1-include ../../Source/make_opts
22
3DEBUG=3DEBUG=
4FF=$(FC) $(FFLAGS) $(DEBUG)4FF=$(FC) $(FFLAGS) $(DEBUG) -std=legacy
5CPP=$(CXX) $(CXXFLAGS) $(DEBUG) $(INCLOPTION)5CPP=$(CXX) $(CXXFLAGS) $(DEBUG) $(INCLOPTION)
6CC=$(CXX) $(CFLAGS) $(DEBUG) $(INCLOPTION)6CC=$(CXX) $(CFLAGS) $(DEBUG) $(INCLOPTION)
77
88
=== modified file 'Template/NLO/SubProcesses/MC_integer.f'
--- Template/NLO/SubProcesses/MC_integer.f 2015-02-09 18:56:09 +0000
+++ Template/NLO/SubProcesses/MC_integer.f 2018-12-20 08:18:55 +0000
@@ -68,6 +68,7 @@
68 common/to_readgrid/flat_grid !Tells if grid read from file68 common/to_readgrid/flat_grid !Tells if grid read from file
69 if (this_dim.lt.1.or.this_dim.gt.maxdim) then69 if (this_dim.lt.1.or.this_dim.gt.maxdim) then
70 write (*,*) 'Increase maxdim in MC_integer.f',maxdim,this_dim70 write (*,*) 'Increase maxdim in MC_integer.f',maxdim,this_dim
71 stop 1
71 endif72 endif
72c Set the number of intervales for all the dimensions to zero the very73c Set the number of intervales for all the dimensions to zero the very
73c first time this subroutine is called74c first time this subroutine is called
@@ -283,12 +284,8 @@
283 enddo284 enddo
284c285c
285c Reset the accumulated results because we start new iteration.286c Reset the accumulated results because we start new iteration.
286 do this_dim=1,maxdim287 call empty_MC_integer
287 do i=0,nintervals(this_dim)288c
288 acc(i,this_dim)=0d0
289 ncall(i,this_dim)=0
290 enddo
291 enddo
292 return289 return
293 999 write (*,*) 'Cannot open "grid.MC_integer" file'290 999 write (*,*) 'Cannot open "grid.MC_integer" file'
294 stop291 stop
295292
=== modified file 'Template/NLO/SubProcesses/analysis_lhe.f'
--- Template/NLO/SubProcesses/analysis_lhe.f 2017-06-19 16:05:33 +0000
+++ Template/NLO/SubProcesses/analysis_lhe.f 2018-12-20 08:18:55 +0000
@@ -7,9 +7,10 @@
7 implicit none7 implicit none
8 integer nwgt8 integer nwgt
9 character*(*) weights_info(*)9 character*(*) weights_info(*)
10 logical passed_unwgt
10 integer nwgts,nevents11 integer nwgts,nevents
11 double precision sum_of_wgts12 double precision sum_of_wgts
12 common/to_lhe_analysis/sum_of_wgts,nwgts,nevents13 common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
13 logical lopen14 logical lopen
14 nwgts=nwgt15 nwgts=nwgt
15 sum_of_wgts=0d016 sum_of_wgts=0d0
@@ -23,9 +24,10 @@
23 subroutine analysis_end(xnorm)24 subroutine analysis_end(xnorm)
24 implicit none25 implicit none
25 double precision xnorm26 double precision xnorm
27 logical passed_unwgt
26 integer nwgts,nevents28 integer nwgts,nevents
27 double precision sum_of_wgts29 double precision sum_of_wgts
28 common/to_lhe_analysis/sum_of_wgts,nwgts,nevents30 common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
29 logical lopen31 logical lopen
30 integer npoints32 integer npoints
31 double precision cross_section33 double precision cross_section
@@ -55,9 +57,10 @@
55 double precision p(0:4,nexternal)57 double precision p(0:4,nexternal)
56 double precision wgts(*)58 double precision wgts(*)
57 integer ibody59 integer ibody
60 logical passed_unwgt
58 integer nwgts,nevents61 integer nwgts,nevents
59 double precision sum_of_wgts62 double precision sum_of_wgts
60 common/to_lhe_analysis/sum_of_wgts,nwgts,nevents63 common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
61c64c
62 integer i,j,npart65 integer i,j,npart
63 double precision ran266 double precision ran2
@@ -95,23 +98,33 @@
95 integer npoints98 integer npoints
96 double precision cross_section99 double precision cross_section
97 common /for_FixedOrder_lhe/ cross_section,npoints100 common /for_FixedOrder_lhe/ cross_section,npoints
101 double precision,allocatable :: wwgts(:)
98c --- do the partial unweighting (don't do it when to close to singular102c --- do the partial unweighting (don't do it when to close to singular
99c --- region)103c --- region)
104 if (.not. allocated(wwgts)) then
105 allocate(wwgts(nwgts))
106 endif
107 passed_unwgt=.true.
100 twgt=abs(cross_section)*FO_LHE_weight_ratio !/npoints108 twgt=abs(cross_section)*FO_LHE_weight_ratio !/npoints
101 if(abs(wgts(1)).lt.abs(twgt) .and. xi_i_fks_ev.gt.1d-3 .and.109 if(abs(wgts(1)).lt.abs(twgt) .and. xi_i_fks_ev.gt.1d-3 .and.
102 $ 1d0-y_ij_fks_ev.gt.1d-2)then110 $ 1d0-y_ij_fks_ev.gt.1d-2)then
103 R = ran2()*abs(twgt)111 R = ran2()*abs(twgt)
104 if (R.gt.abs(wgts(1)))then112 if (R.gt.abs(wgts(1)))then
113 passed_unwgt=.false.
105 return114 return
106 else115 else
107 do i=2,nwgts116 do i=2,nwgts
108 wgts(i) = wgts(i)*abs(twgt/wgts(1))117 wwgts(i) = wgts(i)*abs(twgt/wgts(1))
109 enddo118 enddo
110 wgts(1) = sign(twgt,wgts(1))119 wwgts(1) = sign(twgt,wgts(1))
111 endif120 endif
121 else
122 do i=1,nwgts
123 wwgts(i)=wgts(i)
124 enddo
112 endif125 endif
113c --- accumulate the total weights of all the events in the event file126c --- accumulate the total weights of all the events in the event file
114 sum_of_wgts=sum_of_wgts+wgts(1)127 sum_of_wgts=sum_of_wgts+wwgts(1)
115128
116c --- fill the multi-weight common blocks with the scale and PDF129c --- fill the multi-weight common blocks with the scale and PDF
117c --- variation weights130c --- variation weights
@@ -122,12 +135,12 @@
122 do ii=1,nint(scalevarF(0))135 do ii=1,nint(scalevarF(0))
123 do jj=1,nint(scalevarR(0))136 do jj=1,nint(scalevarR(0))
124 i_wgt=i_wgt+1137 i_wgt=i_wgt+1
125 wgtxsecmu(jj,ii,kk)= wgts(i_wgt)138 wgtxsecmu(jj,ii,kk)= wwgts(i_wgt)
126 enddo139 enddo
127 enddo140 enddo
128 else141 else
129 i_wgt=i_wgt+1142 i_wgt=i_wgt+1
130 wgtxsecmu(1,1,kk)= wgts(i_wgt)143 wgtxsecmu(1,1,kk)= wwgts(i_wgt)
131 endif144 endif
132 enddo145 enddo
133 endif146 endif
@@ -136,11 +149,11 @@
136 if (lpdfvar(nn)) then149 if (lpdfvar(nn)) then
137 do n=0,nmemPDF(nn)150 do n=0,nmemPDF(nn)
138 i_wgt=i_wgt+1151 i_wgt=i_wgt+1
139 wgtxsecPDF(n,nn) = wgts(i_wgt)152 wgtxsecPDF(n,nn) = wwgts(i_wgt)
140 enddo153 enddo
141 else154 else
142 i_wgt=i_wgt+1155 i_wgt=i_wgt+1
143 wgtxsecPDF(0,nn) = wgts(i_wgt)156 wgtxsecPDF(0,nn) = wwgts(i_wgt)
144 endif157 endif
145 enddo158 enddo
146 endif159 endif
@@ -187,7 +200,7 @@
187200
188c --- write the event201c --- write the event
189 call write_lhef_event(41,202 call write_lhef_event(41,
190 & npart,IDPRUP,wgts(1),0d0,0d0,0d0,203 & npart,IDPRUP,wwgts(1),0d0,0d0,0d0,
191 & IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)204 & IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
192205
193 201 format(a9,1x,i1,4(1x,i2),2(1x,d14.8),2x,i2,2(1x,i2),5(1x,d14.8))206 201 format(a9,1x,i1,4(1x,i2),2(1x,d14.8),2x,i2,2(1x,i2),5(1x,d14.8))
@@ -196,12 +209,15 @@
196c This we can use for the event grouping! 209c This we can use for the event grouping!
197 subroutine HwU_add_points210 subroutine HwU_add_points
198 implicit none211 implicit none
212 logical passed_unwgt
199 integer nwgts,nevents213 integer nwgts,nevents
200 double precision sum_of_wgts214 double precision sum_of_wgts
201 common/to_lhe_analysis/sum_of_wgts,nwgts,nevents215 common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
202 nevents=nevents+1216 if (passed_unwgt) then
203 write (41,'(a)') '</eventgroup>'217 nevents=nevents+1
204 write (41,'(a)') '<eventgroup>'218 write (41,'(a)') '</eventgroup>'
219 write (41,'(a)') '<eventgroup>'
220 endif
205 end221 end
206222
207c Dummy routines223c Dummy routines
208224
=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f 2017-07-31 08:48:20 +0000
+++ Template/NLO/SubProcesses/check_poles.f 2018-12-20 08:18:55 +0000
@@ -18,7 +18,7 @@
18 include 'genps.inc'18 include 'genps.inc'
19 include 'nexternal.inc'19 include 'nexternal.inc'
20 include 'nFKSconfigs.inc'20 include 'nFKSconfigs.inc'
21 double precision p(0:3, nexternal), prambo(0:3, nexternal)21 double precision p(0:3, nexternal), prambo(0:3,100)
22 double precision p_born(0:3,nexternal-1)22 double precision p_born(0:3,nexternal-1)
23 common/pborn/p_born23 common/pborn/p_born
24 double precision pswgt24 double precision pswgt
@@ -48,7 +48,7 @@
48 include 'coupl.inc'48 include 'coupl.inc'
49 include 'q_es.inc'49 include 'q_es.inc'
50 integer nsqso,MLResArrayDim50 integer nsqso,MLResArrayDim
51 double precision pmass(nexternal), pmass_rambo(nexternal)51 double precision pmass(nexternal), pmass_rambo(100)
52 integer nfail52 integer nfail
53 logical first_time53 logical first_time
54 data first_time/.TRUE./54 data first_time/.TRUE./
5555
=== modified file 'Template/NLO/SubProcesses/combine_plots_FO.sh'
--- Template/NLO/SubProcesses/combine_plots_FO.sh 2013-05-25 12:52:41 +0000
+++ Template/NLO/SubProcesses/combine_plots_FO.sh 2018-12-20 08:18:55 +0000
@@ -1,6 +1,4 @@
1#!/bin/bash1#!/bin/bash
2TD=td_mac_intel
3PSOPEN=open
42
5# $? is the value of last executed command. A call to this function3# $? is the value of last executed command. A call to this function
6# after a failure will cause the program to quit the script4# after a failure will cause the program to quit the script
@@ -15,126 +13,74 @@
15fi13fi
16}14}
1715
18thisdir=`pwd`16function combine {
19if [ -f $thisdir/MADatNLO.top ]17counter=0
20then18i=1
21rm -f $thisdir/MADatNLO*.top
22fi
23make read40
24EXENAME=$thisdir/read40
25echo -n "" > dir19echo -n "" > dir
26counterp=020for plot_file in "${@:2}" ; do
27for p in P* ; do21 counter=$(($counter + 1))
28 cd $thisdir22 echo $plot_file >> dir
29 cd $p23 if [[ $(($counter % 40)) == 0 ]]; then
30 rm -f MADatNLO*24 echo -n "" > dir2
31 rm -f *read40*.out25 echo "40" > dir2
32 echo -n "" > dir
33 counter=0
34 i=1
35 for g in $* ; do
36 if [ ! -f $g/MADatNLO*.top ]
37 then
38 echo "No topdrawer file in $p/$g"
39 else
40 counter=$(($counter + 1))
41 echo $g"/MADatNLO.top" >> dir
42 if [[ $(($counter % 40)) == 0 ]]; then
43 echo -n "" > dir2
44 echo "40" > dir2
45 cat dir >> dir2
46 $EXENAME < dir2
47 teststatus Failure in step 1
48 rm -f dir
49 rm -f dir2
50 mv fort.88 MADatNLO$i\.top
51 mv read40.out "S1read40_"$i\.out
52 i=$(($i + 1))
53 echo -n "" > dir
54 fi
55 fi
56 done
57 if [[ $(($counter % 40)) != 0 ]] ; then
58 echo $(($counter % 40)) > dir2
59 cat dir >> dir226 cat dir >> dir2
60 while [[ $(($counter % 40)) -ne 0 ]]; do27 ./read40 < dir2
61 counter=$[$counter-1]28 teststatus Failure in step 1
62 done
63 $EXENAME < dir2
64 teststatus Failure in step 2
65 rm -f dir29 rm -f dir
66 rm -f dir230 rm -f dir2
67 mv fort.88 MADatNLO$i\.top 31 mv fort.88 $1_$i\.top
68 mv read40.out "S2read40_"$i\.out32 mv read40.out "S1read40_"$1_$i\.out
69 fi33 i=$(($i + 1))
70 34 echo -n "" > dir
71 counter=035 fi
72 for m in MADatNLO*.top ; do36done
73 counter=$[$counter+1]37if [[ $(($counter % 40)) != 0 ]] ; then
74 echo $m >> dir38 echo $(($counter % 40)) > dir2
75 done39 cat dir >> dir2
76 if [[ $counter -gt 1 ]] ; then40 ./read40 < dir2
77 echo $counter > dir241 teststatus Failure in step 2
78 cat dir >> dir2
79 $EXENAME < dir2
80 teststatus Failure in step 3
81 rm -f dir2
82 mv fort.88 MADatNLO.top
83 mv read40.out S3read40.out
84 else
85 mv MADatNLO1.top MADatNLO.top
86 fi
87 rm -f dir
88 cd $thisdir
89 counterp=$(($counterp + 1))
90 echo $p"/MADatNLO.top" >> dir
91 echo $counterp $p "done"
92done
93
94cd $thisdir
95make split40
96./split40
97
98i=0
99#pdir* are created by split
100for p in pdir* ; do
101 i=$(($i+1))
102 $EXENAME < $p
103 teststatus Failure in step 4
104 mv fort.88 MADatNLO$i\.top
105 mv read40.out "S4read40_"$i\.out
106 rm $p
107 echo $i $p "done"
108done
109
110if [[ $i -gt 1 ]] ; then
111 echo $i > dir2
112 for p in MADatNLO*.top ; do
113 echo $p >>dir2
114 done
115 $EXENAME < dir2
116 teststatus Failure in step 5
117 rm -f dir42 rm -f dir
118 rm -f dir243 rm -f dir2
119 mv fort.88 MADatNLO.top 44 mv fort.88 $1_$i\.top
120 mv read40.out S5read40.out45 mv read40.out "S2read40_"$1_$i\.out
46fi
47}
48
49make read40
50rm -f MADatNLO*.top
51
52if [ "$#" -gt 64000 ]; then
53 combine "MADatNLO_A" "$@"
54 combine "MADatNLO_B" MADatNLO_A*.top
55 combine "MADatNLO_C" MADatNLO_B*.top
56 combine "MADatNLO_D" MADatNLO_C*.top
57 mv MADatNLO_D_1.top MADatNLO.top
58elif [ "$#" -gt 1600 ]; then
59 combine "MADatNLO_A" "$@"
60 combine "MADatNLO_B" MADatNLO_A*.top
61 combine "MADatNLO_C" MADatNLO_B*.top
62 mv MADatNLO_C_1.top MADatNLO.top
63elif [ "$#" -gt 40 ]; then
64 combine "MADatNLO_A" "$@"
65 combine "MADatNLO_B" MADatNLO_A*.top
66 mv MADatNLO_B_1.top MADatNLO.top
67elif [ "$#" -gt 1 ]; then
68 combine "MADatNLO_A" "$@"
69 mv MADatNLO_A_1.top MADatNLO.top
121else70else
122 mv MADatNLO1.top MADatNLO.top71 cp $1 MADatNLO.top
123fi72fi
12473
125if [ -f $thisdir/read40.errors ]74
75if [ -f read40.errors ]
126then76then
127rm -f $thisdir/read40.errors77rm -f read40.errors
128fi78fi
129find . -name S"*"read40"*".out -exec fgrep -il "READ40 ERROR READ40 ERROR" {} \; > $thisdir/read40.errors79
13080find . -name S"*"read40"*".out -exec fgrep -il "READ40 ERROR READ40 ERROR" {} \; > read40.errors
131if [ ! -s $thisdir/read40.errors ] ; then81
82if [ ! -s read40.errors ] ; then
132 echo 'no errors found'83 echo 'no errors found'
133else84else
134 echo 'ERRORS! see read40.errors for details'85 echo 'ERRORS! see read40.errors for details'
135fi86fi
136
137#./NLO_Born3.py
138#$TD MADatNLO_combined.top "POSTSCR,ORIENT=3"
139#$PSOPEN MADatNLO_combined.ps
140##td2pdf MADatNLO_con
14187
=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f 2018-04-16 14:08:47 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f 2018-12-20 08:18:55 +0000
@@ -48,7 +48,7 @@
48 include "mint.inc"48 include "mint.inc"
49 integer nhits_in_grids(maxchannels)49 integer nhits_in_grids(maxchannels)
50 real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals50 real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
51 $ ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals51 $ ,ndimmax,maxchannels),ymax_virt(0:maxchannels),ans(nintegrals
52 $ ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals52 $ ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
53 $ ,0:maxchannels),x(ndimmax),itmax_fl53 $ ,0:maxchannels),x(ndimmax),itmax_fl
54 integer ixi_i,iphi_i,iy_ij,vn54 integer ixi_i,iphi_i,iy_ij,vn
@@ -402,9 +402,10 @@
402 double precision xx(ndimmax),vegas_wgt,f(nintegrals),jac,p(0:3402 double precision xx(ndimmax),vegas_wgt,f(nintegrals),jac,p(0:3
403 $ ,nexternal),rwgt,vol,sig,x(99),MC_int_wgt403 $ ,nexternal),rwgt,vol,sig,x(99),MC_int_wgt
404 integer ifl,nFKS_born,nFKS_picked,iFKS,nFKS_min,iamp404 integer ifl,nFKS_born,nFKS_picked,iFKS,nFKS_min,iamp
405 $ ,nFKS_max,izero,ione,itwo,mohdr405 $ ,nFKS_max,izero,ione,itwo,mohdr,i,iran_picked
406 parameter (izero=0,ione=1,itwo=2,mohdr=-100)406 parameter (izero=0,ione=1,itwo=2,mohdr=-100)
407 logical passcuts,passcuts_nbody,passcuts_n1body,sum407 logical passcuts,passcuts_nbody,passcuts_n1body,sum,firsttime
408 data firsttime/.true./
408 external passcuts409 external passcuts
409 integer ini_fin_fks(maxchannels)410 integer ini_fin_fks(maxchannels)
410 common/fks_channels/ini_fin_fks411 common/fks_channels/ini_fin_fks
@@ -431,6 +432,16 @@
431 common /for_applgrid/ iappl432 common /for_applgrid/ iappl
432 double precision wgt_ME_born,wgt_ME_real433 double precision wgt_ME_born,wgt_ME_real
433 common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real434 common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
435 integer ini_fin_fks_map(0:2,0:fks_configs)
436 save ini_fin_fks_map
437 if (firsttime) then
438 firsttime=.false.
439 call setup_ini_fin_fks_map(ini_fin_fks_map)
440 write (*,*) 'initial-final FKS maps:'
441 write (*,*) 0 ,':',ini_fin_fks_map(0,:)
442 write (*,*) 1 ,':',ini_fin_fks_map(1,:)
443 write (*,*) 2 ,':',ini_fin_fks_map(2,:)
444 endif
434 if (ifl.ne.0) then445 if (ifl.ne.0) then
435 write (*,*) 'ERROR ifl not equal to zero in sigint',ifl446 write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
436 stop 1447 stop 1
@@ -453,20 +464,9 @@
453 if (ickkw.eq.-1) H1_factor_virt=0d0464 if (ickkw.eq.-1) H1_factor_virt=0d0
454 if (ickkw.eq.3) call set_FxFx_scale(0,p)465 if (ickkw.eq.3) call set_FxFx_scale(0,p)
455 call update_vegas_x(xx,x)466 call update_vegas_x(xx,x)
456 call get_MC_integer(max(ini_fin_fks(ichan),1),fks_configs467 call get_MC_integer(max(ini_fin_fks(ichan),1)
457 $ ,nFKS_picked,vol)468 $ ,ini_fin_fks_map(ini_fin_fks(ichan),0),iran_picked,vol)
458 if (ini_fin_fks(ichan).eq.1) then469 nFKS_picked=ini_fin_fks_map(ini_fin_fks(ichan),iran_picked)
459 do while (fks_j_d(nFKS_picked).le.nincoming)
460 call get_MC_integer(ini_fin_fks(ichan),fks_configs
461 $ ,nFKS_picked,vol)
462 enddo
463 elseif (ini_fin_fks(ichan).eq.2) then
464 do while (fks_j_d(nFKS_picked).gt.nincoming)
465 call get_MC_integer(ini_fin_fks(ichan),fks_configs
466 $ ,nFKS_picked,vol)
467 enddo
468 endif
469
470 470
471c The nbody contributions471c The nbody contributions
472 if (abrv.eq.'real') goto 11472 if (abrv.eq.'real') goto 11
@@ -503,14 +503,15 @@
503 nbody=.false.503 nbody=.false.
504 if (sum) then504 if (sum) then
505 nFKS_min=1505 nFKS_min=1
506 nFKS_max=fks_configs506 nFKS_max=ini_fin_fks_map(ini_fin_fks(ichan),0)
507 MC_int_wgt=1d0507 MC_int_wgt=1d0
508 else508 else
509 nFKS_min=nFKS_picked509 nFKS_min=iran_picked
510 nFKS_max=nFKS_picked510 nFKS_max=iran_picked
511 MC_int_wgt=1d0/vol511 MC_int_wgt=1d0/vol
512 endif512 endif
513 do iFKS=nFKS_min,nFKS_max513 do i=nFKS_min,nFKS_max
514 iFKS=ini_fin_fks_map(ini_fin_fks(ichan),i)
514 calculatedBorn=.false. 515 calculatedBorn=.false.
515 ! MZ this is a temporary fix for processes without516 ! MZ this is a temporary fix for processes without
516 ! soft singularities associated to the initial state517 ! soft singularities associated to the initial state
@@ -566,11 +567,11 @@
566c Importance sampling for FKS configurations567c Importance sampling for FKS configurations
567 if (sum) then568 if (sum) then
568 call get_wgt_nbody(sig)569 call get_wgt_nbody(sig)
569 call fill_MC_integer(max(ini_fin_fks(ichan),1),nFKS_picked570 call fill_MC_integer(max(ini_fin_fks(ichan),1),iran_picked
570 $ ,abs(sig))571 $ ,abs(sig))
571 else572 else
572 call get_wgt_no_nbody(sig)573 call get_wgt_no_nbody(sig)
573 call fill_MC_integer(max(ini_fin_fks(ichan),1),nFKS_picked574 call fill_MC_integer(max(ini_fin_fks(ichan),1),iran_picked
574 $ ,abs(sig)*vol)575 $ ,abs(sig)*vol)
575 endif576 endif
576577
@@ -595,6 +596,35 @@
595 return596 return
596 end597 end
597 598
599 subroutine setup_ini_fin_FKS_map(ini_fin_FKS_map)
600 implicit none
601 include 'nexternal.inc'
602 include 'nFKSconfigs.inc'
603 include 'fks_info.inc'
604 integer ini_fin_FKS_map(0:2,0:fks_configs),iFKS
605 ini_fin_FKS_map(0,0)=0
606 ini_fin_FKS_map(1,0)=0
607 ini_fin_FKS_map(2,0)=0
608 do iFKS=1,fks_configs
609 ini_fin_FKS_map(0,0)=ini_fin_FKS_map(0,0)+1
610 ini_fin_FKS_map(0,ini_fin_FKS_map(0,0))=iFKS
611 if (fks_j_d(iFKS).le.nincoming .and.
612 $ fks_j_d(iFKS).gt.0) then
613 ini_fin_FKS_map(2,0)=ini_fin_FKS_map(2,0)+1
614 ini_fin_FKS_map(2,ini_fin_FKS_map(2,0))=iFKS
615 elseif (fks_j_d(iFKS).gt.nincoming .and.
616 $ fks_j_d(iFKS).le.nexternal) then
617 ini_fin_FKS_map(1,0)=ini_fin_FKS_map(1,0)+1
618 ini_fin_FKS_map(1,ini_fin_FKS_map(1,0))=iFKS
619 else
620 write (*,*) 'ERROR in setup_ini_fin_FKS_map',fks_j_d(iFKS)
621 $ ,nincoming,iFKS
622 stop 1
623 endif
624 enddo
625 return
626 end
627
598 subroutine get_born_nFKSprocess(nFKS_in,nFKS_out)628 subroutine get_born_nFKSprocess(nFKS_in,nFKS_out)
599 implicit none629 implicit none
600 include 'nexternal.inc'630 include 'nexternal.inc'
601631
=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f 2018-04-16 14:08:47 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f 2018-12-20 08:18:55 +0000
@@ -59,7 +59,7 @@
59 common /event_normalisation/event_norm59 common /event_normalisation/event_norm
60c For MINT:60c For MINT:
61 real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals61 real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
62 $ ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals62 $ ,ndimmax,maxchannels),ymax_virt(0:maxchannels),ans(nintegrals
63 $ ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals63 $ ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
64 $ ,0:maxchannels),x(ndimmax)64 $ ,0:maxchannels),x(ndimmax)
65 integer ixi_i,iphi_i,iy_ij,vn,nhits_in_grids(maxchannels)65 integer ixi_i,iphi_i,iy_ij,vn,nhits_in_grids(maxchannels)
@@ -842,7 +842,7 @@
842 firsttime=.false.842 firsttime=.false.
843c Determines the proc_map that sets which FKS configuration can be843c Determines the proc_map that sets which FKS configuration can be
844c summed explicitly and which by MC-ing.844c summed explicitly and which by MC-ing.
845 call setup_proc_map(sum,proc_map)845 call setup_proc_map(sum,proc_map,ini_fin_fks)
846c For the S-events, we can combine processes when they give identical846c For the S-events, we can combine processes when they give identical
847c processes at the Born. Make sure we check that we get indeed identical847c processes at the Born. Make sure we check that we get indeed identical
848c IRPOC's848c IRPOC's
@@ -869,16 +869,6 @@
869 if (ickkw.eq.3) call set_FxFx_scale(0,p)869 if (ickkw.eq.3) call set_FxFx_scale(0,p)
870 call update_vegas_x(xx,x)870 call update_vegas_x(xx,x)
871 call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)871 call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
872 if (ini_fin_fks.eq.1) then
873 do while (fks_j_d(proc_map(proc_map(0,1),1)).le.nincoming)
874 call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
875 enddo
876 elseif (ini_fin_fks.eq.2) then
877 do while (fks_j_d(proc_map(proc_map(0,1),1)).gt.nincoming)
878 call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
879 enddo
880 endif
881
882c The nbody contributions872c The nbody contributions
883 if (abrv.eq.'real') goto 11873 if (abrv.eq.'real') goto 11
884 nbody=.true.874 nbody=.true.
@@ -1043,7 +1033,7 @@
1043 end1033 end
10441034
10451035
1046 subroutine setup_proc_map(sum,proc_map)1036 subroutine setup_proc_map(sum,proc_map,ini_fin_fks)
1047c Determines the proc_map that sets which FKS configuration can be1037c Determines the proc_map that sets which FKS configuration can be
1048c summed explicitly and which by MC-ing.1038c summed explicitly and which by MC-ing.
1049 implicit none1039 implicit none
@@ -1056,7 +1046,7 @@
1056 logical found_ini1,found_ini2,found_fnl1046 logical found_ini1,found_ini2,found_fnl
1057 integer proc_map(0:fks_configs,0:fks_configs)1047 integer proc_map(0:fks_configs,0:fks_configs)
1058 $ ,j_fks_proc(fks_configs),i_fks_pdg_proc(fks_configs)1048 $ ,j_fks_proc(fks_configs),i_fks_pdg_proc(fks_configs)
1059 $ ,j_fks_pdg_proc(fks_configs),i,sum,j1049 $ ,j_fks_pdg_proc(fks_configs),i,sum,j,ini_fin_fks
1060 integer nFKSprocess1050 integer nFKSprocess
1061 common/c_nFKSprocess/nFKSprocess1051 common/c_nFKSprocess/nFKSprocess
1062 INTEGER IPROC1052 INTEGER IPROC
@@ -1097,6 +1087,8 @@
1097c them in the process map1087c them in the process map
1098 do nFKSprocess=1,fks_configs1088 do nFKSprocess=1,fks_configs
1099 call fks_inc_chooser()1089 call fks_inc_chooser()
1090 if (ini_fin_fks.eq.1 .and. j_fks.le.nincoming) cycle
1091 if (ini_fin_fks.eq.2 .and. j_fks.gt.nincoming) cycle
1100 if (need_color_links.or.need_charge_links) then1092 if (need_color_links.or.need_charge_links) then
1101 proc_map(0,0)=proc_map(0,0)+11093 proc_map(0,0)=proc_map(0,0)+1
1102 proc_map(proc_map(0,0),0)=proc_map(proc_map(0,0),0)+11094 proc_map(proc_map(0,0),0)=proc_map(proc_map(0,0),0)+1
@@ -1156,6 +1148,8 @@
1156c gluons splitting1148c gluons splitting
1157 do nFKSprocess=1,fks_configs1149 do nFKSprocess=1,fks_configs
1158 call fks_inc_chooser()1150 call fks_inc_chooser()
1151 if (ini_fin_fks.eq.1 .and. j_fks.le.nincoming) cycle
1152 if (ini_fin_fks.eq.2 .and. j_fks.gt.nincoming) cycle
1159 if (.not.(need_color_links.or.need_charge_links)) then1153 if (.not.(need_color_links.or.need_charge_links)) then
1160 if (j_fks.eq.1 .and. found_ini1) then1154 if (j_fks.eq.1 .and. found_ini1) then
1161 do i=1,proc_map(0,0)1155 do i=1,proc_map(0,0)
11621156
=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f 2018-04-16 14:08:47 +0000
+++ Template/NLO/SubProcesses/fks_singular.f 2018-12-20 08:18:55 +0000
@@ -1115,8 +1115,6 @@
1115 double precision xi_i_fks_ev,y_ij_fks_ev1115 double precision xi_i_fks_ev,y_ij_fks_ev
1116 double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)1116 double precision p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
1117 common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt1117 common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
1118 double precision xi_i_hat_ev,xi_i_hat_cnt(-2:2)
1119 common /cxi_i_hat/xi_i_hat_ev,xi_i_hat_cnt
1120 integer i_fks,j_fks1118 integer i_fks,j_fks
1121 common/fks_indices/i_fks,j_fks1119 common/fks_indices/i_fks,j_fks
1122 double precision xi_i_fks_cnt(-2:2)1120 double precision xi_i_fks_cnt(-2:2)
@@ -1249,7 +1247,7 @@
1249 endif1247 endif
12501248
1251c f_* multiplication factors for real-emission, soft counter, ... etc. 1249c f_* multiplication factors for real-emission, soft counter, ... etc.
1252 prefact=1d0/xi_i_hat_ev/(1-y_ij_fks_ev)1250 prefact=xinorm_ev/xi_i_fks_ev/(1-y_ij_fks_ev)
1253 f_r=prefact*jac_ev*enhance_real*fkssymmetryfactor*vegas_wgt1251 f_r=prefact*jac_ev*enhance_real*fkssymmetryfactor*vegas_wgt
1254 f_MC_S=f_r1252 f_MC_S=f_r
1255 f_MC_H=f_r1253 f_MC_H=f_r
@@ -1736,11 +1734,16 @@
1736c coefficients for PDF and scale computations. 1734c coefficients for PDF and scale computations.
1737 use weight_lines1735 use weight_lines
1738 implicit none1736 implicit none
1739 integer i,j1737 include 'orders.inc'
1738 include 'mint.inc'
1739 integer orders(nsplitorders)
1740 integer i,j,iamp
1741 logical virt_found
1740 double precision bias1742 double precision bias
1741 character*7 event_norm1743 character*7 event_norm
1742 common /event_normalisation/event_norm1744 common /event_normalisation/event_norm
1743 double precision virt_wgt_mint,born_wgt_mint1745 double precision virt_wgt_mint(0:n_ave_virt),
1746 & born_wgt_mint(0:n_ave_virt)
1744 common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint1747 common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
1745c Set the bias_wgt to 1 in case we do not have to do any biassing1748c Set the bias_wgt to 1 in case we do not have to do any biassing
1746 if (event_norm(1:4).ne.'bias') then1749 if (event_norm(1:4).ne.'bias') then
@@ -1749,6 +1752,7 @@
1749 enddo1752 enddo
1750 return1753 return
1751 endif1754 endif
1755 virt_found=.false.
1752c loop over all contributions1756c loop over all contributions
1753 do i=1,icontr1757 do i=1,icontr
1754 if (itype(i).eq.1) then1758 if (itype(i).eq.1) then
@@ -1771,9 +1775,15 @@
1771 do j=1,31775 do j=1,3
1772 wgt(j,i)=wgt(j,i)*bias_wgt(i)1776 wgt(j,i)=wgt(j,i)*bias_wgt(i)
1773 enddo1777 enddo
1774 if (itype(i).eq.14) then1778 if (itype(i).eq.14 .and. .not.virt_found) then
1775 virt_wgt_mint=virt_wgt_mint*bias_wgt(i)1779 virt_found=.true.
1776 born_wgt_mint=born_wgt_mint*bias_wgt(i)1780 virt_wgt_mint(0)=virt_wgt_mint(0)*bias_wgt(i)
1781 born_wgt_mint(0)=born_wgt_mint(0)*bias_wgt(i)
1782 do iamp=1,amp_split_size
1783 call amp_split_pos_to_orders(iamp, orders)
1784 virt_wgt_mint(iamp)=virt_wgt_mint(iamp)*bias_wgt(i)
1785 born_wgt_mint(iamp)=born_wgt_mint(iamp)*bias_wgt(i)
1786 enddo
1777 endif1787 endif
1778 enddo1788 enddo
1779 return1789 return
@@ -3805,7 +3815,7 @@
3805 if (iextra_cnt.gt.0)3815 if (iextra_cnt.gt.0)
3806 1 call extra_cnt(p_born, iextra_cnt, ans_extra_cnt)3816 1 call extra_cnt(p_born, iextra_cnt, ans_extra_cnt)
3807 call AP_reduced(j_type,i_type,ch_j,ch_i,t,z,ap)3817 call AP_reduced(j_type,i_type,ch_j,ch_i,t,z,ap)
3808 call Qterms_reduced_timelike(j_type,i_type,ch_i,ch_j,t,z,Q)3818 call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,t,z,Q)
3809 wgt=0d03819 wgt=0d0
3810 do iord = 1, nsplitorders3820 do iord = 1, nsplitorders
3811 if (.not.split_type(iord) .or. (iord.ne.qed_pos .and.3821 if (.not.split_type(iord) .or. (iord.ne.qed_pos .and.
38123822
=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
--- Template/NLO/SubProcesses/genps_fks.f 2018-01-29 05:41:33 +0000
+++ Template/NLO/SubProcesses/genps_fks.f 2018-12-20 08:18:55 +0000
@@ -2557,7 +2557,7 @@
2557 integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)2557 integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
2558 double precision cBW_mass(-1:1,-nexternal:-1),2558 double precision cBW_mass(-1:1,-nexternal:-1),
2559 & cBW_width(-1:1,-nexternal:-1)2559 & cBW_width(-1:1,-nexternal:-1)
2560 double precision s_mass(-nexternal:-1)2560 double precision s_mass(-nexternal:nexternal)
2561 common/to_phase_space_s_channel/s_mass2561 common/to_phase_space_s_channel/s_mass
2562c Common block with granny information2562c Common block with granny information
2563 logical granny_is_res2563 logical granny_is_res
25642564
=== modified file 'Template/NLO/SubProcesses/makefile'
--- Template/NLO/SubProcesses/makefile 2017-06-15 16:22:23 +0000
+++ Template/NLO/SubProcesses/makefile 2018-12-20 08:18:55 +0000
@@ -7,9 +7,6 @@
7# Files for the read40 combiner of top drawer files7# Files for the read40 combiner of top drawer files
8READ40=read40.o8READ40=read40.o
99
10# Files for the split40 splitter of top drawer files
11SPLIT40=split40.o
12
13# Files to collect all the events in the separate integration channels into a single file10# Files to collect all the events in the separate integration channels into a single file
14COLLECT_EVENTS=collect_events.o handling_lhe_events.o fill_MC_mshell.o11COLLECT_EVENTS=collect_events.o handling_lhe_events.o fill_MC_mshell.o
1512
@@ -21,9 +18,6 @@
21read40: $(READ40)18read40: $(READ40)
22 $(FC) $(LDFLAGS) -o read40 $(READ40)19 $(FC) $(LDFLAGS) -o read40 $(READ40)
2320
24split40: $(SPLIT40)
25 $(FC) $(LDFLAGS) -o split40 $(SPLIT40)
26
27collect_events: $(COLLECT_EVENTS) $(LIBS)21collect_events: $(COLLECT_EVENTS) $(LIBS)
28 $(FC) -o collect_events $(COLLECT_EVENTS) $(LINKLIBS) $(LDFLAGS)22 $(FC) -o collect_events $(COLLECT_EVENTS) $(LINKLIBS) $(LDFLAGS)
29 rm handling_lhe_events.o23 rm handling_lhe_events.o
3024
=== modified file 'Template/NLO/SubProcesses/mint-integrator2.f'
--- Template/NLO/SubProcesses/mint-integrator2.f 2017-07-31 08:48:20 +0000
+++ Template/NLO/SubProcesses/mint-integrator2.f 2018-12-20 08:18:55 +0000
@@ -77,7 +77,7 @@
77 $ ,ans(nintegrals,0:maxchannels),unc(nintegrals,0:maxchannels)77 $ ,ans(nintegrals,0:maxchannels),unc(nintegrals,0:maxchannels)
78 $ ,ans3(nintegrals,3),unc3(nintegrals,3),ans_l3(nintegrals)78 $ ,ans3(nintegrals,3),unc3(nintegrals,3),ans_l3(nintegrals)
79 $ ,unc_l3(nintegrals),chi2_l3(nintegrals),vol_chan,error_virt79 $ ,unc_l3(nintegrals),chi2_l3(nintegrals),vol_chan,error_virt
80 real * 8 xint_virt(maxchannels),ymax_virt(maxchannels)80 real * 8 xint_virt(maxchannels),ymax_virt(0:maxchannels)
81 $ ,ans_chan(0:maxchannels)81 $ ,ans_chan(0:maxchannels)
82 real * 8 x(ndimmax),vol82 real * 8 x(ndimmax),vol
83 real * 8 xacc(0:nintervals,ndimmax,maxchannels)83 real * 8 xacc(0:nintervals,ndimmax,maxchannels)
@@ -216,7 +216,7 @@
216 ans_chan(0)=ans_chan(0)+ans(1,kchan)216 ans_chan(0)=ans_chan(0)+ans(1,kchan)
217 enddo217 enddo
218 endif218 endif
219 cross_section=ans_chan(0)219 cross_section=ans_chan(0) * wgt_mult
220 nit=0220 nit=0
221 nit_included=0221 nit_included=0
222 do i=1,nintegrals222 do i=1,nintegrals
@@ -1031,8 +1031,8 @@
1031 integer ndim,imode1031 integer ndim,imode
1032 include "mint.inc"1032 include "mint.inc"
1033 real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)1033 real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)
1034 $ ,ymax(nintervals,ndimmax,maxchannels),ymax_virt(maxchannels)1034 $ ,ymax(nintervals,ndimmax,maxchannels)
1035 $ ,x(ndimmax)1035 $ ,ymax_virt(0:maxchannels),x(ndimmax)
1036 real * 8 dx(ndimmax),xx(ndimmax),vol_chan,dummy1036 real * 8 dx(ndimmax),xx(ndimmax),vol_chan,dummy
1037 integer icell(ndimmax),ncell(ndimmax),ncell_virt1037 integer icell(ndimmax),ncell(ndimmax),ncell_virt
1038 integer ifold(ndimmax),kfold(ndimmax)1038 integer ifold(ndimmax),kfold(ndimmax)
10391039
=== modified file 'Template/NLO/SubProcesses/montecarlocounter.f'
--- Template/NLO/SubProcesses/montecarlocounter.f 2018-04-16 14:08:47 +0000
+++ Template/NLO/SubProcesses/montecarlocounter.f 2018-12-20 08:18:55 +0000
@@ -1,3231 +1,3259 @@
1 subroutine set_mc_matrices
2 implicit none
3 include "genps.inc"
4 include 'nexternal.inc'
5 include "born_nhel.inc"
6c Nexternal is the number of legs (initial and final) al NLO, while max_bcol
7c is the number of color flows at Born level
8 integer i,j,k,l,k0,mothercol(2),i1(2)
9 integer idup(nexternal-1,maxproc)
10 integer mothup(2,nexternal-1,maxproc)
11 integer icolup(2,nexternal-1,max_bcol)
12 include 'born_leshouche.inc'
13 integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
14 common /MC_info/ ipartners,colorflow
15 integer i_fks,j_fks
16 common/fks_indices/i_fks,j_fks
17 integer fksfather
18 common/cfksfather/fksfather
19 logical notagluon,found
20 common/cnotagluon/notagluon
21 integer nglu,nsngl
22 logical isspecial,isspecial0
23 common/cisspecial/isspecial
24 logical spec_case
25
26 include 'orders.inc'
27 logical split_type(nsplitorders)
28 common /c_split_type/split_type
29 double precision particle_charge(nexternal)
30 common /c_charges/particle_charge
31
32c
33 logical is_leading_cflow(max_bcol)
34 integer num_leading_cflows
35 common/c_leading_cflows/is_leading_cflow,num_leading_cflows
36 include 'born_conf.inc'
37 include 'born_coloramps.inc'
38c
39 ipartners(0)=0
40 do i=1,nexternal-1
41 colorflow(i,0)=0
42 enddo
43
44C What follows is true for QCD-type splittings.
45C For QED-type splittings, ipartner is simply all the charged particles
46C in the event except for FKSfather. In this case, all the born color
47C flows are allowed
48
49c ipartners(0): number of particles that can be colour or anticolour partner
50c of the father, the Born-level particle to which i_fks and j_fks are
51c attached. If one given particle is the colour/anticolour partner of
52c the father in more than one colour flow, it is counted only once
53c in ipartners(0)
54c ipartners(i), 1<=i<=nexternal-1: the label (according to Born-level
55c labelling) of the i^th colour partner of the father
56c
57c colorflow(i,0), 1<=i<=nexternal-1: number of colour flows in which
58c the particle ipartners(i) is a colour partner of the father
59c colorflow(i,j): the actual label (according to born_leshouche.inc)
60c of the j^th colour flow in which the father and ipartners(i) are
61c colour partners
62c
63c Example: in the process q(1) qbar(2) -> g(3) g(4), the two color flows are
64c
65c j=1 i icolup(1) icolup(2) j=2 i icolup(1) icolup(2)
66c 1 500 0 1 500 0
67c 2 0 501 2 0 501
68c 3 500 502 3 502 501
69c 4 502 501 4 500 502
70c
71c and if one fixes for example fksfather=3, then ipartners = 3
72c while colorflow = 0 0 0 1
73c 1 1 0 4
74c 2 1 2 2
75c 1 2 0
76c
77
78 fksfather=min(i_fks,j_fks)
79
80c isspecial will be set equal to .true. only if the father is a gluon,
81c and another gluon will be found which is connected to it by both
82c colour and anticolour
83 isspecial=.false.
84c
85 if (split_type(qcd_pos)) then
86 ! identify the color partners
87c consider only leading colour flows
88 num_leading_cflows=0
89 do i=1,max_bcol
90 is_leading_cflow(i)=.false.
91 do j=1,mapconfig(0)
92 if(icolamp(i,j,1))then
93 is_leading_cflow(i)=.true.
94 num_leading_cflows=num_leading_cflows+1
95 exit
96 endif
97 enddo
98 enddo
99c
100 do i=1,max_bcol
101 if(.not.is_leading_cflow(i))cycle
102c Loop over Born-level colour flows
103 isspecial0=.false.
104c nglu and nsngl are the number of gluons (except for the father) and of
105c colour singlets in the Born process, according to the information
106c stored in ICOLUP
107 nglu=0
108 nsngl=0
109 mothercol(1)=ICOLUP(1,fksfather,i)
110 mothercol(2)=ICOLUP(2,fksfather,i)
111 notagluon=(mothercol(1).eq.0 .or. mothercol(2).eq.0)
112c
113 do j=1,nexternal-1
114c Loop over Born-level particles; j is the possible colour partner of father,
115c and whether this is the case is determined inside this loop
116 if (j.ne.fksfather) then
117c Skip father (it cannot be its own colour partner)
118 if(ICOLUP(1,j,i).eq.0.and.ICOLUP(2,j,i).eq.0)
119 # nsngl=nsngl+1
120 if(ICOLUP(1,j,i).ne.0.and.ICOLUP(2,j,i).ne.0)
121 # nglu=nglu+1
122 if ( (j.le.nincoming.and.fksfather.gt.nincoming) .or.
123 # (j.gt.nincoming.and.fksfather.le.nincoming) ) then
124c father and j not both in the initial or in the final state -- connect
125c colour (1) with colour (i1(1)), and anticolour (2) with anticolour (i1(2))
126 i1(1)=1
127 i1(2)=2
128 else
129c father and j both in the initial or in the final state -- connect
130c colour (1) with anticolour (i1(2)), and anticolour (2) with colour (i1(1))
131 i1(1)=2
132 i1(2)=1
133 endif
134 do l=1,2
135c Loop over colour and anticolour of father
136 found=.false.
137 if( ICOLUP(i1(l),j,i).eq.mothercol(l) .and.
138 & ICOLUP(i1(l),j,i).ne.0 ) then
139c When ICOLUP(i1(l),j,i) = mothercol(l), the colour (if i1(l)=1) or
140c the anticolour (if i1(l)=2) of particle j is connected to the
141c colour (if l=1) or the anticolour (if l=2) of the father
142 k0=-1
143 do k=1,ipartners(0)
144c Loop over previously-found colour/anticolour partners of father
145 if(ipartners(k).eq.j)then
146 if(found)then
147c Safety measure: if this condition is met, it means that there exist
148c k1 and k2 such that ipartners(k1)=ipartners(k2). This is thus a bug,
149c since ipartners() is the list of possible partners of father, where each
150c Born-level particle must appears at most once
151 write(*,*)'Error #1 in set_matrices'
152 write(*,*)i,j,l,k
153 stop
154 endif
155 found=.true.
156 k0=k
157 endif
158 enddo
159 if (.not.found) then
160 ipartners(0)=ipartners(0)+1
161 ipartners(ipartners(0))=j
162 k0=ipartners(0)
163 endif
164c At this point, k0 is the k0^th colour/anticolour partner of father.
165c Therefore, ipartners(k0)=j
166 if(k0.le.0.or.ipartners(k0).ne.j)then
167 write(*,*)'Error #2 in set_matrices'
168 write(*,*)i,j,l,k0,ipartners(k0)
169 stop
170 endif
171 spec_case=l.eq.2 .and. colorflow(k0,0).ge.1 .and.
172 & colorflow(k0,colorflow(k0,0)).eq.i
173 if (.not.spec_case)then
174c Increase by one the number of colour flows in which the father is
175c (anti)colour-connected with its k0^th partner (according to the
176c list defined by ipartners)
177 colorflow(k0,0)=colorflow(k0,0)+1
178c Store the label of the colour flow thus found
179 colorflow(k0,colorflow(k0,0))=i
180 elseif (spec_case)then
181c Special case: father and ipartners(k0) are both gluons, connected
182c by colour AND anticolour: the number of colour flows was overcounted
183c by one unit, so decrease it
184 if( notagluon .or.
185 & ICOLUP(i1(1),j,i).eq.0 .or.
186 & ICOLUP(i1(2),j,i).eq.0 )then
187 write(*,*)'Error #3 in set_matrices'
188 write(*,*)i,j,l,k0,i1(1),i1(2)
189 stop
190 endif
191 colorflow(k0,colorflow(k0,0))=i
192 isspecial0=.true.
193 endif
194 endif
195 enddo
196 endif
197 enddo
198 if( ((nglu+nsngl).gt.(nexternal-2)) .or.
199 # (isspecial0.and.(nglu+nsngl).ne.(nexternal-2)) )then
200 write(*,*)'Error #4 in set_matrices'
201 write(*,*)isspecial0,nglu,nsngl
202 stop
203 endif
204 isspecial=isspecial.or.isspecial0
205 enddo
206
207 else if (split_type(qed_pos)) then
208 ! do nothing, the partner will be assigned at run-time
209 ! (it is kinematics-dependent)
210 continue
211 endif
212 call check_mc_matrices
213 return
214 end
215
216
217
218 subroutine check_mc_matrices
219 implicit none
220 include "nexternal.inc"
221 include "born_nhel.inc"
222c include "fks.inc"
223 integer fks_j_from_i(nexternal,0:nexternal)
224 & ,particle_type(nexternal),pdg_type(nexternal)
225 common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
226 integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
227 common /MC_info/ ipartners,colorflow
228 integer i,j,ipart,iflow,ntot,ithere(1000)
229 integer fksfather
230 common/cfksfather/fksfather
231 logical notagluon
232 common/cnotagluon/notagluon
233 logical isspecial
234 common/cisspecial/isspecial
235
236 include 'orders.inc'
237 logical split_type(nsplitorders)
238 common /c_split_type/split_type
239
240 logical is_leading_cflow(max_bcol)
241 integer num_leading_cflows
242 common/c_leading_cflows/is_leading_cflow,num_leading_cflows
243c
244 if(ipartners(0).gt.nexternal-1)then
245 write(*,*)'Error #1 in check_mc_matrices',ipartners(0)
246 stop
247 endif
248c
249
250 if (split_type(QCD_pos)) then
251 ! these tests only apply for QCD-type splittings
252 do i=1,ipartners(0)
253 ipart=ipartners(i)
254 if( ipart.eq.fksfather .or.
255 # ipart.le.0 .or. ipart.gt.nexternal-1 .or.
256 # ( abs(particle_type(ipart)).ne.3 .and.
257 # particle_type(ipart).ne.8 ) )then
258 write(*,*)'Error #2 in check_mc_matrices',i,ipart,
259 # particle_type(ipart)
260 stop
261 endif
262 enddo
263c
264 do i=1,nexternal-1
265 ithere(i)=1
266 enddo
267 do i=1,ipartners(0)
268 ipart=ipartners(i)
269 ithere(ipart)=ithere(ipart)-1
270 if(ithere(ipart).lt.0)then
271 write(*,*)'Error #3 in check_mc_matrices',i,ipart
272 stop
273 endif
274 enddo
275c
276c ntot is the total number of colour plus anticolour partners of father
277 ntot=0
278 do i=1,ipartners(0)
279 ntot=ntot+colorflow(i,0)
280c
281 if( colorflow(i,0).le.0 .or.
282 # colorflow(i,0).gt.max_bcol )then
283 write(*,*)'Error #4 in check_mc_matrices',i,colorflow(i,0)
284 stop
285 endif
286c
287 do j=1,max_bcol
288 ithere(j)=1
289 enddo
290 do j=1,colorflow(i,0)
291 iflow=colorflow(i,j)
292 ithere(iflow)=ithere(iflow)-1
293 if(ithere(iflow).lt.0)then
294 write(*,*)'Error #5 in check_mc_matrices',i,j,iflow
295 stop
296 endif
297 enddo
298c
299 enddo
300c
301 if( (notagluon.and.ntot.ne.num_leading_cflows) .or.
302 # ( (.not.notagluon).and.
303 # ( (.not.isspecial).and.ntot.ne.(2*num_leading_cflows) .or.
304 # (isspecial.and.ntot.ne.num_leading_cflows) ) ) )then
305 write(*,*)'Error #6 in check_mc_matrices',
306 # notagluon,ntot,num_leading_cflows,max_bcol
307 stop
308 endif
309c
310 if(num_leading_cflows.gt.max_bcol)then
311 write(*,*)'Error #7 in check_mc_matrices',
312 # num_leading_cflows,max_bcol
313 stop
314 endif
315
316 else if (split_type(QED_pos)) then
317 ! write here possible checks for QED-type splittings
318 continue
319 endif
320 return
321 end
322
323
324 subroutine find_ipartner_QED(pp, nofpartners)
325 implicit none
326 include 'nexternal.inc'
327 double precision pp(0:3, nexternal)
328 integer nofpartners
329
330 integer fks_j_from_i(nexternal,0:nexternal)
331 & ,particle_type(nexternal),pdg_type(nexternal)
332 common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
333 double precision particle_charge(nexternal)
334 common /c_charges/particle_charge
335
336 integer i_fks,j_fks,fksfather
337 common/fks_indices/i_fks,j_fks
338 double precision pmass(nexternal)
339 double precision zero
340 parameter (zero=0d0)
341
342 include 'genps.inc'
343 include "born_nhel.inc"
344 integer idup(nexternal-1,maxproc)
345 integer mothup(2,nexternal-1,maxproc)
346 integer icolup(2,nexternal-1,max_bcol)
347 include 'born_leshouche.inc'
348 include 'coupl.inc'
349 integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
350 common /MC_info/ ipartners,colorflow
351c
352c Shower MonteCarlo
353c
354 character*10 shower_mc
355 common /cMonteCarloType/shower_mc
356
357 logical found
358 logical same_state
359 double precision ppmin, ppnow
360 integer partner
361 integer i,j
362 double precision chargeprod
363 double precision dot
364
365 include 'pmass.inc'
366
367 found=.false.
368 ppmin=1d99
369 fksfather=min(i_fks,j_fks)
370
371 if (shower_mc.eq.'PYTHIA8') then
372 ! this should follow what is done in TimeShower::setupQEDdip
373 ! first, look for the lowest-mass same- (opposite-)flavour pair of
374 ! particles in the opposite (same) state of the system
375 do j=1,nexternal
376 if (j.ne.fksfather.and.j.ne.i_fks) then
377 same_state = (j.gt.nincoming.and.fksfather.gt.nincoming).or.
378 $ (j.le.nincoming.and.fksfather.le.nincoming)
379
380 if ((pdg_type(j).eq.pdg_type(fksfather).and..not.same_state).or.
381 $ (pdg_type(j).eq.-pdg_type(fksfather).and.same_state)) then
382
383 ppnow=dot(pp(0,fksfather),pp(0,j)) - pmass(fksfather)*pmass(j)
384 if (ppnow.lt.ppmin) then
385 found=.true.
386 partner=j
387 endif
388 endif
389 endif
390 enddo
391
392 ! if no partner has been found, then look for the
393 ! lowest-mass/chargeprod pair
394 if (.not.found) then
395 do j=1,nexternal
396 if (j.ne.fksfather.and.j.ne.i_fks) then
397 if (particle_charge(fksfather).ne.0d0.and.particle_charge(j).ne.0d0) then
398 ppnow=dot(pp(0,fksfather),pp(0,j)) - pmass(fksfather)*pmass(j) /
399 $ (particle_charge(fksfather) * particle_charge(j))
400 if (ppnow.lt.ppmin) then
401 found=.true.
402 partner=j
403 endif
404 endif
405 endif
406 enddo
407 endif
408
409 ! if no partner has been found, then look for the
410 ! lowest-mass pair
411 if (.not.found) then
412 do j=1,nexternal
413 if (j.ne.fksfather.and.j.ne.i_fks) then
414 ppnow=dot(pp(0,fksfather),pp(0,j)) - pmass(fksfather)*pmass(j)
415 if (ppnow.lt.ppmin) then
416 found=.true.
417 partner=j
418 endif