Index: trunk/omega/tests/Makefile.am =================================================================== --- trunk/omega/tests/Makefile.am (revision 8494) +++ trunk/omega/tests/Makefile.am (revision 8495) @@ -1,1063 +1,1064 @@ # Makefile.am -- Makefile for O'Mega within and without WHIZARD ## ## Process this file with automake to produce Makefile.in ## ######################################################################## # # Copyright (C) 1999-2021 by # Wolfgang Kilian # Thorsten Ohl # Juergen Reuter # with contributions from # cf. main AUTHORS file # # WHIZARD is free software; you can redistribute it and/or modify it # under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2, or (at your option) # any later version. # # WHIZARD is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. # ######################################################################## SUBDIRS = UFO DIST_SUBDIRS = UFO # OMEGA_SPLIT = -target:single_function OMEGA_SPLIT = -target:split_function 10 # OMEGA_SPLIT = -target:split_module 10 # OMEGA_SPLIT = -target:split_file 10 OMEGA_QED = $(top_builddir)/omega/bin/omega_QED$(OCAML_NATIVE_EXT) OMEGA_QED_OPTS = $(OMEGA_SPLIT) -target:parameter_module parameters_QED OMEGA_QCD = $(top_builddir)/omega/bin/omega_QCD$(OCAML_NATIVE_EXT) OMEGA_QCD_OPTS = $(OMEGA_SPLIT) -target:parameter_module parameters_QCD OMEGA_SYM = $(top_builddir)/omega/bin/omega_SYM$(OCAML_NATIVE_EXT) OMEGA_SYM_OPTS = $(OMEGA_SPLIT) -target:parameter_module parameters_SYM OMEGA_SM = $(top_builddir)/omega/bin/omega_SM$(OCAML_NATIVE_EXT) OMEGA_SM_OPTS = $(OMEGA_SPLIT) -target:parameter_module parameters_SM OMEGA_SM_CKM = $(top_builddir)/omega/bin/omega_SM_CKM$(OCAML_NATIVE_EXT) OMEGA_SM_Higgs = $(top_builddir)/omega/bin/omega_SM_Higgs$(OCAML_NATIVE_EXT) OMEGA_THDM = $(top_builddir)/omega/bin/omega_THDM$(OCAML_NATIVE_EXT) OMEGA_THDM_CKM = $(top_builddir)/omega/bin/omega_THDM_CKM$(OCAML_NATIVE_EXT) OMEGA_HSExt = $(top_builddir)/omega/bin/omega_HSExt$(OCAML_NATIVE_EXT) OMEGA_Zprime = $(top_builddir)/omega/bin/omega_Zprime$(OCAML_NATIVE_EXT) OMEGA_SM_top_anom = $(top_builddir)/omega/bin/omega_SM_top_anom$(OCAML_NATIVE_EXT) OMEGA_SM_top_anom_OPTS = $(OMEGA_SPLIT) -target:parameter_module parameters_SM_top_anom OMEGA_UFO = $(top_builddir)/omega/bin/omega_UFO$(OCAML_NATIVE_EXT) OMEGA_UFO_MAJORANA = \ $(top_builddir)/omega/bin/omega_UFO_Majorana$(OCAML_NATIVE_EXT) OMEGA_UFO_OPTS = -target:parameter_module parameters_UFO OMEGA_UFO_PATH = $(top_srcdir)/omega/tests/UFO OMEGA_XXX = $(top_builddir)/omega/bin/omega_%%%$(OCAML_NATIVE_EXT) OMEGA_XXX_OPTS = -target:parameter_module parameters_%%% OMEGA_UFO_XXX_OPTS = \ "-model:UFO_dir $(top_srcdir)/omega/tests/UFO/%%%/ -model:exec" OMEGA_XXX_MAJORANA = \ $(top_builddir)/omega/bin/omega_%%%_Majorana$(OCAML_NATIVE_EXT) OMEGA_XXX_MAJORANA_LEGACY = \ $(top_builddir)/omega/bin/omega_%%%_Majorana_legacy$(OCAML_NATIVE_EXT) OMEGA_QED_VM = $(top_builddir)/omega/bin/omega_QED_VM$(OCAML_NATIVE_EXT) OMEGA_QCD_VM = $(top_builddir)/omega/bin/omega_QCD_VM$(OCAML_NATIVE_EXT) OMEGA_SM_VM = $(top_builddir)/omega/bin/omega_SM_VM$(OCAML_NATIVE_EXT) OMEGA_SM_CKM_VM = $(top_builddir)/omega/bin/omega_SM_CKM_VM$(OCAML_NATIVE_EXT) OMEGA_THDM_VM = $(top_builddir)/omega/bin/omega_THDM_VM$(OCAML_NATIVE_EXT) OMEGA_THDM_CKM_VM = $(top_builddir)/omega/bin/omega_THDM_CKM_VM$(OCAML_NATIVE_EXT) OMEGA_HSExt_VM = $(top_builddir)/omega/bin/omega_HSExt_VM$(OCAML_NATIVE_EXT) OMEGA_Zprime_VM = $(top_builddir)/omega/bin/omega_Zprime_VM$(OCAML_NATIVE_EXT) OMEGA_SM_Higgs_VM = $(top_builddir)/omega/bin/omega_SM_Higgs_VM$(OCAML_NATIVE_EXT) OMEGA_XXX_VM = $(top_builddir)/omega/bin/omega_%%%_VM$(OCAML_NATIVE_EXT) OMEGA_XXX_VM_PARAMS_OPTS = -params -target:parameter_module_external \ parameters_%%% -target:wrapper_module %% -target:bytecode_file % AM_FCFLAGS = -I$(top_builddir)/omega/src AM_LDFLAGS = ######################################################################## ## Default Fortran compiler options ## OpenMP if FC_USE_OPENMP AM_FCFLAGS += $(FCFLAGS_OPENMP) AM_TESTS_ENVIRONMENT = \ export OMP_NUM_THREADS=1; endif ######################################################################## TESTS = XFAIL_TESTS = EXTRA_PROGRAMS = EXTRA_DIST = ######################################################################## include $(top_srcdir)/omega/src/Makefile.ocaml if OCAML_AVAILABLE OCAMLFLAGS += -I $(top_builddir)/omega/src OMEGA_CORE = $(top_builddir)/omega/src/omega_core.cmxa OMEGA_MODELS = $(top_builddir)/omega/src/omega_models.cmxa TESTS += omega_unit EXTRA_PROGRAMS += omega_unit omega_unit_SOURCES = omega_unit.ml omega_unit: $(OMEGA_CORE) omega_unit.cmx @if $(AM_V_P); then :; else echo " OCAMLOPT " $@; fi $(AM_V_at)$(OCAMLOPT) $(OCAMLFLAGS) $(OCAMLOPTFLAGS) -o omega_unit \ unix.cmxa $(OMEGA_CORE) omega_unit.cmx omega_unit.cmx: omega_unit.ml omega_unit.cmx: $(OMEGA_CORE) endif ######################################################################## KINDS = $(top_builddir)/omega/src/kinds.lo TESTS += test_omega95 test_omega95_bispinors EXTRA_PROGRAMS += test_omega95 test_omega95_bispinors test_omega95_SOURCES = test_omega95.f90 omega_testtools.f90 test_omega95_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la test_omega95_bispinors_SOURCES = test_omega95_bispinors.f90 omega_testtools.f90 test_omega95_bispinors_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la test_omega95.o test_omega95_bispinors.o: omega_testtools.o if NOWEB_AVAILABLE test_omega95.f90: $(top_srcdir)/omega/src/omegalib.nw $(NOTANGLE) -R[[$@]] $< | $(CPIF) $@ test_omega95_bispinors.f90: $(top_srcdir)/omega/src/omegalib.nw $(NOTANGLE) -R[[$@]] $< | $(CPIF) $@ omega_testtools.f90: $(top_srcdir)/omega/src/omegalib.nw $(NOTANGLE) -R[[$@]] $< | $(CPIF) $@ endif NOWEB_AVAILABLE ######################################################################## if OCAML_AVAILABLE TESTS += test_qed_eemm EXTRA_PROGRAMS += test_qed_eemm test_qed_eemm_SOURCES = test_qed_eemm.f90 parameters_QED.f90 nodist_test_qed_eemm_SOURCES = amplitude_qed_eemm.f90 test_qed_eemm_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la amplitude_qed_eemm.f90: $(OMEGA_QED) Makefile $(OMEGA_QED) $(OMEGA_QED_OPTS) -target:module amplitude_qed_eemm \ -scatter "e+ e- -> m+ m-" > $@ test_qed_eemm.o: amplitude_qed_eemm.o test_qed_eemm.o: parameters_QED.o amplitude_qed_eemm.o: parameters_QED.o endif ######################################################################## EXTENDED_COLOR_TESTS = \ $(srcdir)/fc_s.ects \ $(srcdir)/fc_a.ects $(srcdir)/cf_a.ects $(srcdir)/fa_f.ects \ $(srcdir)/ca_c.ects $(srcdir)/af_f.ects $(srcdir)/ac_c.ects \ $(srcdir)/aa_a.ects \ $(srcdir)/fc_fc.ects \ $(srcdir)/aa_s.ects $(srcdir)/as_a.ects $(srcdir)/sa_a.ects TESTS += ects EXTRA_PROGRAMS += ects EXTRA_DIST += ects_driver.sh $(EXTENDED_COLOR_TESTS) # Explicitly state dependence on model files ects.f90: $(OMEGA_QCD) $(OMEGA_SYM) $(OMEGA_SM) ects.f90: ects_driver.sh $(EXTENDED_COLOR_TESTS) @if $(AM_V_P); then :; else echo " ECTS_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/ects_driver.sh \ $(OMEGA_XXX) $(EXTENDED_COLOR_TESTS) > $@ ects_SOURCES = color_test_lib.f90 \ parameters_SM.f90 parameters_QED.f90 parameters_QCD.f90 parameters_SYM.f90 nodist_ects_SOURCES = ects.f90 ects_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la ######################################################################## TESTS += cascade # if there is some debugging output ... # XFAIL_TESTS += cascade CASCADE_TESTS = \ bhabha-s-channel.cascade bhabha-t-channel.cascade bhabha-full.cascade \ ww-onlycc.cascade ww-notgc.cascade \ jjj-notgc.cascade \ vbf-noh.cascade cascade: cascade_driver.sh Makefile $(SED) -e 's|%%cascade_tests%%|$(CASCADE_TESTS)|' \ -e 's|%%srcdir%%|$(srcdir)|' \ -e 's|%%SED%%|$(SED)|' \ -e 's|%%top_builddir%%|$(top_builddir)|' \ -e 's|%%OCAML_NATIVE_EXT%%|$(OCAML_NATIVE_EXT)|' $< >$@ chmod +x $@ EXTRA_DIST += cascade_driver.sh $(CASCADE_TESTS) ######################################################################## TESTS += phase_space PHASE_SPACE_TESTS = eeee.phs qqggg.phs phase_space: phase_space_driver.sh Makefile $(SED) -e 's|%%phase_space_tests%%|$(PHASE_SPACE_TESTS)|' \ -e 's|%%srcdir%%|$(srcdir)|' \ -e 's|%%SED%%|$(SED)|' \ -e 's|%%top_builddir%%|$(top_builddir)|' \ -e 's|%%OCAML_NATIVE_EXT%%|$(OCAML_NATIVE_EXT)|' $< >$@ chmod +x $@ EXTRA_DIST += phase_space_driver.sh $(PHASE_SPACE_TESTS) ######################################################################## TESTS += fermi # XFAIL_TESTS += fermi EXTRA_PROGRAMS += fermi EXTRA_DIST += fermi_driver.sh EXTRA_DIST += fermi.list FERMI_SUPPORT_F90 = \ omega_interface.f90 omega_testtools.f90 tao_random_numbers.f90 \ parameters_QED.f90 parameters_QCD.f90 parameters_SYM.f90 \ parameters_SM.f90 parameters_MSSM.f90 parameters_SM_top_anom.f90 FERMI_SUPPORT_O = $(FERMI_SUPPORT_F90:.f90=.o) fermi_lib.o: $(FERMI_SUPPORT_O) FERMI_LIB_F90 = fermi_lib.f90 $(FERMI_SUPPORT_F90) FERMI_LIB_O = $(FERMI_LIB_F90:.f90=.o) run_fermi: fermi ./fermi fermi.f90: fermi_driver.sh $(OMEGA_QED) $(OMEGA_QCD) $(OMEGA_SYM) fermi.f90: $(OMEGA_SM) $(OMEGA_SM_top_anom) fermi.f90: fermi.list @if $(AM_V_P); then :; else echo " FERMI_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/fermi_driver.sh \ $(OMEGA_XXX) $(OMEGA_SPLIT) < $< > $@ fermi_SOURCES = $(FERMI_LIB_F90) nodist_fermi_SOURCES = fermi.f90 fermi_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la fermi.o: $(FERMI_LIB_O) ######################################################################## TESTS += ward EXTRA_PROGRAMS += ward EXTRA_DIST += ward_driver.sh EXTRA_DIST += ward_identities.list WARD_SUPPORT_F90 = \ omega_interface.f90 omega_testtools.f90 tao_random_numbers.f90 \ parameters_QED.f90 parameters_QCD.f90 parameters_SYM.f90 \ parameters_SM.f90 parameters_SM_top_anom.f90 WARD_SUPPORT_O = $(WARD_SUPPORT_F90:.f90=.o) ward_lib.o: $(WARD_SUPPORT_O) WARD_LIB_F90 = ward_lib.f90 $(WARD_SUPPORT_F90) WARD_LIB_O = $(WARD_LIB_F90:.f90=.o) run_ward: ward ./ward ward.f90: ward_driver.sh $(OMEGA_QED) $(OMEGA_QCD) $(OMEGA_SYM) ward.f90: $(OMEGA_SM) $(OMEGA_SM_top_anom) ward.f90: ward_identities.list @if $(AM_V_P); then :; else echo " WARD_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/ward_driver.sh \ $(OMEGA_XXX) $(OMEGA_SPLIT) < $< > $@ ward_SOURCES = $(WARD_LIB_F90) nodist_ward_SOURCES = ward.f90 ward_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la ward.o: $(WARD_LIB_O) ######################################################################## EXTRA_PROGRAMS += ward_long EXTRA_DIST += ward_identities_long.list run_ward_long: ward_long ./ward_long ward_long.f90: ward_driver.sh ward_long.f90: ward_identities_long.list @if $(AM_V_P); then :; else echo " WARD_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/ward_driver.sh \ $(OMEGA_XXX) $(OMEGA_SPLIT) < $< > $@ ward_long_SOURCES = $(WARD_LIB_F90) nodist_ward_long_SOURCES = ward_long.f90 ward_long_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la # ward_long.o: ward_long.f90 # $(FCCOMPILE) -c -o $@ $(FCFLAGS_f90) -O0 $< ward_long.o: $(WARD_LIB_O) ######################################################################## EXTRA_PROGRAMS += ward_fail EXTRA_DIST += ward_identities_fail.list run_ward_fail: ward_fail ./ward_fail ward_fail.f90: ward_driver.sh ward_fail.f90: ward_identities_fail.list @if $(AM_V_P); then :; else echo " WARD_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/ward_driver.sh \ $(OMEGA_XXX) $(OMEGA_SPLIT) < $< > $@ ward_fail_SOURCES = $(WARD_LIB_F90) nodist_ward_fail_SOURCES = ward_fail.f90 ward_fail_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la ward_fail.o: ward_fail.f90 $(FCCOMPILE) -c -o $@ $(FCFLAGS_f90) -O0 $< ward_fail.o: $(WARD_LIB_O) ######################################################################## TESTS += compare_split_function compare_split_module EXTRA_PROGRAMS += compare_split_function compare_split_module EXTRA_DIST += compare_driver.sh EXTRA_DIST += comparisons.list COMPARE_SUPPORT_F90 = $(WARD_SUPPORT_F90) COMPARE_SUPPORT_O = $(WARD_SUPPORT_O) compare_lib.o: $(COMPARE_SUPPORT_O) COMPARE_LIB_F90 = compare_lib.f90 $(COMPARE_SUPPORT_F90) COMPARE_LIB_O = $(COMPARE_LIB_F90:.f90=.o) run_compare: compare_split_function compare_split_module ./compare_split_function ./compare_split_module compare_split_function.f90: comparisons.list @if $(AM_V_P); then :; else echo " COMPARE_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/compare_driver.sh SF \ "$(OMEGA_XXX) -target:single_function" \ "$(OMEGA_XXX) -target:split_function 10" < $< > $@ compare_split_module.f90: comparisons.list @if $(AM_V_P); then :; else echo " COMPARE_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/compare_driver.sh SM \ "$(OMEGA_XXX) -target:single_function" \ "$(OMEGA_XXX) -target:split_module 10" < $< > $@ compare_split_function.f90 compare_split_module.f90: \ compare_driver.sh $(OMEGA_QCD) $(OMEGA_SM) compare_split_function_SOURCES = $(COMPARE_LIB_F90) nodist_compare_split_function_SOURCES = compare_split_function.f90 compare_split_function_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la compare_split_module_SOURCES = $(COMPARE_LIB_F90) nodist_compare_split_module_SOURCES = compare_split_module.f90 compare_split_module_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la compare_split_function.o compare_split_module.o: $(COMPARE_LIB_O) ######################################################################## if OCAML_AVAILABLE TESTS += compare_majorana compare_majorana_legacy compare_majorana_UFO # XFAIL_TESTS += compare_majorana_UFO EXTRA_PROGRAMS += compare_majorana compare_majorana_legacy compare_majorana_UFO EXTRA_DIST += compare_driver_majorana.sh compare_driver_majorana_UFO.sh EXTRA_DIST += comparisons_majorana.list comparisons_majorana_legacy.list \ comparisons_majorana_UFO.list compare_majorana.f90: comparisons_majorana.list @if $(AM_V_P); then :; else echo " COMPARE_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/compare_driver_majorana.sh Maj \ "$(OMEGA_XXX)" "$(OMEGA_XXX_MAJORANA)" < $< > $@ compare_majorana_legacy.f90: comparisons_majorana_legacy.list @if $(AM_V_P); then :; else echo " COMPARE_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/compare_driver_majorana.sh MajL \ "$(OMEGA_XXX)" "$(OMEGA_XXX_MAJORANA_LEGACY)" < $< > $@ compare_majorana_UFO.f90: comparisons_majorana_UFO.list @if $(AM_V_P); then :; else echo " COMPARE_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/compare_driver_majorana_UFO.sh MajU \ "$(OMEGA_UFO)" "$(OMEGA_UFO_MAJORANA)" "$(OMEGA_UFO_PATH)" < $< > $@ compare_majorana.f90 compare_majorana_legacy.f90 compare_majorana_UFO.f90: \ compare_driver_majorana.sh $(OMEGA_UFO) $(OMEGA_UFO_MAJORANA) compare_majorana_SOURCES = $(COMPARE_LIB_F90) nodist_compare_majorana_SOURCES = compare_majorana.f90 compare_majorana_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la compare_majorana_legacy_SOURCES = $(COMPARE_LIB_F90) nodist_compare_majorana_legacy_SOURCES = compare_majorana_legacy.f90 compare_majorana_legacy_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la compare_majorana_UFO_SOURCES = $(COMPARE_LIB_F90) parameters_SM_UFO.f90 nodist_compare_majorana_UFO_SOURCES = compare_majorana_UFO.f90 compare_majorana_UFO_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la compare_majorana.o compare_majorana_legacy.o compare_majorana_UFO.o: $(COMPARE_LIB_O) +compare_majorana_UFO.o: parameters_SM_UFO.o endif ######################################################################## if OCAML_AVAILABLE # At quadruple or extended precision, these tests take waaaaaayyyy too long! if FC_PREC else TESTS += compare_amplitude_UFO # XFAIL_TESTS += compare_amplitude_UFO EXTRA_PROGRAMS += compare_amplitude_UFO EXTRA_DIST += compare_driver_UFO.sh EXTRA_DIST += comparisons_UFO.list compare_amplitude_UFO_SOURCES = \ parameters_SM_from_UFO.f90 compare_lib.f90 \ omega_interface.f90 omega_testtools.f90 tao_random_numbers.f90 compare_amplitude_UFO.f90: comparisons_UFO.list compare_driver_UFO.sh $(OMEGA_UFO) @if $(AM_V_P); then :; else echo " COMPARE_DRIVER_UFO"; fi $(AM_V_at)$(SHELL) $(srcdir)/compare_driver_UFO.sh UFO \ "$(OMEGA_XXX) -model:constant_width" \ "$(OMEGA_UFO) -model:UFO_dir $(top_srcdir)/omega/tests/UFO/%%%/ -model:exec" \ < $< > $@ # -model:long_flavors nodist_compare_amplitude_UFO_SOURCES = \ compare_amplitude_UFO.f90 parameters_SM_UFO.f90 compare_amplitude_UFO_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la parameters_SM_from_UFO.o: parameters_SM_UFO.o compare_amplitude_UFO.o: parameters_SM_UFO.o parameters_SM_from_UFO.o compare_amplitude_UFO.o: $(COMPARE_LIB_O) endif parameters_SM_UFO.f90: $(OMEGA_UFO) $(OMEGA_UFO) \ -model:UFO_dir $(OMEGA_UFO_PATH)/SM/ -model:exec \ -target:parameter_module parameters_sm_ufo -params > $@ endif ######################################################################## if OCAML_AVAILABLE # At quadruple or extended precision, these tests take waaaaaayyyy too long! if FC_PREC else TESTS += fermi_UFO # XFAIL_TESTS += fermi_UFO # We need more work on the parameters to pass the tests # at quadruple or extended precision. if FC_PREC XFAIL_TESTS += fermi_UFO endif EXTRA_PROGRAMS += fermi_UFO EXTRA_DIST += fermi_driver_UFO.sh EXTRA_DIST += fermi_UFO.list FERMI_UFO_SUPPORT_F90 = \ omega_interface.f90 omega_testtools.f90 tao_random_numbers.f90 FERMI_UFO_SUPPORT_O = $(FERMI_UFO_SUPPORT_F90:.f90=.o) fermi_UFO_lib.o: $(FERMI_SUPPORT_O) FERMI_UFO_LIB_F90 = fermi_lib.f90 $(FERMI_UFO_SUPPORT_F90) FERMI_UFO_LIB_O = $(FERMI_UFO_LIB_F90:.f90=.o) run_fermi_UFO: fermi_UFO ./fermi_UFO fermi_UFO.f90: fermi_UFO.list fermi_driver_UFO.sh $(OMEGA_UFO) @if $(AM_V_P); then :; else echo " FERMI_UFO_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/fermi_driver_UFO.sh \ $(OMEGA_UFO) $(OMEGA_UFO_MAJORANA) $(OMEGA_UFO_PATH) \ $(OMEGA_SPLIT) < $< > $@ fermi_UFO_SOURCES = $(FERMI_UFO_LIB_F90) nodist_fermi_UFO_SOURCES = fermi_UFO.f90 parameters_SM_UFO.f90 fermi_UFO_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la -fermi_UFO.o: $(FERMI_UFO_LIB_O) +fermi_UFO.o: $(FERMI_UFO_LIB_O) parameters_SM_UFO.o endif endif ######################################################################## if OCAML_AVAILABLE # At quadruple or extended precision, these tests take waaaaaayyyy too long! if FC_PREC else TESTS += ward_UFO # We need more work on the parameters to pass the tests # at quadruple or extended precision. if FC_PREC XFAIL_TESTS += ward_UFO endif EXTRA_PROGRAMS += ward_UFO EXTRA_DIST += ward_driver_UFO.sh EXTRA_DIST += ward_identities_UFO.list WARD_UFO_SUPPORT_F90 = \ omega_interface.f90 omega_testtools.f90 tao_random_numbers.f90 WARD_UFO_SUPPORT_O = $(WARD_UFO_SUPPORT_F90:.f90=.o) ward_UFO_lib.o: $(WARD_SUPPORT_O) WARD_UFO_LIB_F90 = ward_lib.f90 $(WARD_UFO_SUPPORT_F90) WARD_UFO_LIB_O = $(WARD_UFO_LIB_F90:.f90=.o) run_ward_UFO: ward_UFO ./ward_UFO ward_UFO.f90: ward_identities_UFO.list ward_driver_UFO.sh $(OMEGA_UFO) @if $(AM_V_P); then :; else echo " WARD_UFO_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/ward_driver_UFO.sh \ $(OMEGA_UFO) -model:UFO_dir $(top_srcdir)/omega/tests/UFO/SM/ \ $(OMEGA_SPLIT) < $< > $@ ward_UFO_SOURCES = $(WARD_UFO_LIB_F90) nodist_ward_UFO_SOURCES = ward_UFO.f90 parameters_SM_UFO.f90 ward_UFO_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la -ward_UFO.o: $(WARD_UFO_LIB_O) +ward_UFO.o: $(WARD_UFO_LIB_O) parameters_SM_UFO.o endif endif ######################################################################## TESTS += compare_amplitude_VM EXTRA_PROGRAMS += compare_amplitude_VM EXTRA_DIST += compare_driver_VM.sh compare_driver_VM_wrappers.sh EXTRA_DIST += comparisons_VM.list compare_amplitude_VM.f90: comparisons_VM.list comparisons_VM.wrappers.o @if $(AM_V_P); then :; else echo " COMPARE_DRIVER_VM"; fi $(AM_V_at)$(SHELL) $(srcdir)/compare_driver_VM.sh \ "$(OMEGA_XXX) " "$(OMEGA_XXX_VM) " "$(OMEGA_XXX_VM_PARAMS_OPTS)" < $< > $@ comparisons_VM.wrappers.f90: comparisons_VM.list @if $(AM_V_P); then :; else echo " COMPARE_DRIVER_VM_WRAPPERS"; fi $(AM_V_at)$(SHELL) $(srcdir)/compare_driver_VM_wrappers.sh \ "$(OMEGA_XXX) " "$(OMEGA_XXX_VM) " "$(OMEGA_XXX_VM_PARAMS_OPTS)" < $< > $@ # Explicitly state dependence on model files compare_amplitude_VM.f90: compare_driver_VM.sh \ $(OMEGA_QED) $(OMEGA_QED_VM) \ $(OMEGA_QCD) $(OMEGA_QCD_VM) \ $(OMEGA_SM) $(OMEGA_SM_VM) \ $(OMEGA_SM_CKM) $(OMEGA_SM_CKM_VM) \ $(OMEGA_SM_Higgs) $(OMEGA_SM_Higgs_VM) \ $(OMEGA_THDM) $(OMEGA_THDM_VM) \ $(OMEGA_THDM_CKM) $(OMEGA_THDM_CKM_VM) \ $(OMEGA_HSExt) $(OMEGA_HSExt_VM) \ $(OMEGA_Zprime) $(OMEGA_Zprime_VM) COMPARE_EXTRA_MODELS = parameters_SM_CKM.f90 parameters_SM_Higgs.f90 \ parameters_THDM.f90 parameters_THDM_CKM.f90 parameters_HSExt.f90 \ parameters_Zprime.f90 compare_amplitude_VM_SOURCES = $(COMPARE_LIB_F90) $(COMPARE_EXTRA_MODELS) nodist_compare_amplitude_VM_SOURCES = compare_amplitude_VM.f90 comparisons_VM.wrappers.f90 compare_amplitude_VM_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la compare_amplitude_VM.o: $(COMPARE_LIB_O) ######################################################################## if FC_USE_OPENMP TESTS += test_openmp EXTRA_PROGRAMS += test_openmp TESTOPENMP_SUPPORT_F90 = $(WARD_SUPPORT_F90) TESTOPENMP_SUPPORT_O = $(WARD_SUPPORT_O) test_openmp_SOURCES = test_openmp.f90 $(TESTOPENMP_SUPPORT_F90) nodist_test_openmp_SOURCES = amplitude_openmp.f90 test_openmp_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la amplitude_openmp.f90: $(OMEGA_QCD) Makefile $(OMEGA_QCD) $(OMEGA_QCD_OPTS) -target:openmp \ -target:module amplitude_openmp -scatter "gl gl -> gl gl gl" > $@ test_openmp.o: amplitude_openmp.o test_openmp.o: $(TESTOPENMP_SUPPORT_O) amplitude_openmp.o: parameters_QCD.o endif ######################################################################## EXTRA_PROGRAMS += benchmark_VM_vs_Fortran EXTRA_DIST += benchmark_VM_vs_Fortran_driver.sh BENCHMARK_LIB_F90 = benchmark_lib.f90 $(WARD_SUPPORT_F90) BENCHMARK_LIB_O = $(BENCHMARK_LIB_F90:.f90=.o) benchmark_VM_vs_Fortran.f90: benchmark_processes.list benchmark_processes.wrappers.o @if $(AM_V_P); then :; else echo " BENCHMARK_VM_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/benchmark_VM_vs_Fortran_driver.sh \ "$(OMEGA_XXX) " "$(OMEGA_XXX_VM) " "$(OMEGA_XXX_VM_PARAMS_OPTS)" < $< > $@ benchmark_processes.wrappers.f90: benchmark_processes.list @if $(AM_V_P); then :; else echo " BENCHMARK_DRIVER_WRAPPERS"; fi $(AM_V_at)$(SHELL) $(srcdir)/benchmark_driver_wrappers.sh \ "$(OMEGA_XXX) " "$(OMEGA_XXX_VM) " "$(OMEGA_XXX_VM_PARAMS_OPTS)" < $< > $@ # Explicitly state dependence on model files benchmark_VM_vs_Fortran.f90: benchmark_VM_vs_Fortran_driver.sh \ $(OMEGA_QED) $(OMEGA_QED_VM) \ $(OMEGA_QCD) $(OMEGA_QCD_VM) \ $(OMEGA_SM) $(OMEGA_SM_VM) benchmark_VM_vs_Fortran_SOURCES = $(BENCHMARK_LIB_F90) nodist_benchmark_VM_vs_Fortran_SOURCES = benchmark_VM_vs_Fortran.f90 benchmark_processes.wrappers.f90 benchmark_VM_vs_Fortran_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la benchmark_VM_vs_Fortran.o: $(BENCHMARK_LIB_O) ######################################################################## if FC_USE_OPENMP EXTRA_PROGRAMS += benchmark_amp_parallel benchmark_amp_parallel.f90: benchmark_processes.list benchmark_processes.wrappers.o @if $(AM_V_P); then :; else echo " BENCHMARK_PARALLEL_DRIVER"; fi $(AM_V_at)$(SHELL) $(srcdir)/benchmark_amp_parallel_driver.sh \ "$(OMEGA_XXX) " "$(OMEGA_XXX_VM) " "$(OMEGA_XXX_VM_PARAMS_OPTS)" < $< > $@ # Explicitly state dependence on model files benchmark_amp_parallel.f90: benchmark_amp_parallel_driver.sh \ $(OMEGA_QED) $(OMEGA_QED_VM) \ $(OMEGA_QCD) $(OMEGA_QCD_VM) \ $(OMEGA_SM) $(OMEGA_SM_VM) benchmark_amp_parallel_SOURCES = $(BENCHMARK_LIB_F90) nodist_benchmark_amp_parallel_SOURCES = benchmark_amp_parallel.f90 benchmark_processes.wrappers.f90 benchmark_amp_parallel_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la benchmark_amp_parallel.o: $(BENCHMARK_LIB_O) endif ######################################################################## EXTRA_PROGRAMS += benchmark run_benchmark: benchmark ./benchmark BENCHMARK_PROCESS = -scatter "gl gl -> gl gl gl" BENCHMARK_SPLIT_SIZE = 10 benchmark_SOURCES = benchmark.f90 parameters_QCD.f90 nodist_benchmark_SOURCES = \ amplitude_benchmark_v1.f90 amplitude_benchmark_v2.f90 \ amplitude_benchmark_v3.f90 # amplitude_benchmark_v4.f90 benchmark_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la amplitude_benchmark_v1.f90: $(OMEGA_QCD) Makefile $(OMEGA_QCD) $(OMEGA_QCD_OPTS) -target:module amplitude_benchmark_v1 \ $(BENCHMARK_PROCESS) -target:single_function > $@ amplitude_benchmark_v2.f90: $(OMEGA_QCD) Makefile $(OMEGA_QCD) $(OMEGA_QCD_OPTS) -target:module amplitude_benchmark_v2 \ $(BENCHMARK_PROCESS) -target:split_function $(BENCHMARK_SPLIT_SIZE) > $@ amplitude_benchmark_v3.f90: $(OMEGA_QCD) Makefile $(OMEGA_QCD) $(OMEGA_QCD_OPTS) -target:module amplitude_benchmark_v3 \ $(BENCHMARK_PROCESS) -target:split_module $(BENCHMARK_SPLIT_SIZE) > $@ amplitude_benchmark_v4.f90: $(OMEGA_QCD) Makefile $(OMEGA_QCD) $(OMEGA_QCD_OPTS) -target:module amplitude_benchmark_v4 \ $(BENCHMARK_PROCESS) -target:split_file $(BENCHMARK_SPLIT_SIZE) > $@ benchmark.o: \ amplitude_benchmark_v1.o amplitude_benchmark_v2.o \ amplitude_benchmark_v3.o # amplitude_benchmark_v4.o benchmark.o: parameters_QCD.o amplitude_benchmark_v1.o amplitude_benchmark_v2.o \ amplitude_benchmark_v3.o amplitude_benchmark_v4.o: parameters_QCD.o ######################################################################## EXTRA_PROGRAMS += benchmark_UFO_SM run_benchmark_UFO_SM: benchmark_UFO_SM ./benchmark_UFO_SM # NB: This IS portable ... UFO_SM = $(OMEGA_UFO_PATH)/SM/ BENCHMARK_UFO_SM_PROCESS = -scatter "e+ e- -> W+ W- Z Z" benchmark_UFO_SM_SOURCES = \ benchmark_UFO_SM.f90 parameters_SM_from_UFO.f90 nodist_benchmark_UFO_SM_SOURCES = \ amplitude_benchmark_UFO_SM.f90 \ amplitude_benchmark_UFO_SM_classic.f90 \ parameters_SM_UFO.f90 benchmark_UFO_SM_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la amplitude_benchmark_UFO_SM_classic.f90: $(OMEGA_SM) Makefile $(OMEGA_SM) -target:module amplitude_benchmark_UFO_SM_classic \ -target:parameter_module parameters_SM_from_UFO \ $(BENCHMARK_UFO_SM_PROCESS) > $@ amplitude_benchmark_UFO_SM.f90: $(OMEGA_UFO) Makefile $(OMEGA_UFO) -model:UFO_dir $(UFO_SM) -model:exec \ -target:module amplitude_benchmark_UFO_SM \ -target:parameter_module parameters_SM_UFO \ $(BENCHMARK_UFO_SM_PROCESS) > $@ benchmark_UFO_SM.o: \ amplitude_benchmark_UFO_SM.o amplitude_benchmark_UFO_SM_classic.o benchmark_UFO_SM.o: parameters_SM_UFO.o parameters_SM_from_UFO.o amplitude_benchmark_UFO_SM_classic.o: parameters_SM_from_UFO.o amplitude_benchmark_UFO_SM.o: parameters_SM_UFO.o ######################################################################## EXTRA_PROGRAMS += benchmark_UFO_SMEFT run_benchmark_UFO_SMEFT: benchmark_UFO_SMEFT ./benchmark_UFO_SMEFT # NB: This is NOT portable ... UFO_SMEFT = /home/ohl/physics/SMEFT_mW_UFO/ BENCHMARK_UFO_SMEFT_PROCESS = -scatter "e+ e- -> W+ W- Z" benchmark_UFO_SMEFT_SOURCES = benchmark_UFO_SMEFT.f90 nodist_benchmark_UFO_SMEFT_SOURCES = \ amplitude_benchmark_UFO_SMEFT.f90 \ amplitude_benchmark_UFO_SMEFT_opt.f90 \ parameters_UFO_SMEFT.f90 benchmark_UFO_SMEFT_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la amplitude_benchmark_UFO_SMEFT.f90: $(OMEGA_UFO) Makefile $(OMEGA_UFO) -model:UFO_dir $(UFO_SMEFT) -model:exec \ -target:module amplitude_benchmark_UFO_SMEFT \ -target:parameter_module parameters_UFO_SMEFT \ $(BENCHMARK_UFO_SMEFT_PROCESS) | $(SED) 's/g == 0/.false./' > $@ amplitude_benchmark_UFO_SMEFT_opt.f90: $(OMEGA_UFO) Makefile $(OMEGA_UFO) -model:UFO_dir $(UFO_SMEFT) -model:exec \ -target:module amplitude_benchmark_UFO_SMEFT_opt \ -target:parameter_module parameters_UFO_SMEFT \ $(BENCHMARK_UFO_SMEFT_PROCESS) > $@ benchmark_UFO_SMEFT.o: \ amplitude_benchmark_UFO_SMEFT.o amplitude_benchmark_UFO_SMEFT_opt.o benchmark_UFO_SMEFT.o: parameters_UFO_SMEFT.o amplitude_benchmark_UFO_SMEFT.o amplitude_benchmark_UFO_SMEFT_opt.o: \ parameters_UFO_SMEFT.o parameters_UFO_SMEFT.f90: $(OMEGA_UFO) $(OMEGA_UFO) -model:UFO_dir $(UFO_SMEFT) -model:exec \ -target:parameter_module parameters_UFO_SMEFT -params > $@ ######################################################################## if OCAML_AVAILABLE TESTS += vertex_unit EXTRA_PROGRAMS += vertex_unit vertex_unit_SOURCES = vertex_unit.ml vertex_unit: $(OMEGA_CORE) vertex_unit.cmx @if $(AM_V_P); then :; else echo " OCAMLOPT " $@; fi $(AM_V_at)$(OCAMLOPT) $(OCAMLFLAGS) $(OCAMLOPTFLAGS) -o vertex_unit \ unix.cmxa $(OMEGA_CORE) $(OMEGA_MODELS) vertex_unit.cmx vertex_unit.cmx: vertex_unit.ml vertex_unit.cmx: $(OMEGA_CORE) $(OMEGA_MODELS) endif ######################################################################## if OCAML_AVAILABLE TESTS += ufo_unit EXTRA_PROGRAMS += ufo_unit ufo_unit_SOURCES = ufo_unit.ml ufo_unit: $(OMEGA_CORE) ufo_unit.cmx @if $(AM_V_P); then :; else echo " OCAMLOPT " $@; fi $(AM_V_at)$(OCAMLOPT) $(OCAMLFLAGS) $(OCAMLOPTFLAGS) -o ufo_unit \ unix.cmxa $(OMEGA_CORE) $(OMEGA_MODELS) ufo_unit.cmx ufo_unit.cmx: ufo_unit.ml ufo_unit.cmx: $(OMEGA_CORE) $(OMEGA_MODELS) endif ######################################################################## if OCAML_AVAILABLE TESTS += keystones_omegalib keystones_UFO TESTS += keystones_omegalib_bispinors keystones_UFO_bispinors # XFAIL_TESTS += keystones_UFO # XFAIL_TESTS += keystones_UFO_bispinors EXTRA_PROGRAMS += keystones_omegalib keystones_UFO EXTRA_PROGRAMS += keystones_omegalib_bispinors keystones_UFO_bispinors keystones_omegalib_SOURCES = omega_testtools.f90 keystones_tools.f90 nodist_keystones_omegalib_SOURCES = keystones_omegalib.f90 keystones_omegalib_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la keystones_UFO_SOURCES = omega_testtools.f90 keystones_tools.f90 nodist_keystones_UFO_SOURCES = keystones_UFO.f90 keystones_UFO_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la keystones_omegalib_bispinors_SOURCES = omega_testtools.f90 keystones_tools.f90 nodist_keystones_omegalib_bispinors_SOURCES = keystones_omegalib_bispinors.f90 keystones_omegalib_bispinors_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la keystones_UFO_bispinors_SOURCES = omega_testtools.f90 keystones_tools.f90 nodist_keystones_UFO_bispinors_SOURCES = keystones_UFO_bispinors.f90 keystones_UFO_bispinors_LDADD = $(KINDS) $(top_builddir)/omega/src/libomega_core.la EXTRA_PROGRAMS += keystones_omegalib_generate keystones_UFO_generate EXTRA_PROGRAMS += keystones_omegalib_bispinors_generate keystones_UFO_bispinors_generate keystones_omegalib_generate_SOURCES = \ keystones.ml keystones.mli keystones_omegalib_generate.ml keystones_UFO_generate_SOURCES = \ keystones.ml keystones.mli keystones_UFO_generate.ml keystones_omegalib_bispinors_generate_SOURCES = \ keystones.ml keystones.mli keystones_omegalib_bispinors_generate.ml keystones_UFO_bispinors_generate_SOURCES = \ keystones.ml keystones.mli keystones_UFO_bispinors_generate.ml keystones_omegalib.f90: keystones_omegalib_generate ./keystones_omegalib_generate -cat > $@ keystones_UFO.f90: keystones_UFO_generate ./keystones_UFO_generate -cat > $@ keystones_omegalib_bispinors.f90: keystones_omegalib_bispinors_generate ./keystones_omegalib_bispinors_generate -cat > $@ keystones_UFO_bispinors.f90: keystones_UFO_bispinors_generate ./keystones_UFO_bispinors_generate -cat > $@ keystones_omegalib_generate: $(OMEGA_CORE) keystones_omegalib_generate.cmx @if $(AM_V_P); then :; else echo " OCAMLOPT " $@; fi $(AM_V_at)$(OCAMLOPT) $(OCAMLFLAGS) $(OCAMLOPTFLAGS) \ -o keystones_omegalib_generate \ unix.cmxa $(OMEGA_CORE) $(OMEGA_MODELS) \ keystones.cmx keystones_omegalib_generate.cmx keystones_UFO_generate: $(OMEGA_CORE) keystones_UFO_generate.cmx @if $(AM_V_P); then :; else echo " OCAMLOPT " $@; fi $(AM_V_at)$(OCAMLOPT) $(OCAMLFLAGS) $(OCAMLOPTFLAGS) \ -o keystones_UFO_generate \ unix.cmxa $(OMEGA_CORE) $(OMEGA_MODELS) \ keystones.cmx keystones_UFO_generate.cmx keystones_omegalib_bispinors_generate: $(OMEGA_CORE) keystones_omegalib_bispinors_generate.cmx @if $(AM_V_P); then :; else echo " OCAMLOPT " $@; fi $(AM_V_at)$(OCAMLOPT) $(OCAMLFLAGS) $(OCAMLOPTFLAGS) \ -o keystones_omegalib_bispinors_generate \ unix.cmxa $(OMEGA_CORE) $(OMEGA_MODELS) \ keystones.cmx keystones_omegalib_bispinors_generate.cmx keystones_UFO_bispinors_generate: $(OMEGA_CORE) keystones_UFO_bispinors_generate.cmx @if $(AM_V_P); then :; else echo " OCAMLOPT " $@; fi $(AM_V_at)$(OCAMLOPT) $(OCAMLFLAGS) $(OCAMLOPTFLAGS) \ -o keystones_UFO_bispinors_generate \ unix.cmxa $(OMEGA_CORE) $(OMEGA_MODELS) \ keystones.cmx keystones_UFO_bispinors_generate.cmx keystones_omegalib_generate.cmx: \ keystones.cmi keystones.cmx keystones_omegalib_generate.ml keystones_omegalib_generate.cmx: $(OMEGA_CORE) $(OMEGA_MODELS) keystones_UFO_generate.cmx: \ keystones.cmi keystones.cmx keystones_UFO_generate.ml keystones_UFO_generate.cmx: $(OMEGA_CORE) $(OMEGA_MODELS) keystones_omegalib_bispinors_generate.cmx: \ keystones.cmi keystones.cmx keystones_omegalib_bispinors_generate.ml keystones_omegalib_bispinors_generate.cmx: $(OMEGA_CORE) $(OMEGA_MODELS) keystones_UFO_bispinors_generate.cmx: \ keystones.cmi keystones.cmx keystones_UFO_bispinors_generate.ml keystones_UFO_bispinors_generate.cmx: $(OMEGA_CORE) $(OMEGA_MODELS) keystones.cmx: keystones.ml keystones.cmi keystones.cmx: $(OMEGA_CORE) $(OMEGA_MODELS) keystones.cmi: keystones.mli $(OMEGA_CORE) endif ######################################################################## if RECOLA_AVAILABLE TESTS += compare_amplitude_recola # We need more work on the parameters to pass the tests # at quadruple or extended precision if FC_PREC XFAIL_TESTS += compare_amplitude_recola endif EXTRA_PROGRAMS += compare_amplitude_recola AM_FCFLAGS += $(RECOLA_INCLUDES) compare_amplitude_recola_SOURCES = \ parameters_SM_Higgs_recola.f90 \ omega_interface.f90 compare_lib.f90 compare_lib_recola.f90 \ omega_testtools.f90 tao_random_numbers.f90 nodist_compare_amplitude_recola_SOURCES = compare_amplitude_recola.f90 compare_amplitude_recola.f90: comparisons_recola.list compare_driver_recola.sh @if $(AM_V_P); then :; else echo " COMPARE_DRIVER_RECOLA"; fi $(AM_V_at)$(SHELL) $(srcdir)/compare_driver_recola.sh \ "$(OMEGA_XXX) -model:constant_width" < $< > $@ compare_amplitude_recola.o: \ omega_testtools.f90 compare_lib.o compare_lib_recola.o \ tao_random_numbers.o \ parameters_SM_Higgs_recola.o compare_lib_recola.o: \ omega_testtools.f90 compare_lib.o tao_random_numbers.o \ parameters_SM_Higgs_recola.o compare_amplitude_recola_LDADD = \ $(LDFLAGS_RECOLA) \ $(KINDS) $(top_builddir)/omega/src/libomega_core.la run_compare_recola: compare_amplitude_recola ./compare_amplitude_recola endif ######################################################################## installcheck-local: PATH=$(DESTDIR)$(bindir):$$PATH; export PATH; \ LD_LIBRARY_PATH=$(DESTDIR)$(libdir):$(DESTDIR)$(pkglibdir):$$LD_LIBRARY_PATH; \ export LD_LIBRARY_PATH; \ omega_QED.opt $(OMEGA_QED_OPTS) -scatter "e+ e- -> m+ m-" \ -target:module amplitude_qed_eemm > amplitude_qed_eemm.f90; \ $(FC) $(AM_FCFLAGS) $(FCFLAGS) -I$(pkgincludedir) \ -L$(DESTDIR)$(libdir) -L$(DESTDIR)$(pkglibdir) \ $(srcdir)/parameters_QED.f90 amplitude_qed_eemm.f90 \ $(srcdir)/test_qed_eemm.f90 -lomega_core; \ ./a.out ######################################################################## ### Remove DWARF debug information on MAC OS X clean-macosx: -rm -rf a.out.dSYM -rm -rf compare_amplitude_UFO.dSYM -rm -rf compare_amplitude_VM.dSYM -rm -rf compare_split_function.dSYM -rm -rf compare_split_module.dSYM -rm -rf ects.dSYM -rm -rf test_omega95.dSYM -rm -rf test_omega95_bispinors.dSYM -rm -rf test_qed_eemm.dSYM -rm -rf ward.dSYM .PHONY: clean-macosx clean-local: clean-macosx rm -f a.out gmon.out *.$(FCMOD) \ *.o *.cmi *.cmo *.cmx amplitude_*.f90 \ $(EXTRA_PROGRAMS) ects.f90 ward.f90 ward_UFO.f90 \ fermi.f90 fermi_UFO.f90 compare_*.f90 \ parameters_SM_UFO.f90 keystones_omegalib.f90 keystones_UFO.f90 \ keystones_UFO_bispinors.f90 keystones_omegalib_bispinors.f90 \ omega_testtools.f90 test_omega95*.f90 benchmark*.f90 \ parameters_UFO_SMEFT.f90 \ *.hbc *wrappers.f90 cascade phase_space \ output.rcl recola.log rm -fr output_cll if FC_SUBMODULES -rm -f *.smod endif ######################################################################## ## The End. ######################################################################## Index: trunk/omega/src/UFO.ml =================================================================== --- trunk/omega/src/UFO.ml (revision 8494) +++ trunk/omega/src/UFO.ml (revision 8495) @@ -1,2915 +1,2909 @@ (* UFO.ml -- Copyright (C) 1999-2021 by Wolfgang Kilian Thorsten Ohl Juergen Reuter with contributions from Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) (* Unfortunately, \texttt{ocamlweb} will not typeset all multi character operators nicely. E.\,g.~\verb+f @< g+ comes out as [f @< g]. *) let (@@) f g x = f (g x) let (@@@) f g x y = f (g x y) module SMap = Map.Make (struct type t = string let compare = compare end) module SSet = Sets.String module CMap = Map.Make (struct type t = string let compare = ThoString.compare_caseless end) module CSet = Sets.String_Caseless let error_in_string text start_pos end_pos = let i = start_pos.Lexing.pos_cnum and j = end_pos.Lexing.pos_cnum in String.sub text i (j - i) let error_in_file name start_pos end_pos = Printf.sprintf "%s:%d.%d-%d.%d" name start_pos.Lexing.pos_lnum (start_pos.Lexing.pos_cnum - start_pos.Lexing.pos_bol) end_pos.Lexing.pos_lnum (end_pos.Lexing.pos_cnum - end_pos.Lexing.pos_bol) let parse_string text = try UFO_parser.file UFO_lexer.token (UFO_lexer.init_position "" (Lexing.from_string text)) with | UFO_tools.Lexical_Error (msg, start_pos, end_pos) -> invalid_arg (Printf.sprintf "lexical error (%s) at: `%s'" msg (error_in_string text start_pos end_pos)) | UFO_syntax.Syntax_Error (msg, start_pos, end_pos) -> invalid_arg (Printf.sprintf "syntax error (%s) at: `%s'" msg (error_in_string text start_pos end_pos)) | Parsing.Parse_error -> invalid_arg ("parse error: " ^ text) let parse_file name = let ic = open_in name in let result = begin try UFO_parser.file UFO_lexer.token (UFO_lexer.init_position name (Lexing.from_channel ic)) with | UFO_tools.Lexical_Error (msg, start_pos, end_pos) -> begin close_in ic; invalid_arg (Printf.sprintf "%s: lexical error (%s)" (error_in_file name start_pos end_pos) msg) end | UFO_syntax.Syntax_Error (msg, start_pos, end_pos) -> begin close_in ic; invalid_arg (Printf.sprintf "%s: syntax error (%s)" (error_in_file name start_pos end_pos) msg) end | Parsing.Parse_error -> begin close_in ic; invalid_arg ("parse error: " ^ name) end end in close_in ic; result (* These are the contents of the Python files after lexical analysis as context-free variable declarations, before any semantic interpretation. *) module type Files = sig type t = private { particles : UFO_syntax.t; couplings : UFO_syntax.t; coupling_orders : UFO_syntax.t; vertices : UFO_syntax.t; lorentz : UFO_syntax.t; parameters : UFO_syntax.t; propagators : UFO_syntax.t; decays : UFO_syntax.t } val parse_directory : string -> t end module Files : Files = struct type t = { particles : UFO_syntax.t; couplings : UFO_syntax.t; coupling_orders : UFO_syntax.t; vertices : UFO_syntax.t; lorentz : UFO_syntax.t; parameters : UFO_syntax.t; propagators : UFO_syntax.t; decays : UFO_syntax.t } let parse_directory dir = let parse stem = parse_file (Filename.concat dir (stem ^ ".py")) in { particles = parse "particles"; couplings = parse "couplings"; coupling_orders = (try parse "coupling_orders" with _ -> []); vertices = parse "vertices"; lorentz = parse "lorentz"; parameters = parse "parameters"; propagators = parse "propagators"; (* [(try parse "propagators" with _ -> []);] *) decays = (try parse "decays" with _ -> []) } end let dump_file pfx f = List.iter (fun s -> print_endline (pfx ^ ": " ^ s)) (UFO_syntax.to_strings f) type charge = | Q_Integer of int | Q_Fraction of int * int let charge_to_string = function | Q_Integer i -> Printf.sprintf "%d" i | Q_Fraction (n, d) -> Printf.sprintf "%d/%d" n d module S = UFO_syntax let find_attrib name attribs = try (List.find (fun a -> name = a.S.a_name) attribs).S.a_value with | Not_found -> failwith ("UFO.find_attrib: \"" ^ name ^ "\" not found") let find_attrib name attribs = (List.find (fun a -> name = a.S.a_name) attribs).S.a_value let name_to_string ?strip name = let stripped = begin match strip, List.rev name with | Some pfx, head :: tail -> if pfx = head then tail else failwith ("UFO.name_to_string: expected prefix '" ^ pfx ^ "', got '" ^ head ^ "'") | _, name -> name end in String.concat "." stripped let name_attrib ?strip name attribs = match find_attrib name attribs with | S.Name n -> name_to_string ?strip n | _ -> invalid_arg ("UFO.name_attrib: " ^ name) let integer_attrib name attribs = match find_attrib name attribs with | S.Integer i -> i | _ -> invalid_arg ("UFO.integer_attrib: " ^ name) let charge_attrib name attribs = match find_attrib name attribs with | S.Integer i -> Q_Integer i | S.Fraction (n, d) -> Q_Fraction (n, d) | _ -> invalid_arg ("UFO.charge_attrib: " ^ name) let string_attrib name attribs = match find_attrib name attribs with | S.String s -> s | _ -> invalid_arg ("UFO.string_attrib: " ^ name) let string_expr_attrib name attribs = match find_attrib name attribs with | S.Name n -> [S.Macro n] | S.String s -> [S.Literal s] | S.String_Expr e -> e | _ -> invalid_arg ("UFO.string_expr_attrib: " ^ name) let boolean_attrib name attribs = try match ThoString.lowercase (name_attrib name attribs) with | "true" -> true | "false" -> false | _ -> invalid_arg ("UFO.boolean_attrib: " ^ name) with | Not_found -> false type value = | Integer of int | Fraction of int * int | Float of float | Expr of UFOx.Expr.t | Name of string list let map_expr f default = function | Integer _ | Fraction (_, _) | Float _ | Name _ -> default | Expr e -> f e let variables = map_expr UFOx.Expr.variables CSet.empty let functions = map_expr UFOx.Expr.functions CSet.empty let add_to_set_in_map key element map = let set = try CMap.find key map with Not_found -> CSet.empty in CMap.add key (CSet.add element set) map (* Add all variables in [value] to the [map] from variables to the names in which they appear, indicating that [name] depends on these variables. *) let dependency name value map = CSet.fold (fun variable acc -> add_to_set_in_map variable name acc) (variables value) map let dependencies name_value_list = List.fold_left (fun acc (name, value) -> dependency name value acc) CMap.empty name_value_list let dependency_to_string (variable, appearences) = Printf.sprintf "%s -> {%s}" variable (String.concat ", " (CSet.elements appearences)) let dependencies_to_strings map = List.map dependency_to_string (CMap.bindings map) let expr_to_string = UFOx.Value.to_string @@ UFOx.Value.of_expr let value_to_string = function | Integer i -> Printf.sprintf "%d" i | Fraction (n, d) -> Printf.sprintf "%d/%d" n d | Float x -> string_of_float x | Expr e -> "'" ^ expr_to_string e ^ "'" | Name n -> name_to_string n let value_to_expr substitutions = function | Integer i -> Printf.sprintf "%d" i | Fraction (n, d) -> Printf.sprintf "%d/%d" n d | Float x -> string_of_float x | Expr e -> expr_to_string (substitutions e) | Name n -> name_to_string n let value_to_coupling substitutions atom = function | Integer i -> Coupling.Integer i | Fraction (n, d) -> Coupling.Quot (Coupling.Integer n, Coupling.Integer d) | Float x -> Coupling.Float x | Expr e -> UFOx.Value.to_coupling atom (UFOx.Value.of_expr (substitutions e)) | Name n -> failwith "UFO.value_to_coupling: Name not supported yet!" let value_to_numeric = function | Integer i -> Printf.sprintf "%d" i | Fraction (n, d) -> Printf.sprintf "%g" (float n /. float d) | Float x -> Printf.sprintf "%g" x | Expr e -> invalid_arg ("UFO.value_to_numeric: expr = " ^ (expr_to_string e)) | Name n -> invalid_arg ("UFO.value_to_numeric: name = " ^ name_to_string n) let value_to_float = function | Integer i -> float i | Fraction (n, d) -> float n /. float d | Float x -> x | Expr e -> invalid_arg ("UFO.value_to_float: string = " ^ (expr_to_string e)) | Name n -> invalid_arg ("UFO.value_to_float: name = " ^ name_to_string n) let value_attrib name attribs = match find_attrib name attribs with | S.Integer i -> Integer i | S.Fraction (n, d) -> Fraction (n, d) | S.Float x -> Float x | S.String s -> Expr (UFOx.Expr.of_string s) | S.Name n -> Name n | _ -> invalid_arg ("UFO.value_attrib: " ^ name) let string_list_attrib name attribs = match find_attrib name attribs with | S.String_List l -> l | _ -> invalid_arg ("UFO.string_list_attrib: " ^ name) let name_list_attrib ~strip name attribs = match find_attrib name attribs with | S.Name_List l -> List.map (name_to_string ~strip) l | _ -> invalid_arg ("UFO.name_list_attrib: " ^ name) let integer_list_attrib name attribs = match find_attrib name attribs with | S.Integer_List l -> l | _ -> invalid_arg ("UFO.integer_list_attrib: " ^ name) let order_dictionary_attrib name attribs = match find_attrib name attribs with | S.Order_Dictionary d -> d | _ -> invalid_arg ("UFO.order_dictionary_attrib: " ^ name) let coupling_dictionary_attrib ~strip name attribs = match find_attrib name attribs with | S.Coupling_Dictionary d -> List.map (fun (i, j, c) -> (i, j, name_to_string ~strip c)) d | _ -> invalid_arg ("UFO.coupling_dictionary_attrib: " ^ name) let decay_dictionary_attrib name attribs = match find_attrib name attribs with | S.Decay_Dictionary d -> List.map (fun (p, w) -> (List.map List.hd p, w)) d | _ -> invalid_arg ("UFO.decay_dictionary_attrib: " ^ name) (*i The following doesn't typecheck in applications, even with type annotations ... let attrib_handlers : type attribs value. string -> string -> attribs -> ((string -> attribs -> value) -> string -> value) * ((string -> attribs -> value) -> string -> value -> value) = fun kind symbol attribs -> let required query name = try query name attribs with | Not_found -> invalid_arg (Printf.sprintf "fatal UFO error: mandatory attribute `%s' missing for %s `%s'!" name kind symbol) and optional query name default = try query name attribs with | Not_found -> default in (required, optional) i*) let required_handler kind symbol attribs query name = try query name attribs with | Not_found -> invalid_arg (Printf.sprintf "fatal UFO error: mandatory attribute `%s' missing for %s `%s'!" name kind symbol) let optional_handler attribs query name default = try query name attribs with | Not_found -> default (* The UFO paper~\cite{Degrande:2011ua} is not clear on the question whether the \texttt{name} attribute of an instance must match its Python name. While the examples appear to imply this, there are examples of UFO files in the wild that violate this constraint. *) let warn_symbol_name file symbol name = if name <> symbol then Printf.eprintf "UFO: warning: symbol '%s' <> name '%s' in %s.py: \ while legal in UFO, it is unusual and can cause problems!\n" symbol name file let valid_fortran_id kind name = if not (ThoString.valid_fortran_id name) then invalid_arg (Printf.sprintf "fatal UFO error: the %s `%s' is not a valid fortran id!" kind name) let map_to_alist map = SMap.fold (fun key value acc -> (key, value) :: acc) map [] let keys map = SMap.fold (fun key _ acc -> key :: acc) map [] let keys_caseless map = CMap.fold (fun key _ acc -> key :: acc) map [] let values map = SMap.fold (fun _ value acc -> value :: acc) map [] module SKey = struct type t = string let hash = Hashtbl.hash let equal = (=) end module SHash = Hashtbl.Make (SKey) module type Particle = sig type t = private { pdg_code : int; name : string; antiname : string; spin : UFOx.Lorentz.r; color : UFOx.Color.r; mass : string; width : string; propagator : string option; texname : string; antitexname : string; charge : charge; ghost_number : int; lepton_number : int; y : charge; goldstone : bool; propagating : bool; (* NOT HANDLED YET! *) line : string option; (* NOT HANDLED YET! *) is_anti : bool } val of_file : S.t -> t SMap.t val to_string : string -> t -> string val conjugate : t -> t val force_spinor : t -> t val force_conjspinor : t -> t val force_majorana : t -> t val is_majorana : t -> bool val is_ghost : t -> bool val is_goldstone : t -> bool val is_physical : t -> bool val filter : (t -> bool) -> t SMap.t -> t SMap.t end module Particle : Particle = struct type t = { pdg_code : int; name : string; antiname : string; spin : UFOx.Lorentz.r; color : UFOx.Color.r; mass : string; width : string; propagator : string option; texname : string; antitexname : string; charge : charge; ghost_number : int; lepton_number : int; y : charge; goldstone : bool; propagating : bool; (* NOT HANDLED YET! *) line : string option; (* NOT HANDLED YET! *) is_anti : bool } let to_string symbol p = Printf.sprintf "particle: %s => [pdg = %d, name = '%s'/'%s', \ spin = %s, color = %s, \ mass = %s, width = %s,%s \ Q = %s, G = %d, L = %d, Y = %s, \ TeX = '%s'/'%s'%s]" symbol p.pdg_code p.name p.antiname (UFOx.Lorentz.rep_to_string p.spin) (UFOx.Color.rep_to_string p.color) p.mass p.width (match p.propagator with | None -> "" | Some p -> " propagator = " ^ p ^ ",") (charge_to_string p.charge) p.ghost_number p.lepton_number (charge_to_string p.y) p.texname p.antitexname (if p.goldstone then ", GB" else "") let conjugate_charge = function | Q_Integer i -> Q_Integer (-i) | Q_Fraction (n, d) -> Q_Fraction (-n, d) let is_neutral p = (p.name = p.antiname) (* We \emph{must not} mess with [pdg_code] and [color] if the particle is neutral! *) let conjugate p = if is_neutral p then p else { pdg_code = - p.pdg_code; name = p.antiname; antiname = p.name; spin = UFOx.Lorentz.rep_conjugate p.spin; color = UFOx.Color.rep_conjugate p.color; mass = p.mass; width = p.width; propagator = p.propagator; texname = p.antitexname; antitexname = p.texname; charge = conjugate_charge p.charge; ghost_number = p.ghost_number; lepton_number = p.lepton_number; y = p.y; goldstone = p.goldstone; propagating = p.propagating; line = p.line; is_anti = not p.is_anti } let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Particle" ], attribs -> let required query name = required_handler "particle" symbol attribs query name and optional query name default = optional_handler attribs query name default in let name = required string_attrib "name" and antiname = required string_attrib "antiname" in let neutral = (name = antiname) in SMap.add symbol { (* The required attributes per UFO docs. *) pdg_code = required integer_attrib "pdg_code"; name; antiname; spin = UFOx.Lorentz.rep_of_int neutral (required integer_attrib "spin"); color = UFOx.Color.rep_of_int neutral (required integer_attrib "color"); mass = required (name_attrib ~strip:"Param") "mass"; width = required (name_attrib ~strip:"Param") "width"; texname = required string_attrib "texname"; antitexname = required string_attrib "antitexname"; charge = required charge_attrib "charge"; (* The optional attributes per UFO docs. *) ghost_number = optional integer_attrib "GhostNumber" 0; lepton_number = optional integer_attrib "LeptonNumber" 0; y = optional charge_attrib "Y" (Q_Integer 0); goldstone = optional boolean_attrib "goldstone" false; propagating = optional boolean_attrib "propagating" true; line = (try Some (name_attrib "line" attribs) with _ -> None); (* Undocumented extensions. *) propagator = (try Some (name_attrib ~strip:"Prop" "propagator" attribs) with _ -> None); (* O'Mega extensions. *) is_anti = false } map | [ "anti"; p ], [] -> begin try SMap.add symbol (conjugate (SMap.find p map)) map with | Not_found -> invalid_arg ("Particle.of_file: " ^ p ^ ".anti() not yet defined!") end | _ -> invalid_arg ("Particle.of_file: " ^ name_to_string d.S.kind) let of_file particles = List.fold_left of_file1 SMap.empty particles let is_spinor p = match UFOx.Lorentz.omega p.spin with | Coupling.Spinor | Coupling.ConjSpinor | Coupling.Majorana -> true | _ -> false (* \begin{dubious} TODO: this is a bit of a hack: try to expose the type [UFOx.Lorentz_Atom'.r] instead. \end{dubious} *) let force_spinor p = if is_spinor p then { p with spin = UFOx.Lorentz.rep_of_int false 2 } else p let force_conjspinor p = if is_spinor p then { p with spin = UFOx.Lorentz.rep_of_int false (-2) } else p let force_majorana p = if is_spinor p then { p with spin = UFOx.Lorentz.rep_of_int true 2 } else p let is_majorana p = match UFOx.Lorentz.omega p.spin with | Coupling.Majorana | Coupling.Vectorspinor | Coupling.Maj_Ghost -> true | _ -> false let is_ghost p = p.ghost_number <> 0 let is_goldstone p = p.goldstone let is_physical p = not (is_ghost p || is_goldstone p) let filter predicate map = SMap.filter (fun symbol p -> predicate p) map end module type UFO_Coupling = sig type t = private { name : string; value : UFOx.Expr.t; order : (string * int) list } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module UFO_Coupling : UFO_Coupling = struct type t = { name : string; value : UFOx.Expr.t; order : (string * int) list } let order_to_string orders = String.concat ", " (List.map (fun (s, i) -> Printf.sprintf "'%s':%d" s i) orders) let to_string symbol c = Printf.sprintf "coupling: %s => [name = '%s', value = '%s', order = [%s]]" symbol c.name (expr_to_string c.value) (order_to_string c.order) let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Coupling" ], attribs -> let required query name = required_handler "coupling" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "couplings" symbol name; valid_fortran_id "coupling" name; SMap.add symbol { name; value = UFOx.Expr.of_string (required string_attrib "value"); order = required order_dictionary_attrib "order" } map | _ -> invalid_arg ("UFO_Coupling.of_file: " ^ name_to_string d.S.kind) let of_file couplings = List.fold_left of_file1 SMap.empty couplings end module type Coupling_Order = sig type t = private { name : string; expansion_order : int; hierarchy : int } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module Coupling_Order : Coupling_Order = struct type t = { name : string; expansion_order : int; hierarchy : int } let to_string symbol c = Printf.sprintf "coupling_order: %s => [name = '%s', \ expansion_order = '%d', \ hierarchy = %d]" symbol c.name c.expansion_order c.hierarchy let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "CouplingOrder" ], attribs -> let required query name = required_handler "coupling order" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "coupling_orders" symbol name; SMap.add symbol { name; expansion_order = required integer_attrib "expansion_order"; hierarchy = required integer_attrib "hierarchy" } map | _ -> invalid_arg ("Coupling_order.of_file: " ^ name_to_string d.S.kind) let of_file coupling_orders = List.fold_left of_file1 SMap.empty coupling_orders end module type Lorentz_UFO = sig (* If the \texttt{name} attribute of a \texttt{Lorentz} object does \emph{not} match the the name of the object, we need the latter for weeding out unused Lorentz structures (see [Vertex.contains] below). Therefore, we keep it around. *) type t = private { name : string; symbol : string; spins : int list; structure : UFOx.Lorentz.t } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module Lorentz_UFO : Lorentz_UFO = struct type t = { name : string; symbol : string; spins : int list; structure : UFOx.Lorentz.t } let to_string symbol l = Printf.sprintf "lorentz: %s => [name = '%s', spins = [%s], \ structure = %s]" symbol l.name (String.concat ", " (List.map string_of_int l.spins)) (UFOx.Lorentz.to_string l.structure) let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Lorentz" ], attribs -> let required query name = required_handler "lorentz" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "lorentz" symbol name; valid_fortran_id "lorentz" symbol; SMap.add symbol { name; symbol; spins = required integer_list_attrib "spins"; structure = UFOx.Lorentz.of_string (required string_attrib "structure") } map | _ -> invalid_arg ("Lorentz.of_file: " ^ name_to_string d.S.kind) let of_file lorentz = List.fold_left of_file1 SMap.empty lorentz end module type Vertex = sig type lcc = private (* Lorentz-color-coupling *) { lorentz : string; color : UFOx.Color.t; coupling : string } type t = private { name : string; particles : string array; lcc : lcc list } val of_file : Particle.t SMap.t -> S.t -> t SMap.t val to_string : string -> t -> string val to_string_expanded : Lorentz_UFO.t SMap.t -> UFO_Coupling.t SMap.t -> t -> string val contains : Particle.t SMap.t -> (Particle.t -> bool) -> t -> bool val filter : (t -> bool) -> t SMap.t -> t SMap.t end module Vertex : Vertex = struct type lcc = { lorentz : string; color : UFOx.Color.t; coupling : string } type t = { name : string; particles : string array; lcc : lcc list } let to_string symbol c = Printf.sprintf "vertex: %s => [name = '%s', particles = [%s], \ lorentz-color-couplings = [%s]" symbol c.name (String.concat ", " (Array.to_list c.particles)) (String.concat ", " (List.map (fun lcc -> Printf.sprintf "%s * %s * %s" lcc.coupling lcc.lorentz (UFOx.Color.to_string lcc.color)) c.lcc)) let to_string_expanded lorentz couplings c = let expand_lorentz s = try UFOx.Lorentz.to_string (SMap.find s lorentz).Lorentz_UFO.structure with | Not_found -> "?" in Printf.sprintf "expanded: [%s] -> { lorentz-color-couplings = [%s] }" (String.concat ", " (Array.to_list c.particles)) (String.concat ", " (List.map (fun lcc -> Printf.sprintf "%s * %s * %s" lcc.coupling (expand_lorentz lcc.lorentz) (UFOx.Color.to_string lcc.color)) c.lcc)) let contains particles predicate v = let p = v.particles in let rec contains' i = if i < 0 then false else if predicate (SMap.find p.(i) particles) then true else contains' (pred i) in contains' (Array.length p - 1) let force_adj_identity1 adj_indices = function | UFOx.Color_Atom.Identity (a, b) as atom -> begin match List.mem a adj_indices, List.mem b adj_indices with | true, true -> UFOx.Color_Atom.Identity8 (a, b) | false, false -> atom | true, false | false, true -> invalid_arg "force_adj_identity: mixed representations!" end | atom -> atom let force_adj_identity adj_indices tensor = UFOx.Color.map_atoms (force_adj_identity1 adj_indices) tensor let find_adj_indices map particles = let adj_indices = ref [] in Array.iteri (fun i p -> (* We must pattern match against the O'Mega representation, because [UFOx.Color.r] is abstract. *) match UFOx.Color.omega (SMap.find p map).Particle.color with | Color.AdjSUN _ -> adj_indices := succ i :: !adj_indices | _ -> ()) particles; !adj_indices let classify_color_indices map particles = let fund_indices = ref [] and conj_indices = ref [] and adj_indices = ref [] in Array.iteri (fun i p -> (* We must pattern match against the O'Mega representation, because [UFOx.Color.r] is abstract. *) match UFOx.Color.omega (SMap.find p map).Particle.color with | Color.SUN n -> if n > 0 then fund_indices := succ i :: !fund_indices else if n < 0 then conj_indices := succ i :: !conj_indices else failwith "classify_color_indices: SU(0)" | Color.AdjSUN n -> if n <> 0 then adj_indices := succ i :: !adj_indices else failwith "classify_color_indices: SU(0)" | _ -> ()) particles; (!fund_indices, !conj_indices, !adj_indices) (* FIXME: would have expected the opposite order \ldots *) let force_identity1 (fund_indices, conj_indices, adj_indices) = function | UFOx.Color_Atom.Identity (a, b) as atom -> if List.mem a fund_indices then begin if List.mem b conj_indices then UFOx.Color_Atom.Identity (b, a) else invalid_arg "force_adj_identity: mixed representations!" end else if List.mem a conj_indices then begin if List.mem b fund_indices then UFOx.Color_Atom.Identity (a, b) else invalid_arg "force_adj_identity: mixed representations!" end else if List.mem a adj_indices then begin if List.mem b adj_indices then UFOx.Color_Atom.Identity8 (a, b) else invalid_arg "force_adj_identity: mixed representations!" end else atom | atom -> atom let force_identity indices tensor = UFOx.Color.map_atoms (force_identity1 indices) tensor (* Here we don't have the Lorentz structures available yet. Thus we set [fermion_lines = []] for now and correct this later. *) let of_file1 particle_map map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Vertex" ], attribs -> let required query name = required_handler "vertex" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "vertices" symbol name; let particles = Array.of_list (required (name_list_attrib ~strip:"P") "particles") in let color = let indices = classify_color_indices particle_map particles in Array.of_list (List.map (force_identity indices @@ UFOx.Color.of_string) (required string_list_attrib "color")) and lorentz = Array.of_list (required (name_list_attrib ~strip:"L") "lorentz") and couplings_alist = required (coupling_dictionary_attrib ~strip:"C") "couplings" in let lcc = List.map (fun (i, j, c) -> { lorentz = lorentz.(j); color = color.(i); coupling = c }) couplings_alist in SMap.add symbol { name; particles; lcc } map | _ -> invalid_arg ("Vertex.of_file: " ^ name_to_string d.S.kind) let of_file particles vertices = List.fold_left (of_file1 particles) SMap.empty vertices let filter predicate map = SMap.filter (fun symbol p -> predicate p) map end module type Parameter = sig type nature = private Internal | External type ptype = private Real | Complex type t = private { name : string; nature : nature; ptype : ptype; value : value; texname : string; lhablock : string option; lhacode : int list option; sequence : int } val of_file : S.t -> t SMap.t val to_string : string -> t -> string val missing : string -> t end module Parameter : Parameter = struct type nature = Internal | External let nature_to_string = function | Internal -> "internal" | External -> "external" let nature_of_string = function | "internal" -> Internal | "external" -> External | s -> invalid_arg ("Parameter.nature_of_string: " ^ s) type ptype = Real | Complex let ptype_to_string = function | Real -> "real" | Complex -> "complex" let ptype_of_string = function | "real" -> Real | "complex" -> Complex | s -> invalid_arg ("Parameter.ptype_of_string: " ^ s) type t = { name : string; nature : nature; ptype : ptype; value : value; texname : string; lhablock : string option; lhacode : int list option; sequence : int } let to_string symbol p = Printf.sprintf "parameter: %s => [#%d, name = '%s', nature = %s, type = %s, \ value = %s, texname = '%s', \ lhablock = %s, lhacode = [%s]]" symbol p.sequence p.name (nature_to_string p.nature) (ptype_to_string p.ptype) (value_to_string p.value) p.texname (match p.lhablock with None -> "???" | Some s -> s) (match p.lhacode with | None -> "" | Some c -> String.concat ", " (List.map string_of_int c)) let of_file1 (map, n) d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Parameter" ], attribs -> let required query name = required_handler "particle" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "parameters" symbol name; valid_fortran_id "parameter" name; (SMap.add symbol { name; nature = nature_of_string (required string_attrib "nature"); ptype = ptype_of_string (required string_attrib "type"); value = required value_attrib "value"; texname = required string_attrib "texname"; lhablock = (try Some (string_attrib "lhablock" attribs) with Not_found -> None); lhacode = (try Some (integer_list_attrib "lhacode" attribs) with Not_found -> None); sequence = n } map, succ n) | _ -> invalid_arg ("Parameter.of_file: " ^ name_to_string d.S.kind) let of_file parameters = let map, _ = List.fold_left of_file1 (SMap.empty, 0) parameters in map let missing name = { name; nature = External; ptype = Real; value = Integer 0; texname = Printf.sprintf "\\texttt{%s}" name; lhablock = None; lhacode = None; sequence = 0 } end (* Macros are encoded as a special [S.declaration] with [S.kind = "$"]. This is slightly hackish, but general enough and the overhead of a special union type is probably not worth the effort. *) module type Macro = sig type t val empty : t (* The domains and codomains are still a bit too much ad hoc, but it does the job. *) val define : t -> string -> S.value -> t val expand_string : t -> string -> S.value val expand_expr : t -> S.string_atom list -> string (* Only for documentation: *) val expand_atom : t -> S.string_atom -> string end module Macro : Macro = struct type t = S.value SMap.t let empty = SMap.empty let define macros name expansion = SMap.add name expansion macros let expand_string macros name = SMap.find name macros let rec expand_atom macros = function | S.Literal s -> s | S.Macro [name] -> begin try begin match SMap.find name macros with | S.String s -> s | S.String_Expr expr -> expand_expr macros expr | _ -> invalid_arg ("expand_atom: not a string: " ^ name) end with | Not_found -> invalid_arg ("expand_atom: not found: " ^ name) end | S.Macro [] -> invalid_arg "expand_atom: empty" | S.Macro name -> invalid_arg ("expand_atom: compound name: " ^ String.concat "." name) and expand_expr macros expr = String.concat "" (List.map (expand_atom macros) expr) end module type Propagator_UFO = sig type t = (* private *) { name : string; numerator : UFOx.Lorentz.t; denominator : UFOx.Lorentz.t } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module Propagator_UFO : Propagator_UFO = struct type t = { name : string; numerator : UFOx.Lorentz.t; denominator : UFOx.Lorentz.t } let to_string symbol p = Printf.sprintf "propagator: %s => [name = '%s', numerator = '%s', \ denominator = '%s']" symbol p.name (UFOx.Lorentz.to_string p.numerator) (UFOx.Lorentz.to_string p.denominator) (* The \texttt{denominator} attribute is optional and there is a default (cf.~\texttt{arXiv:1308.1668}) *) let default_denominator = "P('mu', id) * P('mu', id) \ - Mass(id) * Mass(id) \ + complex(0,1) * Mass(id) * Width(id)" let of_string_with_error_correction symbol num_or_den s = try UFOx.Lorentz.of_string s with | Invalid_argument msg -> begin let fixed = s ^ ")" in try let tensor = UFOx.Lorentz.of_string fixed in Printf.eprintf "UFO.Propagator.of_string: added missing closing parenthesis \ in %s of %s: \"%s\"\n" num_or_den symbol s; tensor with | Invalid_argument _ -> invalid_arg (Printf.sprintf "UFO.Propagator.of_string: %s of %s: %s in \"%s\"\n" num_or_den symbol msg fixed) end let of_file1 (macros, map) d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Propagator" ], attribs -> let required query name = required_handler "particle" symbol attribs query name and optional query name default = optional_handler attribs query name default in let name = required string_attrib "name" in warn_symbol_name "propagators" symbol name; let num_string_expr = required string_expr_attrib "numerator" and den_string = begin match optional find_attrib "denominator" (S.String default_denominator) with | S.String s -> s | S.Name [n] -> begin match Macro.expand_string macros n with | S.String s -> s | _ -> invalid_arg "Propagator.denominator" end | _ -> invalid_arg "Propagator.denominator: " end in let num_string = Macro.expand_expr macros num_string_expr in let numerator = of_string_with_error_correction symbol "numerator" num_string and denominator = of_string_with_error_correction symbol "denominator" den_string in (macros, SMap.add symbol { name; numerator; denominator } map) | [ "$" ], [ macro ] -> begin match macro.S.a_value with | S.String _ as s -> (Macro.define macros symbol s, map); | S.String_Expr expr -> let expanded = S.String (Macro.expand_expr macros expr) in (Macro.define macros symbol expanded, map) | _ -> invalid_arg ("Propagator:of_file: not a string " ^ symbol) end | [ "$" ], [] -> invalid_arg ("Propagator:of_file: empty declaration " ^ symbol) | [ "$" ], _ -> invalid_arg ("Propagator:of_file: multiple declaration " ^ symbol) | _ -> invalid_arg ("Propagator:of_file: " ^ name_to_string d.S.kind) let of_file propagators = let _, propagators' = List.fold_left of_file1 (Macro.empty, SMap.empty) propagators in propagators' end module type Decay = sig type t = private { name : string; particle : string; widths : (string list * string) list } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module Decay : Decay = struct type t = { name : string; particle : string; widths : (string list * string) list } let width_to_string ws = String.concat ", " (List.map (fun (ps, w) -> "(" ^ String.concat ", " ps ^ ") -> '" ^ w ^ "'") ws) let to_string symbol d = Printf.sprintf "decay: %s => [name = '%s', particle = '%s', widths = [%s]]" symbol d.name d.particle (width_to_string d.widths) let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Decay" ], attribs -> let required query name = required_handler "particle" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "decays" symbol name; SMap.add symbol { name; particle = required (name_attrib ~strip:"P") "particle"; widths = required decay_dictionary_attrib "partial_widths" } map | _ -> invalid_arg ("Decay.of_file: " ^ name_to_string d.S.kind) let of_file decays = List.fold_left of_file1 SMap.empty decays end (* We can read the spinor representations off the vertices to check for consistency. *) (* \begin{dubious} Note that we have to conjugate the representations! \end{dubious} *) let collect_spinor_reps_of_vertex particles lorentz v sets = List.fold_left (fun sets' lcc -> let l = (SMap.find lcc.Vertex.lorentz lorentz).Lorentz_UFO.structure in List.fold_left (fun (spinors, conj_spinors as sets'') (i, rep) -> let p = v.Vertex.particles.(pred i) in match UFOx.Lorentz.omega rep with | Coupling.ConjSpinor -> (SSet.add p spinors, conj_spinors) | Coupling.Spinor -> (spinors, SSet.add p conj_spinors) | _ -> sets'') sets' (UFOx.Lorentz.classify_indices l)) sets v.Vertex.lcc let collect_spinor_reps_of_vertices particles lorentz vertices = SMap.fold (fun _ v -> collect_spinor_reps_of_vertex particles lorentz v) vertices (SSet.empty, SSet.empty) let lorentz_reps_of_vertex particles v = ThoList.alist_of_list ~predicate:(not @@ UFOx.Lorentz.rep_trivial) ~offset:1 (List.map (fun p -> (* Why do we need to conjugate??? *) UFOx.Lorentz.rep_conjugate (SMap.find p particles).Particle.spin) (Array.to_list v.Vertex.particles)) let rep_compatible rep_vertex rep_particle = let open UFOx.Lorentz in let open Coupling in match omega rep_vertex, omega rep_particle with | (Spinor | ConjSpinor), Majorana -> true | r1, r2 -> r1 = r2 let reps_compatible reps_vertex reps_particles = List.for_all2 (fun (iv, rv) (ip, rp) -> iv = ip && rep_compatible rv rp) reps_vertex reps_particles let check_lorentz_reps_of_vertex particles lorentz v = let reps_particles = List.sort compare (lorentz_reps_of_vertex particles v) in List.iter (fun lcc -> let l = (SMap.find lcc.Vertex.lorentz lorentz).Lorentz_UFO.structure in let reps_vertex = List.sort compare (UFOx.Lorentz.classify_indices l) in if not (reps_compatible reps_vertex reps_particles) then begin Printf.eprintf "%s <> %s [%s]\n" (UFOx.Index.classes_to_string UFOx.Lorentz.rep_to_string reps_particles) (UFOx.Index.classes_to_string UFOx.Lorentz.rep_to_string reps_vertex) v.Vertex.name (* [(Vertex.to_string v.Vertex.name v)] *); (* [invalid_arg "check_lorentz_reps_of_vertex"] *) () end) v.Vertex.lcc let color_reps_of_vertex particles v = ThoList.alist_of_list ~predicate:(not @@ UFOx.Color.rep_trivial) ~offset:1 (List.map (fun p -> (SMap.find p particles).Particle.color) (Array.to_list v.Vertex.particles)) let check_color_reps_of_vertex particles v = let reps_particles = List.sort compare (color_reps_of_vertex particles v) in List.iter (fun lcc -> let reps_vertex = List.sort compare (UFOx.Color.classify_indices lcc.Vertex.color) in if reps_vertex <> reps_particles then begin Printf.printf "%s <> %s\n" (UFOx.Index.classes_to_string UFOx.Color.rep_to_string reps_particles) (UFOx.Index.classes_to_string UFOx.Color.rep_to_string reps_vertex); invalid_arg "check_color_reps_of_vertex" end) v.Vertex.lcc module P = Permutation.Default module type Lorentz = sig type spins = private | Unused | Unique of Coupling.lorentz array | Ambiguous of Coupling.lorentz array SMap.t type t = private { name : string; n : int; spins : spins; structure : UFO_Lorentz.t; fermion_lines : Coupling.fermion_lines; variables : string list } val required_charge_conjugates : t -> t list val permute : P.t -> t -> t val of_lorentz_UFO : Particle.t SMap.t -> Vertex.t SMap.t -> Lorentz_UFO.t SMap.t -> t SMap.t val lorentz_to_string : Coupling.lorentz -> string val to_string : string -> t -> string end module Lorentz : Lorentz = struct let rec lorentz_to_string = function | Coupling.Scalar -> "Scalar" | Coupling.Spinor -> "Spinor" | Coupling.ConjSpinor -> "ConjSpinor" | Coupling.Majorana -> "Majorana" | Coupling.Maj_Ghost -> "Maj_Ghost" | Coupling.Vector -> "Vector" | Coupling.Massive_Vector -> "Massive_Vector" | Coupling.Vectorspinor -> "Vectorspinor" | Coupling.Tensor_1 -> "Tensor_1" | Coupling.Tensor_2 -> "Tensor_2" | Coupling.BRS l -> "BRS(" ^ lorentz_to_string l ^ ")" (* Unlike UFO, O'Mega distinguishes bewteen spinors and conjugate spinors. However, we can inspect the particles in the vertices in which a Lorentz structure is used to determine the correct quantum numbers. Most model files in the real world contain unused Lorentz structures. This is not a problem, we can just ignore them. *) type spins = | Unused | Unique of Coupling.lorentz array | Ambiguous of Coupling.lorentz array SMap.t (* \begin{dubious} Use [UFO_targets.Fortran.fusion_name] below in order to avoid communication problems. Or even move away from strings alltogether. \end{dubious} *) type t = { name : string; n : int; spins : spins; structure : UFO_Lorentz.t; fermion_lines : Coupling.fermion_lines; variables : string list } (* Add one charge conjugated fermion lines. *) let charge_conjugate1 l (ket, bra as fermion_line) = { name = l.name ^ Printf.sprintf "_c%x%x" ket bra; n = l.n; spins = l.spins; structure = UFO_Lorentz.charge_conjugate fermion_line l.structure; fermion_lines = l.fermion_lines; variables = l.variables } (* Add several charge conjugated fermion lines. *) let charge_conjugate l fermion_lines = List.fold_left charge_conjugate1 l fermion_lines (*i let all_charge_conjugates l = List.map (charge_conjugate l) (ThoList.power l.fermion_lines) i*) (* Add all combinations of charge conjugated fermion lines that don't leave the fusion. *) let required_charge_conjugates l = let saturated_fermion_lines = List.filter (fun (ket, bra) -> ket != 1 && bra != 1) l.fermion_lines in List.map (charge_conjugate l) (ThoList.power saturated_fermion_lines) let permute_spins p = function | Unused -> Unused | Unique s -> Unique (P.array p s) | Ambiguous map -> Ambiguous (SMap.map (P.array p) map) (* Note that we apply the \emph{inverse} permutation to the indices in order to match the permutation of the particles/spins. *) let permute_structure n p (l, f) = let permuted = P.array (P.inverse p) (Array.init n succ) in let permute_index i = if i > 0 then UFOx.Index.map_position (fun pos -> permuted.(pred pos)) i else i in (UFO_Lorentz.map_indices permute_index l, UFO_Lorentz.map_fermion_lines permute_index f) let permute p l = let structure, fermion_lines = permute_structure l.n p (l.structure, l.fermion_lines) in { name = l.name ^ "_p" ^ P.to_string (P.inverse p); n = l.n; spins = permute_spins p l.spins; structure; fermion_lines; variables = l.variables } let omega_lorentz_reps n alist = let reps = Array.make n Coupling.Scalar in List.iter (fun (i, rep) -> reps.(pred i) <- UFOx.Lorentz.omega rep) alist; reps let contained lorentz vertex = List.exists (fun lcc1 -> lcc1.Vertex.lorentz = lorentz.Lorentz_UFO.symbol) vertex.Vertex.lcc (* Find all vertices in with the Lorentz structure [lorentz] is used and build a map from those vertices to the O'Mega Lorentz representations inferred from UFO's Lorentz structure and the [particles] involved. Then scan the bindings and check that we have inferred the same Lorentz representation from all vertices. *) let lorentz_reps_of_structure particles vertices lorentz = let uses = SMap.fold (fun name v acc -> if contained lorentz v then SMap.add name (omega_lorentz_reps (Array.length v.Vertex.particles) (lorentz_reps_of_vertex particles v)) acc else acc) vertices SMap.empty in let variants = ThoList.uniq (List.sort compare (List.map snd (SMap.bindings uses))) in match variants with | [] -> Unused | [s] -> Unique s | _ -> Printf.eprintf "UFO.Lorentz.lorentz_reps_of_structure: AMBIGUOUS!\n"; List.iter (fun variant -> Printf.eprintf "UFO.Lorentz.lorentz_reps_of_structure: %s\n" (ThoList.to_string lorentz_to_string (Array.to_list variant))) variants; Ambiguous uses let of_lorentz_tensor spins lorentz = match spins with | Unique s -> begin try Some (UFO_Lorentz.parse (Array.to_list s) lorentz) with | Failure msg -> begin prerr_endline msg; Some (UFO_Lorentz.dummy) end end | Unused -> Printf.eprintf "UFO.Lorentz: stripping unused structure %s\n" (UFOx.Lorentz.to_string lorentz); None | Ambiguous _ -> invalid_arg "UFO.Lorentz.of_lorentz_tensor: Ambiguous" (* NB: if the \texttt{name} attribute of a \texttt{Lorentz} object does \emph{not} match the the name of the object, the former has a better chance to correspond to a valid Fortran name. Therefore we use it. *) let of_lorentz_UFO particles vertices lorentz_UFO = SMap.fold (fun name l acc -> let spins = lorentz_reps_of_structure particles vertices l in match of_lorentz_tensor spins l.Lorentz_UFO.structure with | None -> acc | Some structure -> SMap.add name { name = l.Lorentz_UFO.symbol; n = List.length l.Lorentz_UFO.spins; spins; structure; fermion_lines = UFO_Lorentz.fermion_lines structure; variables = UFOx.Lorentz.variables l.Lorentz_UFO.structure } acc) lorentz_UFO SMap.empty let to_string symbol l = Printf.sprintf "lorentz: %s => [name = '%s', spins = %s, \ structure = %s, fermion_lines = %s]" symbol l.name (match l.spins with | Unique s -> "[" ^ String.concat ", " (List.map lorentz_to_string (Array.to_list s)) ^ "]" | Ambiguous _ -> "AMBIGUOUS!" | Unused -> "UNUSED!") (UFO_Lorentz.to_string l.structure) (UFO_Lorentz.fermion_lines_to_string l.fermion_lines) end (* According to arxiv:1308:1668, there should not be a factor of~$i$ in the numerators of propagators, but the (unused) \texttt{propagators.py} in most models violate this rule! *) let divide_propagators_by_i = ref false module type Propagator = sig type t = (* private *) { name : string; spins : Coupling.lorentz * Coupling.lorentz; numerator : UFO_Lorentz.t; denominator : UFO_Lorentz.t; variables : string list } val of_propagator_UFO : ?majorana:bool -> Propagator_UFO.t -> t val of_propagators_UFO : ?majorana:bool -> Propagator_UFO.t SMap.t -> t SMap.t val transpose : t -> t val to_string : string -> t -> string end module Propagator : Propagator = struct type t = (* private *) { name : string; spins : Coupling.lorentz * Coupling.lorentz; numerator : UFO_Lorentz.t; denominator : UFO_Lorentz.t; variables : string list } let lorentz_rep_at rep_classes i = try UFOx.Lorentz.omega (List.assoc i rep_classes) with | Not_found -> Coupling.Scalar let imaginary = Algebra.QC.make Algebra.Q.null Algebra.Q.unit let scalars = [Coupling.Scalar; Coupling.Scalar] (* If~$51$ and~$52$ show up as indices, we must map $(1,51)\to(1001,2001)$ and $(2,52)\to(1002,2002)$, as per the UFO conventions for Lorentz structures. *) (* \begin{dubious} This does not work yet, because [UFOx.Lorentz.map_indices] affects also the position argument of [P], [Mass] and [Width]. \end{dubious} *) let contains_51_52 tensor = List.exists (fun (i, _) -> i = 51 || i = 52) (UFOx.Lorentz.classify_indices tensor) let remap_51_52 = function | 1 -> 1001 | 51 -> 2001 | 2 -> 1002 | 52 -> 2002 | i -> i let canonicalize_51_52 tensor = if contains_51_52 tensor then UFOx.Lorentz.rename_indices remap_51_52 tensor else tensor let force_majorana = function | Coupling.Spinor | Coupling.ConjSpinor -> Coupling.Majorana | s -> s let string_list_union l1 l2 = Sets.String.elements (Sets.String.union (Sets.String.of_list l1) (Sets.String.of_list l2)) (* In the current conventions, the factor of~$i$ is not included: *) let of_propagator_UFO ?(majorana=false) p = let numerator = canonicalize_51_52 p.Propagator_UFO.numerator in let lorentz_reps = UFOx.Lorentz.classify_indices numerator in let spin1 = lorentz_rep_at lorentz_reps 1 and spin2 = lorentz_rep_at lorentz_reps 2 in let numerator_sans_i = if !divide_propagators_by_i then UFOx.Lorentz.map_coeff (fun q -> Algebra.QC.div q imaginary) numerator else numerator in { name = p.Propagator_UFO.name; spins = if majorana then (force_majorana spin1, force_majorana spin2) else (spin1, spin2); numerator = UFO_Lorentz.parse ~allow_denominator:true [spin1; spin2] numerator_sans_i; denominator = UFO_Lorentz.parse scalars p.Propagator_UFO.denominator; variables = string_list_union (UFOx.Lorentz.variables p.Propagator_UFO.denominator) (UFOx.Lorentz.variables numerator_sans_i) } let of_propagators_UFO ?majorana propagators_UFO = SMap.fold (fun name p acc -> SMap.add name (of_propagator_UFO ?majorana p) acc) propagators_UFO SMap.empty let permute12 = function | 1 -> 2 | 2 -> 1 | n -> n let transpose_positions t = UFOx.Index.map_position permute12 t let transpose p = { name = p.name; spins = (snd p.spins, fst p.spins); numerator = UFO_Lorentz.map_indices transpose_positions p.numerator; denominator = p.denominator; variables = p.variables } let to_string symbol p = Printf.sprintf "propagator: %s => [name = '%s', spin = '(%s, %s)', numerator/I = '%s', \ denominator = '%s']" symbol p.name (Lorentz.lorentz_to_string (fst p.spins)) (Lorentz.lorentz_to_string (snd p.spins)) (UFO_Lorentz.to_string p.numerator) (UFO_Lorentz.to_string p.denominator) end type t = { particles : Particle.t SMap.t; particle_array : Particle.t array; (* for diagnostics *) couplings : UFO_Coupling.t SMap.t; coupling_orders : Coupling_Order.t SMap.t; vertices : Vertex.t SMap.t; lorentz_UFO : Lorentz_UFO.t SMap.t; lorentz : Lorentz.t SMap.t; parameters : Parameter.t SMap.t; propagators_UFO : Propagator_UFO.t SMap.t; propagators : Propagator.t SMap.t; decays : Decay.t SMap.t; nc : int } let use_majorana_spinors = ref false let fallback_to_majorana_if_necessary particles vertices lorentz_UFO = let majoranas = SMap.fold (fun p particle acc -> if Particle.is_majorana particle then SSet.add p acc else acc) particles SSet.empty in let spinors, conj_spinors = collect_spinor_reps_of_vertices particles lorentz_UFO vertices in let ambiguous = SSet.diff (SSet.inter spinors conj_spinors) majoranas in let no_majoranas = SSet.is_empty majoranas and no_ambiguities = SSet.is_empty ambiguous in if no_majoranas && no_ambiguities && not !use_majorana_spinors then (SMap.mapi (fun p particle -> if SSet.mem p spinors then Particle.force_spinor particle else if SSet.mem p conj_spinors then Particle.force_conjspinor particle else particle) particles, false) else begin if !use_majorana_spinors then Printf.eprintf "O'Mega: Majorana fermions requested.\n"; if not no_majoranas then Printf.eprintf "O'Mega: found Majorana fermions!\n"; if not no_ambiguities then Printf.eprintf "O'Mega: found ambiguous spinor representations for %s!\n" (String.concat ", " (SSet.elements ambiguous)); Printf.eprintf "O'Mega: falling back to the Majorana representation for all fermions.\n"; (SMap.map Particle.force_majorana particles, true) end let nc_of_particles particles = let nc_set = List.fold_left (fun nc_set (_, p) -> match UFOx.Color.omega p.Particle.color with | Color.Singlet -> nc_set | Color.SUN nc -> Sets.Int.add (abs nc) nc_set | Color.AdjSUN nc -> Sets.Int.add (abs nc) nc_set) Sets.Int.empty (SMap.bindings particles) in match Sets.Int.elements nc_set with | [] -> 0 | [n] -> n | nc_list -> invalid_arg ("UFO.Model: more than one value of N_C: " ^ String.concat ", " (List.map string_of_int nc_list)) let of_file u = let particles = Particle.of_file u.Files.particles in let vertices = Vertex.of_file particles u.Files.vertices and lorentz_UFO = Lorentz_UFO.of_file u.Files.lorentz and propagators_UFO = Propagator_UFO.of_file u.Files.propagators in let particles, majorana = fallback_to_majorana_if_necessary particles vertices lorentz_UFO in let particle_array = Array.of_list (values particles) and lorentz = Lorentz.of_lorentz_UFO particles vertices lorentz_UFO and propagators = Propagator.of_propagators_UFO ~majorana propagators_UFO in let model = { particles; particle_array; couplings = UFO_Coupling.of_file u.Files.couplings; coupling_orders = Coupling_Order.of_file u.Files.coupling_orders; vertices; lorentz_UFO; lorentz; parameters = Parameter.of_file u.Files.parameters; propagators_UFO; propagators; decays = Decay.of_file u.Files.decays; nc = nc_of_particles particles } in SMap.iter (fun _ v -> check_color_reps_of_vertex model.particles v; check_lorentz_reps_of_vertex model.particles model.lorentz_UFO v) model.vertices; model let parse_directory dir = of_file (Files.parse_directory dir) let dump model = Printf.printf "NC = %d\n" model.nc; SMap.iter (print_endline @@@ Particle.to_string) model.particles; SMap.iter (print_endline @@@ UFO_Coupling.to_string) model.couplings; SMap.iter (print_endline @@@ Coupling_Order.to_string) model.coupling_orders; (* [SMap.iter (print_endline @@@ Vertex.to_string) model.vertices;] *) SMap.iter (fun symbol v -> (print_endline @@@ Vertex.to_string) symbol v; print_endline (Vertex.to_string_expanded model.lorentz_UFO model.couplings v)) model.vertices; SMap.iter (print_endline @@@ Lorentz_UFO.to_string) model.lorentz_UFO; SMap.iter (print_endline @@@ Lorentz.to_string) model.lorentz; SMap.iter (print_endline @@@ Parameter.to_string) model.parameters; SMap.iter (print_endline @@@ Propagator_UFO.to_string) model.propagators_UFO; SMap.iter (print_endline @@@ Propagator.to_string) model.propagators; SMap.iter (print_endline @@@ Decay.to_string) model.decays; SMap.iter (fun symbol d -> List.iter (fun (_, w) -> ignore (UFOx.Expr.of_string w)) d.Decay.widths) model.decays exception Unhandled of string let unhandled s = raise (Unhandled s) module Model = struct (* NB: we could use [type flavor = Particle.t], but that would be very inefficient, because we will use [flavor] as a key for maps below. *) type flavor = int type constant = string type gauge = unit module M = Modeltools.Mutable (struct type f = flavor type g = gauge type c = constant end) let flavors = M.flavors let external_flavors = M.external_flavors let external_flavors = M.external_flavors let lorentz = M.lorentz let color = M.color let nc = M.nc let propagator = M.propagator let width = M.width let goldstone = M.goldstone let conjugate = M.conjugate let fermion = M.fermion let vertices = M.vertices let fuse2 = M.fuse2 let fuse3 = M.fuse3 let fuse = M.fuse let max_degree = M.max_degree let parameters = M.parameters let flavor_of_string = M.flavor_of_string let flavor_to_string = M.flavor_to_string let flavor_to_TeX = M.flavor_to_TeX let flavor_symbol = M.flavor_symbol let gauge_symbol = M.gauge_symbol let pdg = M.pdg let mass_symbol = M.mass_symbol let width_symbol = M.width_symbol let constant_symbol = M.constant_symbol module Ch = M.Ch let charges = M.charges let rec fermion_of_lorentz = function | Coupling.Spinor -> 1 | Coupling.ConjSpinor -> -1 | Coupling.Majorana -> 2 | Coupling.Maj_Ghost -> 2 | Coupling.Vectorspinor -> 1 | Coupling.Vector | Coupling.Massive_Vector -> 0 | Coupling.Scalar | Coupling.Tensor_1 | Coupling.Tensor_2 -> 0 | Coupling.BRS f -> fermion_of_lorentz f module Q = Algebra.Q module QC = Algebra.QC let dummy_tensor3 = Coupling.Scalar_Scalar_Scalar 1 let dummy_tensor4 = Coupling.Scalar4 1 let triplet p = (p.(0), p.(1), p.(2)) let quartet p = (p.(0), p.(1), p.(2), p.(3)) let half_times q1 q2 = Q.mul (Q.make 1 2) (Q.mul q1 q2) let name g = g.UFO_Coupling.name let fractional_coupling g r = let g = name g in match Q.to_ratio r with | 0, _ -> "0.0_default" | 1, 1 -> g | -1, 1 -> Printf.sprintf "(-%s)" g | n, 1 -> Printf.sprintf "(%d*%s)" n g | 1, d -> Printf.sprintf "(%s/%d)" g d | -1, d -> Printf.sprintf "(-%s/%d)" g d | n, d -> Printf.sprintf "(%d*%s/%d)" n g d let lorentz_of_symbol model symbol = try SMap.find symbol model.lorentz with | Not_found -> invalid_arg ("lorentz_of_symbol: " ^ symbol) let lorentz_UFO_of_symbol model symbol = try SMap.find symbol model.lorentz_UFO with | Not_found -> invalid_arg ("lorentz_UFO_of_symbol: " ^ symbol) let coupling_of_symbol model symbol = try SMap.find symbol model.couplings with | Not_found -> invalid_arg ("coupling_of_symbol: " ^ symbol) let spin_triplet model name = match (lorentz_of_symbol model name).Lorentz.spins with | Lorentz.Unique [|s0; s1; s2|] -> (s0, s1, s2) | Lorentz.Unique _ -> invalid_arg "spin_triplet: wrong number of spins" | Lorentz.Unused -> invalid_arg "spin_triplet: Unused" | Lorentz.Ambiguous _ -> invalid_arg "spin_triplet: Ambiguous" let spin_quartet model name = match (lorentz_of_symbol model name).Lorentz.spins with | Lorentz.Unique [|s0; s1; s2; s3|] -> (s0, s1, s2, s3) | Lorentz.Unique _ -> invalid_arg "spin_quartet: wrong number of spins" | Lorentz.Unused -> invalid_arg "spin_quartet: Unused" | Lorentz.Ambiguous _ -> invalid_arg "spin_quartet: Ambiguous" let spin_multiplet model name = match (lorentz_of_symbol model name).Lorentz.spins with | Lorentz.Unique sarray -> sarray | Lorentz.Unused -> invalid_arg "spin_multiplet: Unused" | Lorentz.Ambiguous _ -> invalid_arg "spin_multiplet: Ambiguous" (* If we have reason to belive that a $\delta_{ab}$-vertex is an effective $\tr(T_aT_b)$-vertex generated at loop level, like~$gg\to H\ldots$ in the SM, we should interpret it as such and use the expression~(6.2) from~\cite{Kilian:2012pz}. *) (* AFAIK, there is no way to distinguish these cases directly in a UFO file. Instead we rely in a heuristic, in which each massless color octet vector particle or ghost is a gluon and colorless scalars are potential Higgses. *) let is_massless p = match ThoString.uppercase p.Particle.mass with | "ZERO" -> true | _ -> false let is_gluon model f = let p = model.particle_array.(f) in match UFOx.Color.omega p.Particle.color, UFOx.Lorentz.omega p.Particle.spin with | Color.AdjSUN _, Coupling.Vector -> is_massless p | Color.AdjSUN _, Coupling.Scalar -> if p.Particle.ghost_number <> 0 then is_massless p else false | _ -> false let is_color_singlet model f = let p = model.particle_array.(f) in match UFOx.Color.omega p.Particle.color with | Color.Singlet -> true | _ -> false let is_higgs_gluon_vertex model p adjoints = if Array.length p > List.length adjoints then List.for_all (fun (i, p) -> if List.mem i adjoints then is_gluon model p else is_color_singlet model p) (ThoList.enumerate 1 (Array.to_list p)) else false let delta8_heuristics model p a b = if is_higgs_gluon_vertex model p [a; b] then Color.Vertex.delta8_loop a b else Color.Vertex.delta8 a b let verbatim_higgs_glue = ref false let translate_color_atom model p = function | UFOx.Color_Atom.Identity (i, j) -> Color.Vertex.delta3 i j | UFOx.Color_Atom.Identity8 (a, b) -> if !verbatim_higgs_glue then Color.Vertex.delta8 a b else delta8_heuristics model p a b | UFOx.Color_Atom.T (a, i, j) -> Color.Vertex.t a i j | UFOx.Color_Atom.F (a, b, c) -> Color.Vertex.f a b c | UFOx.Color_Atom.D (a, b, c) -> Color.Vertex.d a b c | UFOx.Color_Atom.Epsilon (i, j, k) -> Color.Vertex.epsilon i j k | UFOx.Color_Atom.EpsilonBar (i, j, k) -> Color.Vertex.epsilonbar i j k | UFOx.Color_Atom.T6 (a, i, j) -> Color.Vertex.t6 a i j | UFOx.Color_Atom.K6 (i, j, k) -> Color.Vertex.k6 i j k | UFOx.Color_Atom.K6Bar (i, j, k) -> Color.Vertex.k6bar i j k let translate_color_term model p = function | [], q -> Color.Vertex.scale q Color.Vertex.unit | [atom], q -> Color.Vertex.scale q (translate_color_atom model p atom) | atoms, q -> let atoms = List.map (translate_color_atom model p) atoms in Color.Vertex.scale q (Color.Vertex.multiply atoms) let translate_color model p terms = match terms with | [] -> invalid_arg "translate_color: empty" | [ term ] -> translate_color_term model p term | terms -> Color.Vertex.sum (List.map (translate_color_term model p) terms) let translate_coupling_1 model p lcc = let l = lcc.Vertex.lorentz in let s = Array.to_list (spin_multiplet model l) and fl = (SMap.find l model.lorentz).Lorentz.fermion_lines and c = name (coupling_of_symbol model lcc.Vertex.coupling) in match lcc.Vertex.color with | UFOx.Color.Linear color -> let col = translate_color model p color in (Array.to_list p, Coupling.UFO (QC.unit, l, s, fl, col), c) | UFOx.Color.Ratios _ as color -> invalid_arg ("UFO.Model.translate_coupling: invalid color structure" ^ UFOx.Color.to_string color) let translate_coupling model p lcc = List.map (translate_coupling_1 model p) lcc let long_flavors = ref false module type Lookup = sig type f = private { flavors : flavor list; flavor_of_string : string -> flavor; flavor_of_symbol : string -> flavor; particle : flavor -> Particle.t; flavor_symbol : flavor -> string; conjugate : flavor -> flavor } type flavor_format = | Long | Decimal | Hexadecimal val flavor_format : flavor_format ref val of_model : t -> f end module Lookup : Lookup = struct type f = { flavors : flavor list; flavor_of_string : string -> flavor; flavor_of_symbol : string -> flavor; particle : flavor -> Particle.t; flavor_symbol : flavor -> string; conjugate : flavor -> flavor } type flavor_format = | Long | Decimal | Hexadecimal let flavor_format = ref Hexadecimal let conjugate_of_particle_array particles = Array.init (Array.length particles) (fun i -> let f' = Particle.conjugate particles.(i) in match ThoArray.match_all f' particles with | [i'] -> i' | [] -> invalid_arg ("no charge conjugate: " ^ f'.Particle.name) | _ -> invalid_arg ("multiple charge conjugates: " ^ f'.Particle.name)) let invert_flavor_array a = let table = SHash.create 37 in Array.iteri (fun i s -> SHash.add table s i) a; (fun name -> try SHash.find table name with | Not_found -> invalid_arg ("not found: " ^ name)) let digits base n = let rec digits' acc n = if n < 1 then acc else digits' (succ acc) (n / base) in if n < 0 then digits' 1 (-n) else if n = 0 then 1 else digits' 0 n let of_model model = let particle_array = Array.of_list (values model.particles) in let conjugate_array = conjugate_of_particle_array particle_array and name_array = Array.map (fun f -> f.Particle.name) particle_array and symbol_array = Array.of_list (keys model.particles) in let flavor_symbol f = begin match !flavor_format with | Long -> symbol_array.(f) | Decimal -> let w = digits 10 (Array.length particle_array - 1) in Printf.sprintf "%0*d" w f | Hexadecimal -> let w = digits 16 (Array.length particle_array - 1) in Printf.sprintf "%0*X" w f end in { flavors = ThoList.range 0 (Array.length particle_array - 1); flavor_of_string = invert_flavor_array name_array; flavor_of_symbol = invert_flavor_array symbol_array; particle = Array.get particle_array; flavor_symbol = flavor_symbol; conjugate = Array.get conjugate_array } end (* \begin{dubious} We appear to need to conjugate all flavors. Why??? \end{dubious} *) let translate_vertices model tables = let vn = List.fold_left (fun acc v -> let p = Array.map tables.Lookup.flavor_of_symbol v.Vertex.particles and lcc = v.Vertex.lcc in let p = Array.map conjugate p in (* FIXME: why? *) translate_coupling model p lcc @ acc) [] (values model.vertices) in ([], [], vn) let propagator_of_lorentz = function | Coupling.Scalar -> Coupling.Prop_Scalar | Coupling.Spinor -> Coupling.Prop_Spinor | Coupling.ConjSpinor -> Coupling.Prop_ConjSpinor | Coupling.Majorana -> Coupling.Prop_Majorana | Coupling.Maj_Ghost -> invalid_arg "UFO.Model.propagator_of_lorentz: SUSY ghosts do not propagate" | Coupling.Vector -> Coupling.Prop_Feynman | Coupling.Massive_Vector -> Coupling.Prop_Unitarity | Coupling.Tensor_2 -> Coupling.Prop_Tensor_2 | Coupling.Vectorspinor -> invalid_arg "UFO.Model.propagator_of_lorentz: Vectorspinor" | Coupling.Tensor_1 -> invalid_arg "UFO.Model.propagator_of_lorentz: Tensor_1" | Coupling.BRS _ -> invalid_arg "UFO.Model.propagator_of_lorentz: no BRST" let filter_unphysical model = let physical_particles = Particle.filter Particle.is_physical model.particles in let physical_particle_array = Array.of_list (values physical_particles) in let physical_vertices = Vertex.filter (not @@ (Vertex.contains model.particles (not @@ Particle.is_physical))) model.vertices in { model with particles = physical_particles; particle_array = physical_particle_array; vertices = physical_vertices } let whizard_constants = SSet.of_list [ "ZERO" ] let filter_constants parameters = List.filter (fun p -> not (SSet.mem (ThoString.uppercase p.Parameter.name) whizard_constants)) parameters let add_name set parameter = CSet.add parameter.Parameter.name set let hardcoded_parameters = CSet.of_list ["cmath.pi"] let missing_parameters input derived couplings = let input_parameters = List.fold_left add_name hardcoded_parameters input in let all_parameters = List.fold_left add_name input_parameters derived in let derived_dependencies = dependencies (List.map (fun p -> (p.Parameter.name, p.Parameter.value)) derived) in let coupling_dependencies = dependencies (List.map (fun p -> (p.UFO_Coupling.name, Expr p.UFO_Coupling.value)) (values couplings)) in let missing_input = CMap.filter (fun parameter derived_parameters -> not (CSet.mem parameter all_parameters)) derived_dependencies and missing = CMap.filter (fun parameter couplings -> not (CSet.mem parameter all_parameters)) coupling_dependencies in CMap.iter (fun parameter derived_parameters -> Printf.eprintf "UFO warning: undefined input parameter %s appears in derived \ parameters {%s}: will be added to the list of input parameters!\n" parameter (String.concat "; " (CSet.elements derived_parameters))) missing_input; CMap.iter (fun parameter couplings -> Printf.eprintf "UFO warning: undefined parameter %s appears in couplings {%s}: \ will be added to the list of input parameters!\n" parameter (String.concat "; " (CSet.elements couplings))) missing; keys_caseless missing_input @ keys_caseless missing let classify_parameters model = let compare_parameters p1 p2 = compare p1.Parameter.sequence p2.Parameter.sequence in let input, derived = List.fold_left (fun (input, derived) p -> match p.Parameter.nature with | Parameter.Internal -> (input, p :: derived) | Parameter.External -> begin match p.Parameter.ptype with | Parameter.Real -> () | Parameter.Complex -> Printf.eprintf "UFO warning: invalid complex declaration of input \ parameter `%s' ignored!\n" p.Parameter.name end; (p :: input, derived)) ([], []) (filter_constants (values model.parameters)) in let additional = missing_parameters input derived model.couplings in (List.sort compare_parameters input @ List.map Parameter.missing additional, List.sort compare_parameters derived) (*i List.iter (fun line -> Printf.eprintf "par: %s\n" line) (dependencies_to_strings derived_dependencies); List.iter (fun line -> Printf.eprintf "coupling: %s\n" line) (dependencies_to_strings coupling_dependencies); i*) let translate_name map name = try SMap.find name map with Not_found -> name let translate_input map p = (translate_name map p.Parameter.name, value_to_float p.Parameter.value) let alpha_s_half e = UFOx.Expr.substitute "aS" (UFOx.Expr.half "aS") e let alpha_s_half_etc map e = UFOx.Expr.rename (map_to_alist map) (alpha_s_half e) let translate_derived map p = let make_atom s = s in let c = make_atom (translate_name map p.Parameter.name) and v = value_to_coupling (alpha_s_half_etc map) make_atom p.Parameter.value in match p.Parameter.ptype with | Parameter.Real -> (Coupling.Real c, v) | Parameter.Complex -> (Coupling.Complex c, v) let translate_coupling_constant map c = let make_atom s = s in (Coupling.Complex c.UFO_Coupling.name, Coupling.Quot (value_to_coupling (alpha_s_half_etc map) make_atom (Expr c.UFO_Coupling.value), Coupling.I)) module Lowercase_Parameters = struct type elt = string type base = string let compare_elt = compare let compare_base = compare let pi = ThoString.lowercase end module Lowercase_Bundle = Bundle.Make (Lowercase_Parameters) let coupling_names model = SMap.fold (fun _ c acc -> c.UFO_Coupling.name :: acc) model.couplings [] let parameter_names model = SMap.fold (fun _ c acc -> c.Parameter.name :: acc) model.parameters [] let ambiguous_parameters model = let all_names = List.rev_append (coupling_names model) (parameter_names model) in let lc_bundle = Lowercase_Bundle.of_list all_names in let lc_set = List.fold_left (fun acc s -> SSet.add s acc) SSet.empty (Lowercase_Bundle.base lc_bundle) and ambiguities = List.filter (fun (_, names) -> List.length names > 1) (Lowercase_Bundle.fibers lc_bundle) in (lc_set, ambiguities) let disambiguate1 lc_set name = let rec disambiguate1' i = let name' = Printf.sprintf "%s_%d" name i in let lc_name' = ThoString.lowercase name' in if SSet.mem lc_name' lc_set then disambiguate1' (succ i) else (SSet.add lc_name' lc_set, name') in disambiguate1' 1 let disambiguate lc_set names = let _, replacements = List.fold_left (fun (lc_set', acc) name -> let lc_set'', name' = disambiguate1 lc_set' name in (lc_set'', SMap.add name name' acc)) (lc_set, SMap.empty) names in replacements let omegalib_names = ["u"; "ubar"; "v"; "vbar"; "eps"] let translate_parameters model = let lc_set, ambiguities = ambiguous_parameters model in let replacements = disambiguate lc_set (ThoList.flatmap snd ambiguities) in SMap.iter (Printf.eprintf "warning: case sensitive parameter names: renaming '%s' -> '%s'\n") replacements; let replacements = List.fold_left (fun acc name -> SMap.add name ("UFO_" ^ name) acc) replacements omegalib_names in let input_parameters, derived_parameters = classify_parameters model and couplings = values model.couplings in { Coupling.input = List.map (translate_input replacements) input_parameters; Coupling.derived = List.map (translate_derived replacements) derived_parameters @ List.map (translate_coupling_constant replacements) couplings; Coupling.derived_arrays = [] } (* UFO requires us to look up the mass parameter to distinguish between massless and massive vectors. TODO: this is a candidate for another lookup table. *) let lorentz_of_particle p = match UFOx.Lorentz.omega p.Particle.spin with | Coupling.Vector -> begin match ThoString.uppercase p.Particle.mass with | "ZERO" -> Coupling.Vector | _ -> Coupling.Massive_Vector end | s -> s type state = { directory : string; model : t } let initialized = ref None let is_initialized_from dir = match !initialized with | None -> false | Some state -> dir = state.directory let dump_raw = ref false let init dir = let model = filter_unphysical (parse_directory dir) in if !dump_raw then dump model; let tables = Lookup.of_model model in let vertices () = translate_vertices model tables in let particle f = tables.Lookup.particle f in let lorentz f = lorentz_of_particle (particle f) in let propagator f = let p = particle f in match p.Particle.propagator with | None -> propagator_of_lorentz (lorentz_of_particle p) | Some s -> Coupling.Prop_UFO s in let gauge_symbol () = "?GAUGE?" in let constant_symbol s = s in let parameters = translate_parameters model in M.setup ~color:(fun f -> UFOx.Color.omega (particle f).Particle.color) ~nc:(fun () -> model.nc) ~pdg:(fun f -> (particle f).Particle.pdg_code) ~lorentz ~propagator ~width:(fun f -> Coupling.Constant) ~goldstone:(fun f -> None) ~conjugate:tables.Lookup.conjugate ~fermion:(fun f -> fermion_of_lorentz (lorentz f)) ~vertices ~flavors:[("All Flavors", tables.Lookup.flavors)] ~parameters:(fun () -> parameters) ~flavor_of_string:tables.Lookup.flavor_of_string ~flavor_to_string:(fun f -> (particle f).Particle.name) ~flavor_to_TeX:(fun f -> (particle f).Particle.texname) ~flavor_symbol:tables.Lookup.flavor_symbol ~gauge_symbol ~mass_symbol:(fun f -> (particle f).Particle.mass) ~width_symbol:(fun f -> (particle f).Particle.width) ~constant_symbol; initialized := Some { directory = dir; model = model } let ufo_directory = ref Config.default_UFO_dir let load () = if is_initialized_from !ufo_directory then () else init !ufo_directory let include_all_fusions = ref false (* In case of Majorana spinors, also generate all combinations of charge conjugated fermion lines. The naming convention is to append \texttt{\_c}$nm$ if the $\gamma$-matrices of the fermion line $n\to m$ has been charge conjugated (this could become impractical for too many fermions at a vertex, but shouldn't matter in real life). *) - (* \begin{dubious} - To be decided: is it better (in the sense of - \emph{easier to understand}, not efficiency) - to permute first or to charge conjugate first? - \end{dubious} *) - (* Here we alway generate \emph{all} charge conjugations, because we treat \emph{all} fermions as Majorana fermion, if there is at least one Majorana fermion in the model! *) let is_majorana = function | Coupling.Majorana | Coupling.Vectorspinor | Coupling.Maj_Ghost -> true | _ -> false let name_spins_structure spins l = (l.Lorentz.name, spins, l.Lorentz.structure) let fusions_of_model ?only model = let include_fusion = match !include_all_fusions, only with | true, _ | false, None -> (fun name -> true) | false, Some names -> (fun name -> SSet.mem name names) in SMap.fold (fun name l acc -> if include_fusion name then List.fold_left (fun acc p -> let l' = Lorentz.permute p l in match l'.Lorentz.spins with | Lorentz.Unused -> acc | Lorentz.Unique spins -> if ThoArray.exists is_majorana spins then List.map (name_spins_structure spins) (Lorentz.required_charge_conjugates l') @ acc else name_spins_structure spins l' :: acc | Lorentz.Ambiguous _ -> failwith "fusions: Lorentz.Ambiguous") [] (Permutation.Default.cyclic l.Lorentz.n) @ acc else acc) model.lorentz [] let fusions ?only () = match !initialized with | None -> [] | Some { model = model } -> fusions_of_model ?only model let propagators_of_model ?only model = let include_propagator = match !include_all_fusions, only with | true, _ | false, None -> (fun name -> true) | false, Some names -> (fun name -> SSet.mem name names) in SMap.fold (fun name p acc -> if include_propagator name then (name, p) :: acc else acc) model.propagators [] let propagators ?only () = match !initialized with | None -> [] | Some { model = model } -> propagators_of_model ?only model let include_hadrons = ref true let ufo_majorana_warnings = [ "***************************************************"; "* *"; "* CAVEAT: *"; "* *"; "* These amplitudes have been computed for a *"; "* UFO model containing Majorana fermions. *"; "* This version of O'Mega contains some known *"; "* bugs for this case. It was released early at *"; "* the request of the Linear Collider community. *"; "* *"; "* These amplitudes MUST NOT be used for *"; "* publications without prior consulation *"; "* with the WHIZARD authors !!! *"; "* *"; "***************************************************" ] let caveats () = if !use_majorana_spinors then ufo_majorana_warnings else [] module Whizard : sig val write : unit -> unit end = struct let write_header dir = Printf.printf "# WHIZARD Model file derived from UFO directory\n"; Printf.printf "# '%s'\n\n" dir; List.iter (fun s -> Printf.printf "# %s\n" s) (M.caveats ()); Printf.printf "model \"%s\"\n\n" (Filename.basename dir) let write_input_parameters parameters = let open Parameter in Printf.printf "# Independent (input) Parameters\n"; List.iter (fun p -> Printf.printf "parameter %s = %s" p.name (value_to_numeric p.value); begin match p.lhablock, p.lhacode with | None, None -> () | Some name, Some (index :: indices) -> Printf.printf " slha_entry %s %d" name index; List.iter (fun i -> Printf.printf " %d" i) indices | Some name, None -> Printf.eprintf "UFO: parameter %s: slhablock %s without slhacode\n" p.name name | Some name, Some [] -> Printf.eprintf "UFO: parameter %s: slhablock %s with empty slhacode\n" p.name name | None, Some _ -> Printf.eprintf "UFO: parameter %s: slhacode without slhablock\n" p.name end; Printf.printf "\n") parameters; Printf.printf "\n" let write_derived_parameters parameters = let open Parameter in Printf.printf "# Dependent (derived) Parameters\n"; List.iter (fun p -> Printf.printf "derived %s = %s\n" p.name (value_to_expr alpha_s_half p.value)) parameters let write_particles particles = let open Particle in Printf.printf "# Particles\n"; Printf.printf "# NB: hypercharge assignments appear to be unreliable\n"; Printf.printf "# therefore we can't infer the isospin\n"; Printf.printf "# NB: parton-, gauge- & handedness are unavailable\n"; List.iter (fun p -> if not p.is_anti then begin Printf.printf "particle \"%s\" %d ### parton? gauge? left?\n" p.name p.pdg_code; Printf.printf " spin %s charge %s color %s ### isospin?\n" (UFOx.Lorentz.rep_to_string_whizard p.spin) (charge_to_string p.charge) (UFOx.Color.rep_to_string_whizard p.color); Printf.printf " name \"%s\"\n" p.name; if p.antiname <> p.name then Printf.printf " anti \"%s\"\n" p.antiname; Printf.printf " tex_name \"%s\"\n" p.texname; if p.antiname <> p.name then Printf.printf " tex_anti \"%s\"\n" p.antitexname; Printf.printf " mass %s width %s\n\n" p.mass p.width end) (values particles); Printf.printf "\n" let write_hadrons () = Printf.printf "# Hadrons (protons and beam remnants)\n"; Printf.printf "# NB: these are NOT part of the UFO model\n"; Printf.printf "# but added for WHIZARD's convenience!\n"; Printf.printf "particle PROTON 2212\n"; Printf.printf " spin 1/2 charge 1\n"; Printf.printf " name p \"p+\"\n"; Printf.printf " anti pbar \"p-\"\n"; Printf.printf "particle HADRON_REMNANT 90\n"; Printf.printf " name hr\n"; Printf.printf " tex_name \"had_r\"\n"; Printf.printf "particle HADRON_REMNANT_SINGLET 91\n"; Printf.printf " name hr1\n"; Printf.printf " tex_name \"had_r^{(1)}\"\n"; Printf.printf "particle HADRON_REMNANT_TRIPLET 92\n"; Printf.printf " color 3\n"; Printf.printf " name hr3\n"; Printf.printf " tex_name \"had_r^{(3)}\"\n"; Printf.printf " anti hr3bar\n"; Printf.printf " tex_anti \"had_r^{(\\bar 3)}\"\n"; Printf.printf "particle HADRON_REMNANT_OCTET 93\n"; Printf.printf " color 8\n"; Printf.printf " name hr8\n"; Printf.printf " tex_name \"had_r^{(8)}\"\n"; Printf.printf "\n" let write_vertices model vertices = Printf.printf "# Vertices (for phasespace generation only)\n"; Printf.printf "# NB: particles should be sorted increasing in mass.\n"; Printf.printf "# This is NOT implemented yet!\n"; List.iter (fun v -> let particles = String.concat " " (List.map (fun s -> "\"" ^ (SMap.find s model.particles).Particle.name ^ "\"") (Array.to_list v.Vertex.particles)) in Printf.printf "vertex %s\n" particles) (values vertices); Printf.printf "\n" let write () = match !initialized with | None -> failwith "UFO.Whizard.write: UFO model not initialized" | Some { directory = dir; model = model } -> let input_parameters, derived_parameters = classify_parameters model in write_header dir; write_input_parameters input_parameters; write_derived_parameters derived_parameters; write_particles model.particles; if !include_hadrons then write_hadrons (); write_vertices model model.vertices; exit 0 end let options = Options.create [ ("UFO_dir", Arg.String (fun name -> ufo_directory := name), "UFO model directory (default: " ^ !ufo_directory ^ ")"); ("Majorana", Arg.Set use_majorana_spinors, "use Majorana spinors (must come _before_ exec!)"); ("divide_propagators_by_i", Arg.Set divide_propagators_by_i, "divide propagators by I (pre 2013 FeynRules convention)"); ("verbatim_Hg", Arg.Set verbatim_higgs_glue, "don't correct the color flows for effective Higgs Gluon couplings"); ("write_WHIZARD", Arg.Unit Whizard.write, "write the WHIZARD model file (required once per model)"); ("long_flavors", Arg.Unit (fun () -> Lookup.flavor_format := Lookup.Long), "write use the UFO flavor names instead of integers"); ("dump", Arg.Set dump_raw, "dump UFO model for debugging the parser (must come _before_ exec!)"); ("all_fusions", Arg.Set include_all_fusions, "include all fusions in the fortran module"); ("no_hadrons", Arg.Clear include_hadrons, "don't add any particle not in the UFO file"); ("add_hadrons", Arg.Set include_hadrons, "add protons and beam remants for WHIZARD"); ("exec", Arg.Unit load, "load the UFO model files (required _before_ using particles names)"); ("help", Arg.Unit (fun () -> prerr_endline "..."), "print information on the model")] end module type Fortran_Target = sig val fuse : Algebra.QC.t -> string -> Coupling.lorentzn -> Coupling.fermion_lines -> string -> string list -> string list -> Coupling.fusen -> unit val lorentz_module : ?only:SSet.t -> ?name:string -> ?fortran_module:string -> ?parameter_module:string -> Format_Fortran.formatter -> unit -> unit end module Targets = struct module Fortran : Fortran_Target = struct open Format_Fortran let fuse = UFO_targets.Fortran.fuse let lorentz_functions ff fusions () = List.iter (fun (name, s, l) -> UFO_targets.Fortran.lorentz ff name s l) fusions let propagator_functions ff parameter_module propagators () = List.iter (fun (name, p) -> UFO_targets.Fortran.propagator ff name parameter_module p.Propagator.variables p.Propagator.spins p.Propagator.numerator p.Propagator.denominator) propagators let lorentz_module ?only ?(name="omega_amplitude_ufo") ?(fortran_module="omega95") ?(parameter_module="parameter_module") ff () = let printf fmt = fprintf ff fmt and nl = pp_newline ff in printf "module %s" name; nl (); printf " use kinds"; nl (); printf " use %s" fortran_module; nl (); printf " implicit none"; nl (); printf " private"; nl (); let fusions = Model.fusions ?only () and propagators = Model.propagators () in List.iter (fun (name, _, _) -> printf " public :: %s" name; nl ()) fusions; List.iter (fun (name, _) -> printf " public :: pr_U_%s" name; nl ()) propagators; UFO_targets.Fortran.eps4_g4_g44_decl ff (); UFO_targets.Fortran.eps4_g4_g44_init ff (); printf "contains"; nl (); UFO_targets.Fortran.inner_product_functions ff (); lorentz_functions ff fusions (); propagator_functions ff parameter_module propagators (); printf "end module %s" name; nl (); pp_flush ff () end end module type Test = sig val suite : OUnit.test end module Test : Test = struct open OUnit let lexer s = UFO_lexer.token (UFO_lexer.init_position "" (Lexing.from_string s)) let suite_lexer_escapes = "escapes" >::: [ "single-quote" >:: (fun () -> assert_equal (UFO_parser.STRING "a'b'c") (lexer "'a\\'b\\'c'")); "unterminated" >:: (fun () -> assert_raises End_of_file (fun () -> lexer "'a\\'b\\'c")) ] let suite_lexer = "lexer" >::: [suite_lexer_escapes] let suite = "UFO" >::: [suite_lexer] end Index: trunk/omega/src/UFO_targets.ml =================================================================== --- trunk/omega/src/UFO_targets.ml (revision 8494) +++ trunk/omega/src/UFO_targets.ml (revision 8495) @@ -1,1551 +1,1545 @@ (* UFO_targets.ml -- Copyright (C) 1999-2017 by Wolfgang Kilian Thorsten Ohl Juergen Reuter with contributions from Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) let (@@) f g x = f (g x) (* \thocwmodulesection{Generating Code for UFO Lorentz Structures} *) (* O'Caml before 4.02 had a module typing bug that forces us to put this definition outside [Lorentz_Fusion]. *) module Q = Algebra.Q module QC = Algebra.QC module type T = sig (* [lorentz formatter name spins v] writes a representation of the Lorentz structure [v] of particles with the Lorentz representations [spins] as a (Fortran) function [name] to [formatter]. *) val lorentz : Format_Fortran.formatter -> string -> Coupling.lorentz array -> UFO_Lorentz.t -> unit val propagator : Format_Fortran.formatter -> string -> string -> string list -> Coupling.lorentz * Coupling.lorentz -> UFO_Lorentz.t -> UFO_Lorentz.t -> unit val fusion_name : string -> Permutation.Default.t -> Coupling.fermion_lines -> string val fuse : Algebra.QC.t -> string -> Coupling.lorentzn -> Coupling.fermion_lines -> string -> string list -> string list -> Coupling.fusen -> unit val eps4_g4_g44_decl : Format_Fortran.formatter -> unit -> unit val eps4_g4_g44_init : Format_Fortran.formatter -> unit -> unit val inner_product_functions : Format_Fortran.formatter -> unit -> unit module type Test = sig val suite : OUnit.test end module Test : Test end module Fortran : T = struct open Format_Fortran let pp_divide ?(indent=0) ff () = fprintf ff "%*s! %s" indent "" (String.make (70 - indent) '-'); pp_newline ff () let conjugate = function | Coupling.Spinor -> Coupling.ConjSpinor | Coupling.ConjSpinor -> Coupling.Spinor | r -> r let spin_mnemonic = function | Coupling.Scalar -> "phi" | Coupling.Spinor -> "psi" | Coupling.ConjSpinor -> "psibar" | Coupling.Majorana -> "chi" | Coupling.Maj_Ghost -> invalid_arg "UFO_targets: Maj_Ghost" | Coupling.Vector -> "a" | Coupling.Massive_Vector -> "v" | Coupling.Vectorspinor -> "grav" (* itino *) | Coupling.Tensor_1 -> invalid_arg "UFO_targets: Tensor_1" | Coupling.Tensor_2 -> "h" | Coupling.BRS l -> invalid_arg "UFO_targets: BRS" let fortran_type = function | Coupling.Scalar -> "complex(kind=default)" | Coupling.Spinor -> "type(spinor)" | Coupling.ConjSpinor -> "type(conjspinor)" | Coupling.Majorana -> "type(bispinor)" | Coupling.Maj_Ghost -> invalid_arg "UFO_targets: Maj_Ghost" | Coupling.Vector -> "type(vector)" | Coupling.Massive_Vector -> "type(vector)" | Coupling.Vectorspinor -> "type(vectorspinor)" | Coupling.Tensor_1 -> invalid_arg "UFO_targets: Tensor_1" | Coupling.Tensor_2 -> "type(tensor)" | Coupling.BRS l -> invalid_arg "UFO_targets: BRS" (* The \texttt{omegalib} separates time from space. Maybe not a good idea after all. Mend it locally \ldots *) type wf = { pos : int; spin : Coupling.lorentz; name : string; local_array : string option; momentum : string; momentum_array : string; fortran_type : string } let wf_table spins = Array.mapi (fun i s -> let spin = if i = 0 then conjugate s else s in let pos = succ i in let i = string_of_int pos in let name = spin_mnemonic s ^ i in let local_array = begin match spin with | Coupling.Vector | Coupling.Massive_Vector -> Some (name ^ "a") | _ -> None end in { pos; spin; name; local_array; momentum = "k" ^ i; momentum_array = "p" ^ i; fortran_type = fortran_type spin } ) spins module L = UFO_Lorentz (* Format rational ([Q.t]) and complex rational ([QC.t]) numbers as fortran values. *) let format_rational q = if Q.is_integer q then string_of_int (Q.to_integer q) else let n, d = Q.to_ratio q in Printf.sprintf "%d.0_default/%d" n d let format_complex_rational cq = let real = QC.real cq and imag = QC.imag cq in if Q.is_null imag then begin if Q.is_negative real then "(" ^ format_rational real ^ ")" else format_rational real end else if Q.is_integer real && Q.is_integer imag then Printf.sprintf "(%d,%d)" (Q.to_integer real) (Q.to_integer imag) else Printf.sprintf "cmplx(%s,%s,kind=default)" (format_rational real) (format_rational imag) (* Optimize the representation if used as a prefactor of a summand in a sum. *) let format_rational_factor q = if Q.is_unit q then "+ " else if Q.is_unit (Q.neg q) then "- " else if Q.is_negative q then "- " ^ format_rational (Q.neg q) ^ "*" else "+ " ^ format_rational q ^ "*" let format_complex_rational_factor cq = let real = QC.real cq and imag = QC.imag cq in if Q.is_null imag then begin if Q.is_unit real then "+ " else if Q.is_unit (Q.neg real) then "- " else if Q.is_negative real then "- " ^ format_rational (Q.neg real) ^ "*" else "+ " ^ format_rational real ^ "*" end else if Q.is_integer real && Q.is_integer imag then Printf.sprintf "+ (%d,%d)*" (Q.to_integer real) (Q.to_integer imag) else Printf.sprintf "+ cmplx(%s,%s,kind=default)*" (format_rational real) (format_rational imag) (* Append a formatted list of indices to [name]. *) let append_indices name = function | [] -> name | indices -> name ^ "(" ^ String.concat "," (List.map string_of_int indices) ^ ")" (* Dirac string variables and their names. *) type dsv = | Ket of int | Bra of int | Braket of int let dsv_name = function | Ket n -> Printf.sprintf "ket%02d" n | Bra n -> Printf.sprintf "bra%02d" n | Braket n -> Printf.sprintf "bkt%02d" n let dirac_dimension dsv indices = let tail ilist = String.concat "," (List.map (fun _ -> "0:3") ilist) ^ ")" in match dsv, indices with | Braket _, [] -> "" | (Ket _ | Bra _), [] -> ", dimension(1:4)" | Braket _, indices -> ", dimension(" ^ tail indices | (Ket _ | Bra _), indices -> ", dimension(1:4," ^ tail indices (* Write Fortran code to [decl] and [eval]: apply the Dirac matrix [gamma] with complex rational entries to the spinor [ket] from the left. [ket] must be the name of a scalar variable and cannot be an array element. The result is stored in [dsv_name (Ket n)] which can have additional [indices]. Return [Ket n] for further processing. *) let dirac_ket_to_fortran_decl ff n indices = let printf fmt = fprintf ff fmt and nl = pp_newline ff in let dsv = Ket n in printf " @[<2>complex(kind=default)%s ::@ %s@]" (dirac_dimension dsv indices) (dsv_name dsv); nl () let dirac_ket_to_fortran_eval ff n indices gamma ket = let printf fmt = fprintf ff fmt and nl = pp_newline ff in let dsv = Ket n in for i = 0 to 3 do let name = append_indices (dsv_name dsv) (succ i :: indices) in printf " @[<%d>%s = 0" (String.length name + 4) name; for j = 0 to 3 do if not (QC.is_null gamma.(i).(j)) then printf "@ %s%s%%a(%d)" (format_complex_rational_factor gamma.(i).(j)) ket.name (succ j) done; printf "@]"; nl () done; dsv (* The same as [dirac_ket_to_fortran], but apply the Dirac matrix [gamma] to [bra] from the right and return [Bra n]. *) let dirac_bra_to_fortran_decl ff n indices = let printf fmt = fprintf ff fmt and nl = pp_newline ff in let dsv = Bra n in printf " @[<2>complex(kind=default)%s ::@ %s@]" (dirac_dimension dsv indices) (dsv_name dsv); nl () let dirac_bra_to_fortran_eval ff n indices bra gamma = let printf fmt = fprintf ff fmt and nl = pp_newline ff in let dsv = Bra n in for j = 0 to 3 do let name = append_indices (dsv_name dsv) (succ j :: indices) in printf " @[<%d>%s = 0" (String.length name + 4) name; for i = 0 to 3 do if not (QC.is_null gamma.(i).(j)) then printf "@ %s%s%%a(%d)" (format_complex_rational_factor gamma.(i).(j)) bra.name (succ i) done; printf "@]"; nl () done; dsv (* More of the same, but evaluating a spinor sandwich and returning [Braket n]. *) let dirac_braket_to_fortran_decl ff n indices = let printf fmt = fprintf ff fmt and nl = pp_newline ff in let dsv = Braket n in printf " @[<2>complex(kind=default)%s ::@ %s@]" (dirac_dimension dsv indices) (dsv_name dsv); nl () let dirac_braket_to_fortran_eval ff n indices bra gamma ket = let printf fmt = fprintf ff fmt and nl = pp_newline ff in let dsv = Braket n in let name = append_indices (dsv_name dsv) indices in printf " @[<%d>%s = 0" (String.length name + 4) name; for i = 0 to 3 do for j = 0 to 3 do if not (QC.is_null gamma.(i).(j)) then printf "@ %s%s%%a(%d)*%s%%a(%d)" (format_complex_rational_factor gamma.(i).(j)) bra.name (succ i) ket.name (succ j) done done; printf "@]"; nl (); dsv (* Choose among the previous functions according to the position of [bra] and [ket] among the wavefunctions. If any is in the first position evaluate the spinor expression with the corresponding spinor removed, otherwise evaluate the spinir sandwich. *) let dirac_bra_or_ket_to_fortran_decl ff n indices bra ket = if bra = 1 then dirac_ket_to_fortran_decl ff n indices else if ket = 1 then dirac_bra_to_fortran_decl ff n indices else dirac_braket_to_fortran_decl ff n indices let dirac_bra_or_ket_to_fortran_eval ff n indices wfs bra gamma ket = if bra = 1 then dirac_ket_to_fortran_eval ff n indices gamma wfs.(pred ket) else if ket = 1 then dirac_bra_to_fortran_eval ff n indices wfs.(pred bra) gamma else dirac_braket_to_fortran_eval ff n indices wfs.(pred bra) gamma wfs.(pred ket) (* UFO summation indices are negative integers. Derive a valid Fortran variable name. *) let prefix_summation = "mu" let prefix_polarization = "nu" let index_spinor = "alpha" let index_tensor = "nu" let index_variable mu = if mu < 0 then Printf.sprintf "%s%d" prefix_summation (- mu) else if mu == 0 then prefix_polarization else Printf.sprintf "%s%d" prefix_polarization mu let format_indices indices = String.concat "," (List.map index_variable indices) module IntPM = Partial.Make (struct type t = int let compare = compare end) type tensor = | DS of dsv | V of string | T of UFOx.Lorentz_Atom.vector | S of UFOx.Lorentz_Atom.scalar | Inv of UFOx.Lorentz_Atom.scalar (* Transform the Dirac strings if we have Majorana fermions involved, in order to implement the algorithm from JRR's thesis. NB: The following is for reference only, to better understand what JRR was doing\ldots *) (* If the vertex is (suppressing the Lorentz indices of~$\phi_2$ and~$\Gamma$) \begin{equation} \label{eq:FVF-Vertex} \bar\psi \Gamma\phi \psi = \Gamma_{\alpha\beta} \bar\psi_{\alpha} \phi \psi_{\beta} \end{equation} (cf.~[Coupling.FBF] in the hardcoded O'Mega models), then this is the version implemented by [fuse] below. *) let tho_print_dirac_current f c wf1 wf2 fusion = match fusion with | [1; 3] -> printf "%s_ff(%s,%s,%s)" f c wf1 wf2 (* $\Gamma_{\alpha\beta} \bar\psi_{1,\alpha} \psi_{2,\beta}$ *) | [3; 1] -> printf "%s_ff(%s,%s,%s)" f c wf2 wf1 (* $\Gamma_{\alpha\beta} \bar\psi_{1,\alpha} \psi_{2,\beta}$ *) | [2; 3] -> printf "f_%sf(%s,%s,%s)" f c wf1 wf2 (* $\Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta}$ *) | [3; 2] -> printf "f_%sf(%s,%s,%s)" f c wf2 wf1 (* $\Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta}$ *) | [1; 2] -> printf "f_f%s(%s,%s,%s)" f c wf1 wf2 (* $\Gamma_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2$ *) | [2; 1] -> printf "f_f%s(%s,%s,%s)" f c wf2 wf1 (* $\Gamma_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2$ *) | _ -> () (* The corresponding UFO [fuse] exchanges the arguments in the case of two fermions. This is the natural choice for cyclic permutations. *) let tho_print_FBF_current f c wf1 wf2 fusion = match fusion with | [3; 1] -> printf "f%sf_p120(%s,%s,%s)" f c wf1 wf2 (* $\Gamma_{\alpha\beta} \psi_{1,\beta} \bar\psi_{2,\alpha}$ *) | [1; 3] -> printf "f%sf_p120(%s,%s,%s)" f c wf2 wf1 (* $\Gamma_{\alpha\beta} \psi_{1,\beta} \bar\psi_{2,\alpha}$ *) | [2; 3] -> printf "f%sf_p012(%s,%s,%s)" f c wf1 wf2 (* $\Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta}$ *) | [3; 2] -> printf "f%sf_p012(%s,%s,%s)" f c wf2 wf1 (* $\Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta}$ *) | [1; 2] -> printf "f%sf_p201(%s,%s,%s)" f c wf1 wf2 (* $\Gamma_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2$ *) | [2; 1] -> printf "f%sf_p201(%s,%s,%s)" f c wf2 wf1 (* $\Gamma_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2$ *) | _ -> () (* This is how JRR implemented (see subsection~\ref{sec:dirac-matrices-jrr}) the Dirac matrices that don't change sign under $C\Gamma^T C^{-1} = \Gamma$, i.\,e.~$\mathbf{1}$, $\gamma_5$ and~$\gamma_5\gamma_\mu$ (see [Targets.Fortran_Majorana_Fermions.print_fermion_current]) \begin{itemize} \item In the case of two fermions, the second wave function [wf2] is always put into the second slot, as described in JRR's thesis. \label{pg:JRR-Fusions} \item In the case of a boson and a fermion, there is no need for both ["f_%sf"] and ["f_f%s"], since the latter can be obtained by exchanging arguments. \end{itemize} *) let jrr_print_majorana_current_S_P_A f c wf1 wf2 fusion = match fusion with | [1; 3] -> printf "%s_ff(%s,%s,%s)" f c wf1 wf2 (* $(C\Gamma)_{\alpha\beta} \bar\psi_{1,\alpha} \psi_{2,\beta} \cong C\Gamma $ *) | [3; 1] -> printf "%s_ff(%s,%s,%s)" f c wf1 wf2 (* $(C\Gamma)_{\alpha\beta} \psi_{1,\alpha} \bar\psi_{2,\beta} \cong C\Gamma = C\,C\Gamma^T C^{-1} $ *) | [2; 3] -> printf "f_%sf(%s,%s,%s)" f c wf1 wf2 (* $\Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [3; 2] -> printf "f_%sf(%s,%s,%s)" f c wf2 wf1 (* $\Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [1; 2] -> printf "f_%sf(%s,%s,%s)" f c wf2 wf1 (* $\Gamma_{\alpha\beta} \phi_1 \bar\psi_{2,\beta} \cong \Gamma = C\Gamma^T C^{-1} $ *) | [2; 1] -> printf "f_%sf(%s,%s,%s)" f c wf1 wf2 (* $\Gamma_{\alpha\beta} \phi_1 \bar\psi_{2,\beta} \cong \Gamma = C\Gamma^T C^{-1} $ *) | _ -> () (* This is how JRR implemented the Dirac matrices that do change sign under $C\Gamma^T C^{-1} = - \Gamma$, i.\,e.~$\gamma_\mu$ and~$\sigma_{\mu\nu}$ (see [Targets.Fortran_Majorana_Fermions.print_fermion_current_vector]). *) let jrr_print_majorana_current_V f c wf1 wf2 fusion = match fusion with | [1; 3] -> printf "%s_ff( %s,%s,%s)" f c wf1 wf2 (* $ (C\Gamma)_{\alpha\beta} \bar\psi_{1,\alpha} \psi_{2,\beta} \cong C\Gamma $ *) | [3; 1] -> printf "%s_ff(-%s,%s,%s)" f c wf1 wf2 (* $-(C\Gamma)_{\alpha\beta} \psi_{1,\alpha} \bar\psi_{2,\beta} \cong -C\Gamma = C\,C\Gamma^T C^{-1} $ *) | [2; 3] -> printf "f_%sf( %s,%s,%s)" f c wf1 wf2 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [3; 2] -> printf "f_%sf( %s,%s,%s)" f c wf2 wf1 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [1; 2] -> printf "f_%sf(-%s,%s,%s)" f c wf2 wf1 (* $-\Gamma_{\alpha\beta} \phi_1 \bar\psi_{2,\beta} \cong -\Gamma = C\Gamma^T C^{-1} $ *) | [2; 1] -> printf "f_%sf(-%s,%s,%s)" f c wf1 wf2 (* $-\Gamma_{\alpha\beta} \phi_1 \bar\psi_{2,\beta} \cong -\Gamma = C\Gamma^T C^{-1} $ *) | _ -> () (* These two can be unified, if the \texttt{\_c} functions implement~$\Gamma'=C\Gamma^T C^{-1}$, but we \emph{must} make sure that the multiplication with~$C$ from the left happens \emph{after} the transformation~$\Gamma\to\Gamma'$. *) let jrr_print_majorana_current f c wf1 wf2 fusion = match fusion with | [1; 3] -> printf "%s_ff (%s,%s,%s)" f c wf1 wf2 (* $ (C\Gamma)_{\alpha\beta} \bar\psi_{1,\alpha} \psi_{2,\beta} \cong C\Gamma $ *) | [3; 1] -> printf "%s_ff_c(%s,%s,%s)" f c wf1 wf2 (* $(C\Gamma')_{\alpha\beta} \psi_{1,\alpha} \bar\psi_{2,\beta} \cong C\Gamma' = C\,C\Gamma^T C^{-1} $ *) | [2; 3] -> printf "f_%sf (%s,%s,%s)" f c wf1 wf2 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [3; 2] -> printf "f_%sf (%s,%s,%s)" f c wf2 wf1 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [1; 2] -> printf "f_%sf_c(%s,%s,%s)" f c wf2 wf1 (* $\Gamma'_{\alpha\beta} \phi_1 \bar\psi_{2,\beta} \cong \Gamma' = C\Gamma^T C^{-1} $ *) | [2; 1] -> printf "f_%sf_c(%s,%s,%s)" f c wf1 wf2 (* $\Gamma'_{\alpha\beta} \phi_1 \bar\psi_{2,\beta} \cong \Gamma' = C\Gamma^T C^{-1} $ *) | _ -> () (* Since we may assume~$C^{-1}=-C=C^T$, this can be rewritten if the \texttt{\_c} functions implement \begin{equation} \Gamma^{\prime\,T} = \left(C\Gamma^T C^{-1}\right)^T = \left(C^{-1}\right)^T \Gamma C^T = C \Gamma C^{-1} \end{equation} instead. *) let jrr_print_majorana_current_transposing f c wf1 wf2 fusion = match fusion with | [1; 3] -> printf "%s_ff (%s,%s,%s)" f c wf1 wf2 (* $ (C\Gamma)_{\alpha\beta} \bar\psi_{1,\alpha} \psi_{2,\beta} \cong C\Gamma $ *) | [3; 1] -> printf "%s_ff_c(%s,%s,%s)" f c wf2 wf1 (* $(C\Gamma')^T_{\alpha\beta} \bar\psi_{1,\alpha} \psi_{2,\beta} \cong (C\Gamma')^T = - C\Gamma $ *) | [2; 3] -> printf "f_%sf (%s,%s,%s)" f c wf1 wf2 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [3; 2] -> printf "f_%sf (%s,%s,%s)" f c wf2 wf1 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [1; 2] -> printf "f_f%s_c(%s,%s,%s)" f c wf1 wf2 (* $\Gamma^{\prime\,T}_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2 \cong \Gamma^{\prime\,T} = C\Gamma C^{-1}$ *) | [2; 1] -> printf "f_f%s_c(%s,%s,%s)" f c wf2 wf1 (* $\Gamma^{\prime\,T}_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2 \cong \Gamma^{\prime\,T} = C\Gamma C^{-1} $ *) | _ -> () (* where we have used \begin{equation} (C\Gamma')^T = \Gamma^{\prime,T}C^T = C\Gamma C^{-1} C^T = C\Gamma C^{-1} (-C) = - C\Gamma\,. \end{equation} *) (* This puts the arguments in the same slots as [tho_print_dirac_current] above and can be implemented by [fuse], iff we inject the proper transformations in [dennerize] below. We notice that we do \emph{not} need the conjugated version for all combinations, but only for the case of two fermions. In the two cases of one column spinor~$\psi$, only the original version appears and in the two cases of one row spinor~$\bar\psi$, only the conjugated version appears. *) (* Before we continue, we must however generalize from the assumption~\eqref{eq:FVF-Vertex} that the fields in the vertex are always ordered as in~[Coupling.FBF]. First, even in this case the slots of the fermions must be exchanged to accomodate the cyclic permutations. Therefore we exchange the arguments of the [[1; 3]] and [[3; 1]] fusions. *) let jrr_print_majorana_FBF f c wf1 wf2 fusion = match fusion with (* [fline = (3, 1)] *) | [3; 1] -> printf "f%sf_p120_c(%s,%s,%s)" f c wf1 wf2 (* $(C\Gamma')^T_{\alpha\beta} \psi_{1,\beta} \bar\psi_{2,\alpha} \cong (C\Gamma')^T = - C\Gamma $ *) | [1; 3] -> printf "f%sf_p120 (%s,%s,%s)" f c wf2 wf1 (* $ (C\Gamma)_{\alpha\beta} \psi_{1,\beta} \bar\psi_{2,\alpha} \cong C\Gamma $ *) | [2; 3] -> printf "f%sf_p012 (%s,%s,%s)" f c wf1 wf2 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [3; 2] -> printf "f%sf_p012 (%s,%s,%s)" f c wf2 wf1 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [1; 2] -> printf "f%sf_p201 (%s,%s,%s)" f c wf1 wf2 (* $\Gamma^{\prime\,T}_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2 \cong \Gamma^{\prime\,T} = C\Gamma C^{-1}$ *) | [2; 1] -> printf "f%sf_p201 (%s,%s,%s)" f c wf2 wf1 (* $\Gamma^{\prime\,T}_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2 \cong \Gamma^{\prime\,T} = C\Gamma C^{-1} $ *) | _ -> () (* The other two permutations: *) let jrr_print_majorana_FFB f c wf1 wf2 fusion = match fusion with (* [fline = (1, 2)] *) | [3; 1] -> printf "ff%s_p120 (%s,%s,%s)" f c wf1 wf2 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [1; 3] -> printf "ff%s_p120 (%s,%s,%s)" f c wf2 wf1 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [2; 3] -> printf "ff%s_p012 (%s,%s,%s)" f c wf1 wf2 (* $\Gamma^{\prime\,T}_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2 \cong \Gamma^{\prime\,T} = C\Gamma C^{-1}$ *) | [3; 2] -> printf "ff%s_p012 (%s,%s,%s)" f c wf2 wf1 (* $\Gamma^{\prime\,T}_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2 \cong \Gamma^{\prime\,T} = C\Gamma C^{-1} $ *) | [1; 2] -> printf "ff%s_p201 (%s,%s,%s)" f c wf1 wf2 (* $ (C\Gamma)_{\alpha\beta} \psi_{1,\beta} \bar\psi_{2,\alpha} \cong C\Gamma $ *) | [2; 1] -> printf "ff%s_p201_c(%s,%s,%s)" f c wf2 wf1 (* $(C\Gamma')^T_{\alpha\beta} \psi_{1,\beta} \bar\psi_{2,\alpha} \cong (C\Gamma')^T = - C\Gamma $ *) | _ -> () let jrr_print_majorana_BFF f c wf1 wf2 fusion = match fusion with (* [fline = (2, 3)] *) | [3; 1] -> printf "%sff_p120 (%s,%s,%s)" f c wf1 wf2 (* $\Gamma^{\prime\,T}_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2 \cong \Gamma^{\prime\,T} = C\Gamma C^{-1} $ *) | [1; 3] -> printf "%sff_p120 (%s,%s,%s)" f c wf2 wf1 (* $\Gamma^{\prime\,T}_{\alpha\beta} \bar\psi_{1,\alpha} \phi_2 \cong \Gamma^{\prime\,T} = C\Gamma C^{-1}$ *) | [2; 3] -> printf "%sff_p012 (%s,%s,%s)" f c wf1 wf2 (* $ (C\Gamma)_{\alpha\beta} \psi_{1,\beta} \bar\psi_{2,\alpha} \cong C\Gamma $ *) | [3; 2] -> printf "%sff_p012_c(%s,%s,%s)" f c wf2 wf1 (* $(C\Gamma')^T_{\alpha\beta} \psi_{1,\beta} \bar\psi_{2,\alpha} \cong (C\Gamma')^T = - C\Gamma $ *) | [1; 2] -> printf "%sff_p201 (%s,%s,%s)" f c wf1 wf2 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | [2; 1] -> printf "%sff_p201 (%s,%s,%s)" f c wf2 wf1 (* $ \Gamma_{\alpha\beta} \phi_1 \psi_{2,\beta} \cong \Gamma $ *) | _ -> () - (* \begin{dubious} - Now we want to test \emph{only} if the [fermion_line] - matches the inverted [fusion]. Why does [jrr_print_majorana_FBF] - differ from the others? Must be a typo. - \end{dubious} *) - (* In the model, the necessary information is provided as [Coupling.fermion_lines], encoded as [(right,left)] in the usual direction of the lines. E.\,g.~the case of~\eqref{eq:FVF-Vertex} is~[(3,1)]. Equivalent information is available as~[(ket, bra)] in [UFO_Lorentz.dirac_string]. *) let is_majorana = function | Coupling.Majorana | Coupling.Vectorspinor | Coupling.Maj_Ghost -> true | _ -> false let is_dirac = function | Coupling.Spinor | Coupling.ConjSpinor -> true | _ -> false let dennerize ~eval wfs atom = let printf fmt = fprintf eval fmt and nl = pp_newline eval in if is_majorana wfs.(pred atom.L.bra).spin || is_majorana wfs.(pred atom.L.ket).spin then if atom.L.bra = 1 then (* Fusing one or more bosons with a ket like fermion: $\chi \leftarrow \Gamma\chi$. *) (* Don't do anything, as per subsection~\ref{sec:dirac-matrices-jrr}. *) atom else if atom.L.ket = 1 then (* We fuse one or more bosons with a bra like fermion: $\bar\chi \leftarrow \bar\chi\Gamma$. *) (* $\Gamma\to C \Gamma C^{-1}$. *) begin let atom = L.conjugate atom in printf " ! conjugated for Majorana"; nl (); printf " ! %s" (L.dirac_string_to_string atom); nl (); atom end else if not atom.L.conjugated then (* We fuse zero or more bosons with a sandwich of fermions. $\phi \leftarrow \bar\chi\gamma\chi$.*) (* Multiply by~$C$ from the left, as per subsection~\ref{sec:dirac-matrices-jrr}. *) begin let atom = L.cc_times atom in printf " ! multiplied by CC for Majorana"; nl (); printf " ! %s" (L.dirac_string_to_string atom); nl (); atom end else (* Transposed: multiply by~$-C$ from the left. *) begin let atom = L.minus (L.cc_times atom) in printf " ! multiplied by -CC for Majorana"; nl (); printf " ! %s" (L.dirac_string_to_string atom); nl (); atom end else atom (* Write the [i]th Dirac string [ds] as Fortran code to [eval], including a shorthand representation as a comment. Return [ds] with [ds.L.atom] replaced by the dirac string variable, i,\,e.~[DS dsv] annotated with the internal and external indices. In addition write the declaration to [decl]. *) let dirac_string_to_fortran ~decl ~eval i wfs ds = let printf fmt = fprintf eval fmt and nl = pp_newline eval in let bra = ds.L.atom.L.bra and ket = ds.L.atom.L.ket in pp_divide ~indent:4 eval (); printf " ! %s" (L.dirac_string_to_string ds.L.atom); nl (); let atom = dennerize ~eval wfs ds.L.atom in begin match ds.L.indices with | [] -> let gamma = L.dirac_string_to_matrix (fun _ -> 0) atom in dirac_bra_or_ket_to_fortran_decl decl i [] bra ket; let dsv = dirac_bra_or_ket_to_fortran_eval eval i [] wfs bra gamma ket in L.map_atom (fun _ -> DS dsv) ds | indices -> dirac_bra_or_ket_to_fortran_decl decl i indices bra ket; let combinations = Product.power (List.length indices) [0; 1; 2; 3] in let dsv = List.map (fun combination -> let substitution = IntPM.of_lists indices combination in let substitute = IntPM.apply substitution in let indices = List.map substitute indices in let gamma = L.dirac_string_to_matrix substitute atom in dirac_bra_or_ket_to_fortran_eval eval i indices wfs bra gamma ket) combinations in begin match ThoList.uniq (List.sort compare dsv) with | [dsv] -> L.map_atom (fun _ -> DS dsv) ds | _ -> failwith "dirac_string_to_fortran: impossible" end end (* Write the Dirac strings in the list [ds_list] as Fortran code to [eval], including shorthand representations as comments. Return the list of variables and corresponding indices to be contracted. *) let dirac_strings_to_fortran ~decl ~eval wfs last ds_list = List.fold_left (fun (i, acc) ds -> let i = succ i in (i, dirac_string_to_fortran ~decl ~eval i wfs ds :: acc)) (last, []) ds_list (* Perform a nested sum of terms, as printed by [print_term] (which takes the number of spaces to indent as only argument) of the cartesian product of [indices] running from 0 to 3. *) let nested_sums ~decl ~eval initial_indent indices print_term = let rec nested_sums' indent = function | [] -> print_term indent | index :: indices -> let var = index_variable index in fprintf eval "%*s@[<2>do %s = 0, 3@]" indent "" var; pp_newline eval (); nested_sums' (indent + 2) indices; pp_newline eval (); fprintf eval "%*s@[<2>end do@]" indent "" in nested_sums' (initial_indent + 2) indices (* Polarization indices also need to be summed over, but they appear only once. *) let indices_of_contractions contractions = let index_pairs, polarizations = L.classify_indices (ThoList.flatmap (fun ds -> ds.L.indices) contractions) in try ThoList.pairs index_pairs @ ThoList.uniq (List.sort compare polarizations) with | Invalid_argument s -> invalid_arg ("indices_of_contractions: " ^ ThoList.to_string string_of_int index_pairs) (*i Printf.eprintf "indices_of_contractions: %s / %s\n" (ThoList.to_string string_of_int index_pairs) (ThoList.to_string string_of_int polarizations); i*) let format_dsv dsv indices = match dsv, indices with | Braket _, [] -> dsv_name dsv | Braket _, ilist -> Printf.sprintf "%s(%s)" (dsv_name dsv) (format_indices indices) | (Bra _ | Ket _), [] -> Printf.sprintf "%s(%s)" (dsv_name dsv) index_spinor | (Bra _ | Ket _), ilist -> Printf.sprintf "%s(%s,%s)" (dsv_name dsv) index_spinor (format_indices indices) let denominator_name = "denom_" let mass_name = "m_" let width_name = "w_" let format_tensor t = let indices = t.L.indices in match t.L.atom with | DS dsv -> format_dsv dsv indices | V vector -> Printf.sprintf "%s(%s)" vector (format_indices indices) | T UFOx.Lorentz_Atom.P (mu, n) -> Printf.sprintf "p%d(%s)" n (index_variable mu) | T UFOx.Lorentz_Atom.Epsilon (mu1, mu2, mu3, mu4) -> Printf.sprintf "eps4_(%s)" (format_indices [mu1; mu2; mu3; mu4]) | T UFOx.Lorentz_Atom.Metric (mu1, mu2) -> if mu1 > 0 && mu2 > 0 then Printf.sprintf "g44_(%s)" (format_indices [mu1; mu2]) else failwith "format_tensor: compress_metrics has failed!" | S (UFOx.Lorentz_Atom.Mass _) -> mass_name | S (UFOx.Lorentz_Atom.Width _) -> width_name | S (UFOx.Lorentz_Atom.P2 i) -> Printf.sprintf "g2_(p%d)" i | S (UFOx.Lorentz_Atom.P12 (i, j)) -> Printf.sprintf "g12_(p%d,p%d)" i j | Inv (UFOx.Lorentz_Atom.Mass _) -> "1/" ^ mass_name | Inv (UFOx.Lorentz_Atom.Width _) -> "1/" ^ width_name | Inv (UFOx.Lorentz_Atom.P2 i) -> Printf.sprintf "1/g2_(p%d)" i | Inv (UFOx.Lorentz_Atom.P12 (i, j)) -> Printf.sprintf "1/g12_(p%d,p%d)" i j | S (UFOx.Lorentz_Atom.Variable s) -> s | Inv (UFOx.Lorentz_Atom.Variable s) -> "1/" ^ s | S (UFOx.Lorentz_Atom.Coeff c) -> UFOx.Value.to_string c | Inv (UFOx.Lorentz_Atom.Coeff c) -> "1/(" ^ UFOx.Value.to_string c ^ ")" let rec multiply_tensors ~decl ~eval = function | [] -> fprintf eval "1"; | [t] -> fprintf eval "%s" (format_tensor t) | t :: tensors -> fprintf eval "%s@,*" (format_tensor t); multiply_tensors ~decl ~eval tensors let pseudo_wfs_for_denominator = Array.init 2 (fun i -> let ii = string_of_int i in { pos = i; spin = Coupling.Scalar; name = denominator_name; local_array = None; momentum = "k" ^ ii; momentum_array = "p" ^ ii; fortran_type = fortran_type Coupling.Scalar }) let contract_indices ~decl ~eval indent wf_indices wfs (fusion, contractees) = let printf fmt = fprintf eval fmt and nl = pp_newline eval in let sum_var = begin match wf_indices with | [] -> wfs.(0).name | ilist -> let indices = String.concat "," ilist in begin match wfs.(0).local_array with | None -> let component = begin match wfs.(0).spin with | Coupling.Spinor | Coupling.ConjSpinor | Coupling.Majorana -> "a" | Coupling.Tensor_2 -> "t" | Coupling.Vector | Coupling.Massive_Vector -> failwith "contract_indices: expected local_array for vectors" | _ -> failwith "contract_indices: unexpected spin" end in Printf.sprintf "%s%%%s(%s)" wfs.(0).name component indices | Some a -> Printf.sprintf "%s(%s)" a indices end end in let indices = List.filter (fun i -> UFOx.Index.position i <> 1) (indices_of_contractions contractees) in nested_sums ~decl ~eval indent indices (fun indent -> printf "%*s@[<2>%s = %s" indent "" sum_var sum_var; printf "@ %s" (format_complex_rational_factor fusion.L.coeff); List.iter (fun i -> printf "@,g4_(%s)*" (index_variable i)) indices; printf "@,("; multiply_tensors ~decl ~eval contractees; printf ")"; begin match fusion.L.denominator with | [] -> () | d -> printf " / %s" denominator_name end; printf "@]"); printf "@]"; nl () let scalar_expression1 ~decl ~eval fusion = let printf fmt = fprintf eval fmt in match fusion.L.dirac, fusion.L.vector with | [], [] -> let scalars = List.map (fun t -> { L.atom = S t; L.indices = [] }) fusion.L.scalar and inverses = List.map (fun t -> { L.atom = Inv t; L.indices = [] }) fusion.L.inverse in let contractees = scalars @ inverses in printf "@ %s" (format_complex_rational_factor fusion.L.coeff); multiply_tensors ~decl ~eval contractees | _, [] -> invalid_arg "UFO_targets.Fortran.scalar_expression1: unexpected spinor indices" | [], _ -> invalid_arg "UFO_targets.Fortran.scalar_expression1: unexpected vector indices" | _, _ -> invalid_arg "UFO_targets.Fortran.scalar_expression1: unexpected indices" let scalar_expression ~decl ~eval indent name fusions = let printf fmt = fprintf eval fmt and nl = pp_newline eval in let sum_var = name in printf "%*s@[<2>%s =" indent "" sum_var; List.iter (scalar_expression1 ~decl ~eval) fusions; printf "@]"; nl () let local_vector_copies ~decl ~eval wfs = begin match wfs.(0).local_array with | None -> () | Some a -> fprintf decl " @[<2>complex(kind=default),@ dimension(0:3) ::@ %s@]" a; pp_newline decl () end; let n = Array.length wfs in for i = 1 to n - 1 do match wfs.(i).local_array with | None -> () | Some a -> fprintf decl " @[<2>complex(kind=default),@ dimension(0:3) ::@ %s@]" a; pp_newline decl (); fprintf eval " @[<2>%s(0) = %s%%t@]" a wfs.(i).name; pp_newline eval (); fprintf eval " @[<2>%s(1:3) = %s%%x@]" a wfs.(i).name; pp_newline eval () done let return_vector ff wfs = let printf fmt = fprintf ff fmt and nl = pp_newline ff in match wfs.(0).local_array with | None -> () | Some a -> pp_divide ~indent:4 ff (); printf " @[<2>%s%%t = %s(0)@]" wfs.(0).name a; nl (); printf " @[<2>%s%%x = %s(1:3)@]" wfs.(0).name a; nl () let multiply_coupling_and_scalars ff g_opt wfs = let printf fmt = fprintf ff fmt and nl = pp_newline ff in pp_divide ~indent:4 ff (); let g = match g_opt with | None -> "" | Some g -> g ^ "*" in let wfs0name = match wfs.(0).local_array with | None -> wfs.(0).name | Some a -> a in printf " @[<2>%s = %s%s" wfs0name g wfs0name; for i = 1 to Array.length wfs - 1 do match wfs.(i).spin with | Coupling.Scalar -> printf "@,*%s" wfs.(i).name | _ -> () done; printf "@]"; nl () let local_momentum_copies ~decl ~eval wfs = let n = Array.length wfs in fprintf decl " @[<2>real(kind=default),@ dimension(0:3) ::@ %s" wfs.(0).momentum_array; for i = 1 to n - 1 do fprintf decl ",@ %s" wfs.(i).momentum_array; fprintf eval " @[<2>%s(0) = %s%%t@]" wfs.(i).momentum_array wfs.(i).momentum; pp_newline eval (); fprintf eval " @[<2>%s(1:3) = %s%%x@]" wfs.(i).momentum_array wfs.(i).momentum; pp_newline eval () done; fprintf eval " @[<2>%s =" wfs.(0).momentum_array; for i = 1 to n - 1 do fprintf eval "@ - %s" wfs.(i).momentum_array done; fprintf decl "@]"; pp_newline decl (); fprintf eval "@]"; pp_newline eval () let contractees_of_fusion ~decl ~eval wfs (max_dsv, indices_seen, contractees) fusion = let max_dsv', dirac_strings = dirac_strings_to_fortran ~decl ~eval wfs max_dsv fusion.L.dirac and vectors = List.fold_left (fun acc wf -> match wf.spin, wf.local_array with | Coupling.Tensor_2, None -> { L.atom = V (Printf.sprintf "%s%d%%t" (spin_mnemonic wf.spin) wf.pos); L.indices = [UFOx.Index.pack wf.pos 1; UFOx.Index.pack wf.pos 2] } :: acc | _, None -> acc | _, Some a -> { L.atom = V a; L.indices = [wf.pos] } :: acc) [] (List.tl (Array.to_list wfs)) and tensors = List.map (L.map_atom (fun t -> T t)) fusion.L.vector and scalars = List.map (fun t -> { L.atom = S t; L.indices = [] }) fusion.L.scalar and inverses = List.map (fun t -> { L.atom = Inv t; L.indices = [] }) fusion.L.inverse in let contractees' = dirac_strings @ vectors @ tensors @ scalars @ inverses in let indices_seen' = Sets.Int.of_list (indices_of_contractions contractees') in (max_dsv', Sets.Int.union indices_seen indices_seen', (fusion, contractees') :: contractees) let local_name wf = match wf.local_array with | Some a -> a | None -> match wf.spin with | Coupling.Spinor | Coupling.ConjSpinor | Coupling.Majorana -> wf.name ^ "%a" | Coupling.Scalar -> wf.name | Coupling.Tensor_2 -> wf.name ^ "%t" | Coupling.Vector | Coupling.Massive_Vector -> failwith "UFO_targets.Fortran.local_name: unexpected spin 1" | _ -> failwith "UFO_targets.Fortran.local_name: unhandled spin" let external_wf_loop ~decl ~eval ~indent wfs (fusion, _ as contractees) = pp_divide ~indent eval (); fprintf eval "%*s! %s" indent "" (L.to_string [fusion]); pp_newline eval (); pp_divide ~indent eval (); begin match fusion.L.denominator with | [] -> () | denominator -> scalar_expression ~decl ~eval 4 denominator_name denominator end; match wfs.(0).spin with | Coupling.Scalar -> contract_indices ~decl ~eval 2 [] wfs contractees | Coupling.Spinor | Coupling.ConjSpinor | Coupling.Majorana -> let idx = index_spinor in fprintf eval "%*s@[<2>do %s = 1, 4@]" indent "" idx; pp_newline eval (); contract_indices ~decl ~eval 4 [idx] wfs contractees; fprintf eval "%*send do@]" indent ""; pp_newline eval () | Coupling.Vector | Coupling.Massive_Vector -> let idx = index_variable 1 in fprintf eval "%*s@[<2>do %s = 0, 3@]" indent "" idx; pp_newline eval (); contract_indices ~decl ~eval 4 [idx] wfs contractees; fprintf eval "%*send do@]" indent ""; pp_newline eval () | Coupling.Tensor_2 -> let idx1 = index_variable (UFOx.Index.pack 1 1) and idx2 = index_variable (UFOx.Index.pack 1 2) in fprintf eval "%*s@[<2>do %s = 0, 3@]" indent "" idx1; pp_newline eval (); fprintf eval "%*s@[<2>do %s = 0, 3@]" (indent + 2) "" idx2; pp_newline eval (); contract_indices ~decl ~eval 6 [idx1; idx2] wfs contractees; fprintf eval "%*send do@]" (indent + 2) ""; pp_newline eval (); fprintf eval "%*send do@]" indent ""; pp_newline eval () | Coupling.Vectorspinor -> failwith "external_wf_loop: Vectorspinor not supported yet!" | Coupling.Maj_Ghost -> failwith "external_wf_loop: unexpected Maj_Ghost" | Coupling.Tensor_1 -> failwith "external_wf_loop: unexpected Tensor_1" | Coupling.BRS _ -> failwith "external_wf_loop: unexpected BRS" let fusions_to_fortran ~decl ~eval wfs ?(denominator=[]) ?coupling fusions = local_vector_copies ~decl ~eval wfs; local_momentum_copies ~decl ~eval wfs; begin match denominator with | [] -> () | _ -> fprintf decl " @[<2>complex(kind=default) :: %s@]" denominator_name; pp_newline decl () end; let max_dsv, indices_used, contractions = List.fold_left (contractees_of_fusion ~decl ~eval wfs) (0, Sets.Int.empty, []) fusions in Sets.Int.iter (fun index -> fprintf decl " @[<2>integer ::@ %s@]" (index_variable index); pp_newline decl ()) indices_used; begin match wfs.(0).spin with | Coupling.Spinor | Coupling.ConjSpinor | Coupling.Majorana -> fprintf decl " @[<2>integer ::@ %s@]" index_spinor; pp_newline decl () | _ -> () end; pp_divide ~indent:4 eval (); let wfs0name = local_name wfs.(0) in fprintf eval " %s = 0" wfs0name; pp_newline eval (); List.iter (external_wf_loop ~decl ~eval ~indent:4 wfs) contractions; multiply_coupling_and_scalars eval coupling wfs; begin match denominator with | [] -> () | denominator -> pp_divide ~indent:4 eval (); fprintf eval "%*s! %s" 4 "" (L.to_string denominator); pp_newline eval (); scalar_expression ~decl ~eval 4 denominator_name denominator; fprintf eval " @[<2>%s =@ %s / %s@]" wfs0name wfs0name denominator_name; pp_newline eval () end; return_vector eval wfs (* TODO: eventually, we should include the momentum among the arguments only if required. But this can wait for another day. *) let lorentz ff name spins lorentz = let printf fmt = fprintf ff fmt and nl = pp_newline ff in let wfs = wf_table spins in let n = Array.length wfs in printf " @[<4>pure function %s@ (g,@ " name; for i = 1 to n - 2 do printf "%s,@ %s,@ " wfs.(i).name wfs.(i).momentum done; printf "%s,@ %s" wfs.(n - 1).name wfs.(n - 1).momentum; printf ")@ result (%s)@]" wfs.(0).name; nl (); printf " @[<2>%s ::@ %s@]" wfs.(0).fortran_type wfs.(0).name; nl(); printf " @[<2>complex(kind=default),@ intent(in) ::@ g@]"; nl(); for i = 1 to n - 1 do printf " @[<2>%s, intent(in) :: %s@]" wfs.(i).fortran_type wfs.(i).name; nl(); done; printf " @[<2>type(momentum), intent(in) ::@ %s" wfs.(1).momentum; for i = 2 to n - 1 do printf ",@ %s" wfs.(i).momentum done; printf "@]"; nl (); let width = 80 in (* get this from the default formatter instead! *) let decl_buf = Buffer.create 1024 and eval_buf = Buffer.create 1024 in let decl = formatter_of_buffer ~width decl_buf and eval = formatter_of_buffer ~width eval_buf in fusions_to_fortran ~decl ~eval ~coupling:"g" wfs lorentz; pp_flush decl (); pp_flush eval (); pp_divide ~indent:4 ff (); (*i printf " ! %s" (L.to_string lorentz); nl (); pp_divide ~indent:4 ff (); i*) printf "%s" (Buffer.contents decl_buf); pp_divide ~indent:4 ff (); printf " if (g == 0) then"; nl (); printf " call set_zero (%s)" wfs.(0).name; nl (); printf " return"; nl (); printf " end if"; nl (); pp_divide ~indent:4 ff (); printf "%s" (Buffer.contents eval_buf); printf " end function %s@]" name; nl (); Buffer.reset decl_buf; Buffer.reset eval_buf; () let use_variables ff parameter_module variables = let printf fmt = fprintf ff fmt and nl = pp_newline ff in match variables with | [] -> () | v :: v_list -> printf " @[<2>use %s, only: %s" parameter_module v; List.iter (fun s -> printf ", %s" s) v_list; printf "@]"; nl () let propagator ff name parameter_module variables (bra_spin, ket_spin) numerator denominator = let printf fmt = fprintf ff fmt and nl = pp_newline ff in let width = 80 in (* get this from the default formatter instead! *) let wf_name = spin_mnemonic ket_spin and wf_type = fortran_type ket_spin in let wfs = wf_table [| ket_spin; ket_spin |] in printf " @[<4>pure function pr_U_%s@ (k2, %s, %s, %s2)" name mass_name width_name wf_name; printf " result (%s1)@]" wf_name; nl (); use_variables ff parameter_module variables; printf " %s :: %s1" wf_type wf_name; nl (); printf " type(momentum), intent(in) :: k2"; nl (); printf " real(kind=default), intent(in) :: %s, %s" mass_name width_name; nl (); printf " %s, intent(in) :: %s2" wf_type wf_name; nl (); let decl_buf = Buffer.create 1024 and eval_buf = Buffer.create 1024 in let decl = formatter_of_buffer ~width decl_buf and eval = formatter_of_buffer ~width eval_buf in fusions_to_fortran ~decl ~eval wfs ~denominator numerator; pp_flush decl (); pp_flush eval (); pp_divide ~indent:4 ff (); printf "%s" (Buffer.contents decl_buf); pp_divide ~indent:4 ff (); printf "%s" (Buffer.contents eval_buf); printf " end function pr_U_%s@]" name; nl (); Buffer.reset decl_buf; Buffer.reset eval_buf; () let scale_coupling c g = if c = 1 then g else if c = -1 then "-" ^ g else Printf.sprintf "%d*%s" c g let scale_coupling z g = format_complex_rational_factor z ^ g (* As a prototypical example consider the vertex \begin{subequations} \label{eq:cyclic-UFO-fusions} \begin{equation} \bar\psi\fmslash{A}\psi = \tr\left(\psi\otimes\bar\psi\fmslash{A}\right) \end{equation} encoded as \texttt{FFV} in the SM UFO file. This example is useful, because all three fields have different type and we can use the Fortran compiler to check our implementation. In this case we need to generate the following function calls with the arguments in the following order \begin{center} \begin{tabular}{lcl} \texttt{F12}:&$\psi_1\bar\psi_2\to A$& \texttt{FFV\_p201(g,psi1,p1,psibar2,p2)} \\ \texttt{F21}:&$\bar\psi_1\psi_2\to A$& \texttt{FFV\_p201(g,psi2,p2,psibar1,p1)} \\ \texttt{F23}:&$\bar\psi_1 A_2 \to \bar\psi$& \texttt{FFV\_p012(g,psibar1,p1,A2,p2)} \\ \texttt{F32}:&$A_1\bar\psi_2 \to \bar\psi$& \texttt{FFV\_p012(g,psibar2,p2,A1,p1)} \\ \texttt{F31}:&$A_1\psi_2\to \psi$& \texttt{FFV\_p120(g,A1,p1,psi2,p2)} \\ \texttt{F13}:&$\psi_1A_2\to \psi$& \texttt{FFV\_p120(g,A2,p2,psi1,p1)} \end{tabular} \end{center} *) (* Fortunately, all Fermi signs have been taken care of by [Fusions] and we can concentrate on injecting the wave functions into the correct slots. *) (* The other possible cases are \begin{equation} \bar\psi\fmslash{A}\psi \end{equation} which would be encoded as \texttt{FVF} in a UFO file \begin{center} \begin{tabular}{lcl} \texttt{F12}:&$\bar\psi_1 A_2 \to \bar\psi$& \texttt{FVF\_p201(g,psibar1,p1,A2,p2)} \\ \texttt{F21}:&$A_1\bar\psi_2 \to \bar\psi$& \texttt{FVF\_p201(g,psibar2,p2,A1,p1)} \\ \texttt{F23}:&$A_1\psi_2\to \psi$& \texttt{FVF\_p012(g,A1,p1,psi2,p2)} \\ \texttt{F32}:&$\psi_1A_2\to \psi$& \texttt{FVF\_p012(g,A2,p2,psi1,p1)} \\ \texttt{F31}:&$\psi_1\bar\psi_2\to A$& \texttt{FVF\_p120(g,psi1,p1,psibar2,p2)} \\ \texttt{F13}:&$\bar\psi_1\psi_2\to A$& \texttt{FVF\_p120(g,psi2,p2,psibar1,p1)} \end{tabular} \end{center} and \begin{equation} \bar\psi\fmslash{A}\psi = \tr\left(\fmslash{A}\psi\otimes\bar\psi\right)\,, \end{equation} corresponding to \texttt{VFF} \begin{center} \begin{tabular}{lcl} \texttt{F12}:&$A_1\psi_2\to \psi$& \texttt{VFF\_p201(g,A1,p1,psi2,p2)} \\ \texttt{F21}:&$\psi_1A_2\to \psi$& \texttt{VFF\_p201(g,A2,p2,psi1,p1)} \\ \texttt{F23}:&$\psi_1\bar\psi_2\to A$& \texttt{VFF\_p012(g,psi1,p1,psibar2,p2)} \\ \texttt{F32}:&$\bar\psi_1\psi_2\to A$& \texttt{VFF\_p012(g,psi2,p2,psibar1,p1)} \\ \texttt{F31}:&$\bar\psi_1 A_2 \to \bar\psi$& \texttt{VFF\_p120(g,psibar1,p1,A2,p2)} \\ \texttt{F13}:&$A_1\bar\psi_2 \to \bar\psi$& \texttt{VFF\_p120(g,psibar2,p2,A1,p1)} \end{tabular} \end{center} \end{subequations} *) (* \begin{dubious} Once the Majorana code generation is fully debugged, we should replace the lists by reverted lists everywhere in order to become a bit more efficient. \end{dubious} *) module P = Permutation.Default let factor_cyclic f12__n = let f12__, fn = ThoList.split_last f12__n in let cyclic = ThoList.cycle_until fn (List.sort compare f12__n) in (P.of_list (List.map pred cyclic), P.of_lists (List.tl cyclic) f12__) let ccs_to_string ccs = String.concat "" (List.map (fun (f, i) -> Printf.sprintf "_c%x%x" i f) ccs) let fusion_name v perm ccs = Printf.sprintf "%s_p%s%s" v (P.to_string perm) (ccs_to_string ccs) let fuse_dirac c v s fl g wfs ps fusion = let g = scale_coupling c g and cyclic, factor = factor_cyclic fusion in let wfs_ps = List.map2 (fun wf p -> (wf, p)) wfs ps in let args = P.list (P.inverse factor) wfs_ps in let args_string = String.concat "," (List.map (fun (wf, p) -> wf ^ "," ^ p) args) in printf "%s(%s,%s)" (fusion_name v cyclic []) g args_string (* We need to look at the permuted fermion lines in order to decide wether to apply charge conjugations. *) (* It is not enough to look at the cyclic permutation used to move the fields into the correct arguments of the fusions \ldots *) let map_indices perm unit = let pmap = IntPM.of_lists unit (P.list perm unit) in IntPM.apply pmap (* \ldots{} we also need to inspect the full permutation of the fields. *) let map_indices2 perm unit = let pmap = IntPM.of_lists unit (1 :: P.list (P.inverse perm) (List.tl unit)) in IntPM.apply pmap (* This is a more direct implementation of the composition of [map_indices2] and [map_indices], that is used in the unit tests. *) let map_indices_raw fusion = let unit = ThoList.range 1 (List.length fusion) in let f12__, fn = ThoList.split_last fusion in let fusion = fn :: f12__ in let map_index = IntPM.of_lists fusion unit in IntPM.apply map_index (* Map the fermion line indices in [fl] according to [map_index]. *) let map_fermion_lines map_index fl = List.map (fun (i, f) -> (map_index i, map_index f)) fl (* Map the fermion line indices in [fl] according to [map_index], but keep a copy of the original. *) let map_fermion_lines2 map_index fl = List.map (fun (i, f) -> ((i, f), (map_index i, map_index f))) fl let permute_fermion_lines cyclic unit fl = map_fermion_lines (map_indices cyclic unit) fl let permute_fermion_lines2 cyclic factor unit fl = map_fermion_lines2 (map_indices2 factor unit) (map_fermion_lines (map_indices cyclic unit) fl) (* \begin{dubious} TODO: this needs more more work for the fully - general case. + general case with 4-fermion operators involving Majoranas. \end{dubious} *) let charge_conjugations fl2 = ThoList.filtermap (fun ((i, f), (i', f')) -> match (i, f), (i', f') with | (1, 2), _ | (2, 1), _ -> Some (f, i) (* $\chi^T\Gamma'$ *) | _, (2, 3) -> Some (f, i) (* $\chi^T(C\Gamma')\chi$ *) | _ -> None) fl2 (*i let fuse_majorana c v s fl g wfs ps fusion = let g = scale_coupling c g and cyclic, factor = factor_cyclic fusion in let wfs_ps = List.map2 (fun wf p -> (wf, p)) wfs ps in let wfs_ps_string = String.concat "," (List.map (fun (wf, p) -> wf ^ "," ^ p) wfs_ps) in let args = P.list (P.inverse factor) wfs_ps in let args_string = String.concat "," (List.map (fun (wf, p) -> wf ^ "," ^ p) args) in let f12__, fn = ThoList.split_last fusion in Printf.eprintf "fusion : %d < %s\n" fn (ThoList.to_string string_of_int f12__); Printf.eprintf "cyclic : %s\n" (P.to_string cyclic); Printf.eprintf "factor : %s\n" (P.to_string factor); let unit = ThoList.range 1 (List.length fusion) in Printf.eprintf "permutation : %s -> %s\n" (ThoList.to_string string_of_int unit) (ThoList.to_string string_of_int (List.map (map_indices cyclic unit) unit)); Printf.eprintf "permutation raw : %s -> %s\n" (ThoList.to_string string_of_int unit) (ThoList.to_string string_of_int (List.map (map_indices_raw fusion) unit)); Printf.eprintf "fermion lines : %s\n" (ThoList.to_string (fun (i, f) -> Printf.sprintf "%d>%d" i f) fl); let fl2 = permute_fermion_lines2 cyclic factor unit fl in let fl = permute_fermion_lines cyclic unit fl in Printf.eprintf "permuted : %s\n" (ThoList.to_string (fun (i, f) -> Printf.sprintf "%d>%d" i f) fl); Printf.eprintf "arguments : %s\n" wfs_ps_string; Printf.eprintf "permuted : %s\n" args_string; Printf.eprintf ">> %s(%s,%s)\n" (fusion_name v cyclic (charge_conjugations fl2)) g args_string; printf "%s(%s,%s)" (fusion_name v cyclic (charge_conjugations fl2)) g args_string i*) let charge_conjugations fl2 = ThoList.filtermap (fun ((i, f), (i', f')) -> match (i, f), (i', f') with | _, (2, 3) -> Some (f, i) | _ -> None) fl2 let fuse_majorana c v s fl g wfs ps fusion = let g = scale_coupling c g and cyclic, factor = factor_cyclic fusion in let wfs_ps = List.map2 (fun wf p -> (wf, p)) wfs ps in let args = P.list (P.inverse factor) wfs_ps in let args_string = String.concat "," (List.map (fun (wf, p) -> wf ^ "," ^ p) args) in let unit = ThoList.range 1 (List.length fusion) in let ccs = charge_conjugations (permute_fermion_lines2 cyclic factor unit fl) in printf "%s(%s,%s)" (fusion_name v cyclic ccs) g args_string let fuse c v s fl g wfs ps fusion = if List.exists is_majorana s then fuse_majorana c v s fl g wfs ps fusion else fuse_dirac c v s fl g wfs ps fusion let eps4_g4_g44_decl ff () = let printf fmt = fprintf ff fmt and nl = pp_newline ff in printf " @[<2>integer,@ dimension(0:3)"; printf ",@ save,@ private ::@ g4_@]"; nl (); printf " @[<2>integer,@ dimension(0:3,0:3)"; printf ",@ save,@ private ::@ g44_@]"; nl (); printf " @[<2>integer,@ dimension(0:3,0:3,0:3,0:3)"; printf ",@ save,@ private ::@ eps4_@]"; nl () let eps4_g4_g44_init ff () = let printf fmt = fprintf ff fmt and nl = pp_newline ff in printf " @[<2>data g4_@ /@ 1, -1, -1, -1 /@]"; nl (); printf " @[<2>data g44_(0,:)@ /@ 1, 0, 0, 0 /@]"; nl (); printf " @[<2>data g44_(1,:)@ /@ 0, -1, 0, 0 /@]"; nl (); printf " @[<2>data g44_(2,:)@ /@ 0, 0, -1, 0 /@]"; nl (); printf " @[<2>data g44_(3,:)@ /@ 0, 0, 0, -1 /@]"; nl (); for mu1 = 0 to 3 do for mu2 = 0 to 3 do for mu3 = 0 to 3 do printf " @[<2>data eps4_(%d,%d,%d,:)@ /@ " mu1 mu2 mu3; for mu4 = 0 to 3 do if mu4 <> 0 then printf ",@ "; let mus = [mu1; mu2; mu3; mu4] in if List.sort compare mus = [0; 1; 2; 3] then printf "%2d" (Combinatorics.sign mus) else printf "%2d" 0; done; printf " /@]"; nl () done done done let inner_product_functions ff () = let printf fmt = fprintf ff fmt and nl = pp_newline ff in printf " pure function g2_ (p) result (p2)"; nl(); printf " real(kind=default), dimension(0:3), intent(in) :: p"; nl(); printf " real(kind=default) :: p2"; nl(); printf " p2 = p(0)*p(0) - p(1)*p(1) - p(2)*p(2) - p(3)*p(3)"; nl(); printf " end function g2_"; nl(); printf " pure function g12_ (p1, p2) result (p12)"; nl(); printf " real(kind=default), dimension(0:3), intent(in) :: p1, p2"; nl(); printf " real(kind=default) :: p12"; nl(); printf " p12 = p1(0)*p2(0) - p1(1)*p2(1) - p1(2)*p2(2) - p1(3)*p2(3)"; nl(); printf " end function g12_"; nl() module type Test = sig val suite : OUnit.test end module Test : Test = struct open OUnit let assert_mappings fusion = let unit = ThoList.range 1 (List.length fusion) in let cyclic, factor = factor_cyclic fusion in let raw = map_indices_raw fusion and map1 = map_indices cyclic unit and map2 = map_indices2 factor unit in let map i = map2 (map1 i) in assert_equal ~printer:(ThoList.to_string string_of_int) (List.map raw unit) (List.map map unit) let suite_mappings = "mappings" >::: [ "1<-2" >:: (fun () -> List.iter assert_mappings (Combinatorics.permute [1;2;3])); "1<-3" >:: (fun () -> List.iter assert_mappings (Combinatorics.permute [1;2;3;4])) ] let suite = "UFO_targets" >::: [suite_mappings] end end