Index: trunk/synchronize.sh =================================================================== --- trunk/synchronize.sh (revision 8773) +++ trunk/synchronize.sh (revision 8774) @@ -1,60 +1,60 @@ #!/bin/sh ### Consider it safer to explicitly mention all files that contain ### email addresses or copyright tags. OLD_YEAR="Copyright (C) 1999-2020"; NEW_YEAR="Copyright (C) 1999-2021"; OLD_YEAR2="Copyright (C) 2001-2020"; NEW_YEAR2="Copyright (C) 2001-2021"; OLD_YEAR3="Copyright (C) 2019-2020"; NEW_YEAR3="Copyright (C) 2019-2021"; # OLD_ADDRESS="Soyoung Shim " # NEW_ADDRESS="So Young Shim " OLD_ADDRESS="Soyoung Shim" NEW_ADDRESS="So Young Shim" OLD_DATE="Jul 08 2021" NEW_DATE="Nov 23 2021" OLD_VERSION="3.0.2" NEW_VERSION="3.0.2+" #OLD_STATUS="PACKAGE_STATUS=\"alpha\"" #NEW_STATUS="PACKAGE_STATUS=\"beta\"" OLD_STATUS="PACKAGE_STATUS=\"release\"" #NEW_STATUS="PACKAGE_STATUS=\"rc1\"" NEW_STATUS="PACKAGE_STATUS=\"alpha\"" ## We should add an option to add an author here. ## share/doc/manual.tex should be changed manually ## We have to discuss the entries in gamelan/manual ## We have to discuss the entries in src/shower MAIN_FILES="AUTHORS BUGS Makefile.am.in README build_master.sh tests/Makefile.am tests/models/Makefile.am tests/models/UFO/Makefile.am tests/models/UFO/SM/Makefile.am tests/models/UFO/MSSM/Makefile.am tests/functional_tests/Makefile.am tests/ext_tests_mssm/Makefile.am tests/ext_tests_nmssm/Makefile.am tests/ext_tests_ilc/Makefile.am tests/ext_tests_nlo/Makefile.am tests/ext_tests_nlo_add/Makefile.am tests/ext_tests_shower/Makefile.am tests/unit_tests/Makefile.am" CONFIGURE_FILES="configure.ac.in src/noweb-frame/whizard-prelude.nw" VERSION_FILES="NEWS circe2/src/circe2.nw" SCRIPTS_FILES="scripts/Makefile.am scripts/whizard-config.in scripts/whizard-setup.csh.in scripts/whizard-setup.sh.in" SHARE_FILES="share/Makefile.am share/doc/Makefile.am share/doc/custom.hva share/examples/NLO_eettbar_GoSam.sin share/examples/NLO_eettbar_OpenLoops.sin share/examples/HERA_DIS.sin share/examples/LEP_cc10.sin share/examples/LEP_higgs.sin share/examples/W-endpoint.sin share/examples/Z-lineshape.sin share/examples/Zprime.sin share/examples/casc_dec.sin share/examples/circe1.sin share/examples/eeww_polarized.sin share/examples/DrellYanMatchingP.sin share/examples/DrellYanMatchingW.sin share/examples/DrellYanNoMatchingP.sin share/examples/DrellYanNoMatchingW.sin share/examples/EEMatching2P.sin share/examples/EEMatching2W.sin share/examples/EEMatching3P.sin share/examples/EEMatching3W.sin share/examples/EEMatching4P.sin share/examples/EEMatching4W.sin share/examples/EEMatching5P.sin share/examples/EEMatching5W.sin share/examples/EENoMatchingP.sin share/examples/EENoMatchingW.sin share/examples/LHC_VBS_likesign.sin share/tests/Makefile.am" -SRC_FILES="src/Makefile.am src/feynmf/Makefile.am src/hepmc/Makefile.am src/hepmc/HepMCWrap_dummy.f90 src/lcio/Makefile.am src/lcio/LCIOWrap_dummy.f90 src/tauola/Makefile.am src/lhapdf/Makefile.am src/lhapdf/lhapdf.f90 src/lhapdf5/Makefile.am src/pdf_builtin/Makefile.am src/pdf_builtin/pdf_builtin.f90 src/qed_pdf/Makefile.am src/qed_pdf/qed_pdf.nw src/fastjet/Makefile.am src/fastjet/cpp_strings.f90 src/fastjet/fastjet.f90 src/fastjet/Makefile.am src/hoppet/Makefile.am src/hoppet/hoppet.f90 pythia6/Makefile.am tauola/Makefile.am mcfio/Makefile.am stdhep/Makefile.am src/noweb-frame/Makefile.am src/noweb-frame/whizard-prelude.nw src/noweb-frame/whizard-postlude.nw src/utilities/Makefile.am src/matrix_elements/Makefile.am src/matrix_elements/matrix_elements.nw src/mci/Makefile.am src/vegas/Makefile.am src/vegas/vegas.nw src/mci/mci.nw src/utilities/utilities.nw src/testing/Makefile.am src/testing/testing.nw src/system/Makefile.am src/system/system.nw src/system/system_dependencies.f90.in src/system/debug_master.f90.in src/combinatorics/Makefile.am src/combinatorics/combinatorics.nw src/parsing/Makefile.am src/parsing/parsing.nw src/particles/Makefile.am src/particles/particles.nw src/phase_space/Makefile.am src/phase_space/phase_space.nw src/physics/Makefile.am src/physics/physics.nw src/beams/Makefile.am src/beams/beams.nw src/qft/Makefile.am src/qft/qft.nw src/rng/Makefile.am src/rng/rng.nw src/types/Makefile.am src/types/types.nw src/whizard-core/Makefile.am src/whizard-core/whizard.nw src/pythia8/Makefile.am src/shower/Makefile.am src/shower/shower.nw src/muli/Makefile.am src/muli/muli.nw src/model_features/model_features.nw src/model_features/Makefile.am src/me_methods/Makefile.am src/me_methods/me_methods.nw src/gosam/Makefile.am src/gosam/gosam.nw src/fks/Makefile.am src/fks/fks.nw src/expr_base/Makefile.am src/expr_base/expr_base.nw src/events/Makefile.am src/events/events.nw src/blha/Makefile.am src/blha/blha.nw src/variables/Makefile.am src/variables/variables.nw src/xdr/Makefile.am src/xdr/xdr_wo_stdhep.f90 src/looptools/Makefile.am src/process_integration/Makefile.am src/process_integration/process_integration.nw src/matching/Makefile.am src/matching/matching.nw src/openloops/Makefile.am src/openloops/openloops.nw src/recola/Makefile.am src/recola/recola.nw src/transforms/Makefile.am src/transforms/transforms.nw src/threshold/Makefile.am src/threshold/threshold.nw src/api/Makefile.am src/api/api.nw src/main/Makefile.am src/main/main.nw" +SRC_FILES="src/Makefile.am src/feynmf/Makefile.am src/hepmc/Makefile.am src/hepmc/HepMCWrap_dummy.f90 src/lcio/Makefile.am src/lcio/LCIOWrap_dummy.f90 src/tauola/Makefile.am src/lhapdf/Makefile.am src/lhapdf/lhapdf.f90 src/lhapdf5/Makefile.am src/pdf_builtin/Makefile.am src/pdf_builtin/pdf_builtin.f90 src/pdf/pdf_builtin_sub.f90 src/qed_pdf/Makefile.am src/qed_pdf/qed_pdf.nw src/fastjet/Makefile.am src/fastjet/cpp_strings.f90 src/fastjet/fastjet.f90 src/fastjet/Makefile.am src/hoppet/Makefile.am src/hoppet/hoppet.f90 pythia6/Makefile.am tauola/Makefile.am mcfio/Makefile.am stdhep/Makefile.am src/noweb-frame/Makefile.am src/noweb-frame/whizard-prelude.nw src/noweb-frame/whizard-postlude.nw src/utilities/Makefile.am src/matrix_elements/Makefile.am src/matrix_elements/matrix_elements.nw src/mci/Makefile.am src/vegas/Makefile.am src/vegas/vegas.nw src/mci/mci.nw src/utilities/utilities.nw src/testing/Makefile.am src/testing/testing.nw src/system/Makefile.am src/system/system.nw src/system/system_dependencies.f90.in src/system/debug_master.f90.in src/combinatorics/Makefile.am src/combinatorics/combinatorics.nw src/parsing/Makefile.am src/parsing/parsing.nw src/particles/Makefile.am src/particles/particles.nw src/phase_space/Makefile.am src/phase_space/phase_space.nw src/physics/Makefile.am src/physics/physics.nw src/beams/Makefile.am src/beams/beams.nw src/qft/Makefile.am src/qft/qft.nw src/rng/Makefile.am src/rng/rng.nw src/types/Makefile.am src/types/types.nw src/whizard-core/Makefile.am src/whizard-core/whizard.nw src/pythia8/Makefile.am src/shower/Makefile.am src/shower/shower.nw src/muli/Makefile.am src/muli/muli.nw src/model_features/model_features.nw src/model_features/Makefile.am src/me_methods/Makefile.am src/me_methods/me_methods.nw src/gosam/Makefile.am src/gosam/gosam.nw src/fks/Makefile.am src/fks/fks.nw src/expr_base/Makefile.am src/expr_base/expr_base.nw src/events/Makefile.am src/events/events.nw src/blha/Makefile.am src/blha/blha.nw src/variables/Makefile.am src/variables/variables.nw src/xdr/Makefile.am src/xdr/xdr_wo_stdhep.f90 src/looptools/Makefile.am src/process_integration/Makefile.am src/process_integration/process_integration.nw src/matching/Makefile.am src/matching/matching.nw src/openloops/Makefile.am src/openloops/openloops.nw src/recola/Makefile.am src/recola/recola.nw src/transforms/Makefile.am src/transforms/transforms.nw src/threshold/Makefile.am src/threshold/threshold.nw src/api/Makefile.am src/api/api.nw src/main/Makefile.am src/main/main.nw" CIRCE1_FILES="circe1/Makefile.am circe1/share/Makefile.am circe1/share/doc/Makefile.am circe1/src/Makefile.am circe1/src/circe1.nw circe1/minuit/Makefile.am circe1/src/minuit.nw circe1/tools/Makefile.am" CIRCE2_FILES="circe2/Makefile.am circe2/share/Makefile.am circe2/share/doc/Makefile.am circe2/src/Makefile.am circe2/src/Makefile.ocaml circe2/src/circe2.nw circe2/src/Makefile.sources circe2/src/postlude.nw circe2/tests/Makefile.am circe2/src/circe2_tool.ml circe2/src/commands.ml circe2/src/commands.mli circe2/src/diffmap.ml circe2/src/diffmap.mli circe2/src/diffmaps.ml circe2/src/diffmaps.mli circe2/src/division.ml circe2/src/division.mli circe2/src/events.ml circe2/src/events.mli circe2/src/filter.ml circe2/src/filter.mli circe2/src/float.ml circe2/src/float.mli circe2/src/grid.ml circe2/src/grid.mli circe2/src/histogram.mli circe2/src/histogram.ml circe2/src/syntax.ml circe2/src/syntax.mli circe2/src/thoArray.ml circe2/src/thoArray.mli circe2/src/thoMatrix.ml circe2/src/thoMatrix.mli circe2/src/bigarray_module.ml circe2/src/bigarray_library.ml circe2/src/bigarray_compat.mli" SRC_GAMELAN_FILES="src/gamelan/Makefile.am src/gamelan/whizard-gml.in" SRC_BASICS_FILES="src/basics/constants.f90 src/basics/io_units.f90 src/basics/Makefile.am" SRC_MODELS_FILES="src/models/threeshl_bundle/Makefile.am src/models/threeshl_bundle/threeshl_bundle.f90 src/models/threeshl_bundle/threeshl_bundle_lt.f90 src/models/external.Test.f90 src/models/external.Threeshl.f90 src/models/external.SM_tt_threshold.f90 src/models/Makefile.am src/models/parameters.THDM.f90 src/models/parameters.GravTest.f90 src/models/parameters.Littlest.f90 src/models/parameters.Littlest_Eta.f90 src/models/parameters.Littlest_Tpar.f90 src/models/parameters.MSSM.f90 src/models/parameters.MSSM_4.f90 src/models/parameters.MSSM_CKM.f90 src/models/parameters.MSSM_Grav.f90 src/models/parameters.MSSM_Hgg.f90 src/models/parameters.NMSSM.f90 src/models/parameters.NMSSM_CKM.f90 src/models/parameters.NMSSM_Hgg.f90 src/models/parameters.PSSSM.f90 src/models/parameters.QCD.f90 src/models/parameters.QED.f90 src/models/parameters.SM.f90 src/models/parameters.SM_CKM.f90 src/models/parameters.SM_ac.f90 src/models/parameters.SM_ac_CKM.f90 src/models/parameters.SM_dim6.f90 src/models/parameters.SM_rx.f90 src/models/parameters.SM_ul.f90 src/models/parameters.NoH_rx.f90 src/models/parameters.AltH.f90 src/models/parameters.SSC.f90 src/models/parameters.SSC_2.f90 src/models/parameters.SSC_AltT.f90 src/models/parameters.SM_top.f90 src/models/parameters.SM_top_anom.f90 src/models/parameters.SM_Higgs.f90 src/models/parameters.SM_Higgs_CKM.f90 src/models/parameters.SM_tt_threshold.f90 src/models/parameters.Simplest.f90 src/models/parameters.Simplest_univ.f90 src/models/parameters.Template.f90 src/models/parameters.HSExt.f90 src/models/parameters.Test.f90 src/models/parameters.Threeshl.f90 src/models/parameters.UED.f90 src/models/parameters.Xdim.f90 src/models/parameters.Zprime.f90 src/models/parameters.WZW.f90" OMEGA_FILES="omega/Makefile.am omega/share/Makefile.am omega/share/doc/Makefile.am omega/src/Makefile.am omega/src/Makefile.ocaml omega/src/Makefile.sources omega/bin/Makefile.am omega/extensions/Makefile.am omega/extensions/people/Makefile.am omega/extensions/people/jr/Makefile.am omega/extensions/people/jr/f90_SAGT.ml omega/extensions/people/jr/f90_SQED.ml omega/extensions/people/jr/f90_WZ.ml omega/extensions/people/tho/Makefile.am omega/extensions/people/tho/f90_O2.ml omega/lib/Makefile.am omega/models/Makefile.am omega/scripts/Makefile.am omega/scripts/omega-config.in omega/tools/Makefile.am omega/tests/parameters_QED.f90 omega/tests/parameters_QCD.f90 omega/tests/parameters_SM.f90 omega/tests/parameters_SM_CKM.f90 omega/tests/parameters_SM_Higgs.f90 omega/tests/parameters_SM_from_UFO.f90 omega/tests/parameters_SYM.f90 omega/tests/parameters_SM_top_anom.f90 omega/tests/parameters_HSExt.f90 omega/tests/parameters_THDM.f90 omega/tests/parameters_THDM_CKM.f90 omega/tests/parameters_Zprime.f90 omega/tests/test_openmp.f90 omega/tests/tao_random_numbers.f90 omega/tests/test_qed_eemm.f90 omega/tests/Makefile.am omega/tests/benchmark.f90 omega/tests/color_test_lib.f90 omega/tests/omega_interface.f90 omega/tests/ward_lib.f90 omega/tests/omega_unit.ml omega/tests/compare_lib.f90 omega/tests/compare_lib_recola.f90 omega/tests/benchmark_UFO_SM.f90 omega/tests/benchmark_UFO_SMEFT.f90 omega/tests/keystones_omegalib_generate.ml omega/tests/keystones_UFO_generate.ml omega/tests/keystones_omegalib_bispinors_generate.ml omega/tests/keystones_UFO_bispinors_generate.ml omega/tests/keystones.ml omega/tests/keystones.mli omega/tests/keystones_tools.f90 omega/tests/fermi_lib.f90 omega/tests/parameters_SM_Higgs_recola.f90 omega/tests/parameters_MSSM.f90 omega/tests/keystones.mli" OMEGA_SRC_FILES="omega/src/algebra.ml omega/src/algebra.mli omega/src/bundle.ml omega/src/bundle.mli omega/src/cache.ml omega/src/cache.mli omega/src/cascade.ml omega/src/cascade.mli omega/src/cascade_lexer.mll omega/src/cascade_parser.mly omega/src/cascade_syntax.ml omega/src/cascade_syntax.mli omega/src/charges.ml omega/src/charges.mli omega/src/color.ml omega/src/color.mli omega/src/colorize.ml omega/src/colorize.mli omega/src/combinatorics.ml omega/src/combinatorics.mli omega/src/complex.ml omega/src/complex.mli omega/src/config.ml.in omega/src/config.mli omega/src/count.ml omega/src/coupling.mli omega/src/DAG.ml omega/src/DAG.mli omega/src/fusion.ml omega/src/fusion_vintage.ml omega/src/fusion.mli omega/src/fusion_vintage.mli omega/src/linalg.ml omega/src/linalg.mli omega/src/model.mli omega/src/modellib_BSM.ml omega/src/modellib_NoH.ml omega/src/modellib_NoH.mli omega/src/modellib_BSM.mli omega/src/modellib_MSSM.ml omega/src/modellib_MSSM.mli omega/src/modellib_NMSSM.ml omega/src/modellib_NMSSM.mli omega/src/modellib_PSSSM.ml omega/src/modellib_PSSSM.mli omega/src/modellib_SM.ml omega/src/modellib_SM.mli omega/src/modellib_Zprime.mli omega/src/modellib_Zprime.ml omega/src/modellib_WZW.mli omega/src/modellib_WZW.ml omega/src/UFO.ml omega/src/UFO.mli omega/src/UFO_targets.ml omega/src/UFO_Lorentz.ml omega/src/UFO_syntax.ml omega/src/UFO_syntax.mli omega/src/UFOx.ml omega/src/UFOx.mli omega/src/UFO_lexer.mll omega/src/UFO_parser.mly omega/src/UFOx_syntax.ml omega/src/UFOx_syntax.mli omega/src/UFOx_lexer.mll omega/src/UFOx_parser.mly omega/src/omega_UFO.ml omega/src/modeltools.ml omega/src/modeltools.mli omega/src/momentum.ml omega/src/momentum.mli omega/src/OVM.ml omega/src/OVM.mli omega/src/ogiga.ml omega/src/omega.ml omega/src/omega.mli omega/src/omega_THDM.ml omega/src/omega_THDM_VM.ml omega/src/omega_THDM_CKM.ml omega/src/omega_THDM_CKM_VM.ml omega/src/omega_CQED.ml omega/src/omega_GravTest.ml omega/src/omega_Littlest.ml omega/src/omega_Littlest_Eta.ml omega/src/omega_Littlest_Tpar.ml omega/src/omega_Littlest_Zprime.ml omega/src/omega_MSSM.ml omega/src/omega_MSSM_CKM.ml omega/src/omega_MSSM_Grav.ml omega/src/omega_MSSM_Hgg.ml omega/src/omega_NMSSM.ml omega/src/omega_NMSSM_CKM.ml omega/src/omega_NMSSM_Hgg.ml omega/src/omega_PSSSM.ml omega/src/omega_Phi3.ml omega/src/omega_Phi3h.ml omega/src/omega_Phi4.ml omega/src/omega_Phi4h.ml omega/src/omega_QCD.ml omega/src/omega_QCD_VM.ml omega/src/omega_QED.ml omega/src/omega_QED_VM.ml omega/src/omega_SM.ml omega/src/omega_SM_tt_threshold.ml omega/src/omega_SM_VM.ml omega/src/omega_SM_CKM.ml omega/src/omega_SM_CKM_VM.ml omega/src/ovm_SM.ml omega/src/process.ml omega/src/process.mli omega/src/thoFilename.ml omega/src/thoFilename.mli omega/src/omega_SM_Higgs.ml omega/src/omega_SM_Higgs_CKM.ml omega/src/omega_SM_Higgs_VM.ml omega/src/omega_SM_Higgs_CKM_VM.ml omega/src/omega_SM_Rxi.ml omega/src/omega_SM_ac.ml omega/src/omega_SM_ac_CKM.ml omega/src/omega_SM_clones.ml omega/src/omega_SM_rx.ml omega/src/omega_SM_ul.ml omega/src/omega_SM_Majorana_legacy.ml omega/src/omega_SM_Majorana.ml omega/src/omega_NoH_rx.ml omega/src/omega_AltH.ml omega/src/omega_SSC.ml omega/src/omega_SSC_2.ml omega/src/omega_SM_top.ml omega/src/omega_SM_top_anom.ml omega/src/omega_SMh.ml omega/src/omega_SYM.ml omega/src/omega_Simplest.ml omega/src/omega_Simplest_univ.ml omega/src/omega_Template.ml omega/src/omega_HSExt.ml omega/src/omega_HSExt_VM.ml omega/src/omega_Threeshl.ml omega/src/omega_Threeshl_nohf.ml omega/src/omega_UED.ml omega/src/omega_Xdim.ml omega/src/omega_Zprime.ml omega/src/omega_Zprime_VM.ml omega/src/omega_logo.mp omega/src/omega_parameters_tool.nw omega/src/omegalib.nw omega/src/options.ml omega/src/options.mli omega/src/partition.ml omega/src/partition.mli omega/src/phasespace.ml omega/src/phasespace.mli omega/src/pmap.ml omega/src/pmap.mli omega/src/powSet.ml omega/src/powSet.mli omega/src/product.ml omega/src/product.mli omega/src/progress.ml omega/src/progress.mli omega/src/permutation.ml omega/src/permutation.mli omega/src/target.mli omega/src/targets.ml omega/src/targets.mli omega/src/targets_Kmatrix.ml omega/src/targets_Kmatrix.mli omega/src/test_linalg.ml omega/src/thoArray.ml omega/src/thoFilename.ml omega/src/thoArray.mli omega/src/thoGButton.ml omega/src/thoGButton.mli omega/src/thoGDraw.ml omega/src/thoGDraw.mli omega/src/thoGMenu.ml omega/src/thoGMenu.mli omega/src/thoGWindow.ml omega/src/thoGWindow.mli omega/src/thoList.ml omega/src/thoList.mli omega/src/thoString.ml omega/src/thoString.mli omega/src/topology.ml omega/src/topology.mli omega/src/tree.ml omega/src/tree.mli omega/src/tree2.ml omega/src/tree2.mli omega/src/trie.ml omega/src/trie.mli omega/src/tuple.ml omega/src/tuple.mli omega/src/vertex.ml omega/src/vertex.mli omega/src/vertex_lexer.mll omega/src/vertex_parser.mly omega/src/vertex_syntax.ml omega/src/vertex_syntax.mli omega/src/whizard.ml omega/src/whizard.mli omega/src/whizard_tool.ml omega/src/constants.f90 omega/src/sets.mli omega/src/sets.ml omega/src/UFO_tools.ml omega/src/UFO_tools.mli omega/src/fortran_unit.ml omega/src/format_Fortran.ml omega/src/format_Fortran.mli omega/src/omega_UFO_Majorana.ml omega/src/omega_UFO_Dirac.ml" SRC_PDF_BUILTIN_FILES="src/pdf_builtin/pdf_builtin.f90" VAMP_FILES="vamp/Makefile.am vamp/share/Makefile.am vamp/share/doc/Makefile.am vamp/src/Makefile.am vamp/tests/Makefile.am" FILES="$MAIN_FILES $CONFIGURE_FILES $VERSION_FILES $SHARE_FILES $OMEGA_FILES $SCRIPTS_FILES $SRC_FILES $CIRCE1_FILES $CIRCE2_FILES $SRC_GAMELAN_FILES $SRC_PDF_BUILTIN_FILES $VAMP_FILES $SRC_BASICS_FILES $SRC_MODELS_FILES $OMEGA_SRC_FILES" for f in $FILES; do sed -e "s/$OLD_YEAR/$NEW_YEAR/g" -e "s/$OLD_YEAR2/$NEW_YEAR2/g" -e "s/$OLD_YEAR3/$NEW_YEAR3/g" $f > $f.tmp; cp -f $f.tmp $f; rm -f $f.tmp; done CHANGE_FILES="$CONFIGURE_FILES $VERSION_FILES" for f in $CHANGE_FILES; do sed -e "s/$OLD_DATE/$NEW_DATE/g" -e "s/$OLD_VERSION/$NEW_VERSION/g" -e "s/$OLD_STATUS/$NEW_STATUS/g" $f > $f.tmp; cp -f $f.tmp $f; rm -f $f.tmp; done Index: trunk/src/combinatorics/Makefile.am =================================================================== --- trunk/src/combinatorics/Makefile.am (revision 8773) +++ trunk/src/combinatorics/Makefile.am (revision 8774) @@ -1,266 +1,266 @@ ## Makefile.am -- Makefile for 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. # ######################################################################## ## The files in this directory implement standard algorithms for WHIZARD ## We create a library which is still to be combined with auxiliary libs. noinst_LTLIBRARIES = libcombinatorics.la check_LTLIBRARIES = libcombinatorics_ut.la COMMON_F90 = \ bytes.f90 \ hashes.f90 \ md5.f90 \ permutations.f90 \ sorting.f90 \ solver.f90 MPI_F90 = \ grids.f90_mpi \ grids_sub.f90_mpi SERIAL_F90 = \ grids.f90_serial \ grids_sub.f90_serial COMMON_SUBMODULES = \ bytes_sub.f90 \ hashes_sub.f90 \ md5_sub.f90 \ permutations_sub.f90 \ sorting_sub.f90 \ solver_sub.f90 COMBINATORICS_SUBMODULES = \ $(COMMON_SUBMODULES) \ grids_sub.f90 COMBINATORICS_MODULES = \ $(COMMON_F90) \ grids.f90 EXTRA_DIST = \ $(COMMON_F90) \ $(COMBINATORICS_SUBMODULES) \ $(SERIAL_F90) \ $(MPI_F90) nodist_libcombinatorics_la_SOURCES = \ $(COMBINATORICS_MODULES) \ $(COMBINATORICS_SUBMODULES) DISTCLEANFILES = grids.f90 if FC_USE_MPI grids.f90: grids.f90_mpi -cp -f $< $@ grids_sub.f90: grids_sub.f90_mpi -cp -f $< $@ else grids.f90: grids.f90_serial -cp -f $< $@ grids_sub.f90: grids_sub.f90_serial -cp -f $< $@ endif libcombinatorics_ut_la_SOURCES = \ md5_uti.f90 md5_ut.f90 \ sorting_uti.f90 sorting_ut.f90 \ grids_uti.f90 grids_ut.f90 \ solver_uti.f90 solver_ut.f90 ## Omitting this would exclude it from the distribution dist_noinst_DATA = combinatorics.nw # Modules and installation # Dump module names into file Modules execmoddir = $(fmoddir)/whizard nodist_execmod_HEADERS = \ ${COMBINAROTICS_MODULES:.f90=.$(FCMOD)} -#Submodules must not be included here +# Submodules must not be included here libcombinatorics_Modules = \ ${COMBINATORICS_MODULES:.f90=} \ ${libcombinatorics_ut_la_SOURCES:.f90=} Modules: Makefile @for module in \ $(libcombinatorics_Modules); do \ echo $$module >> $@.new; \ done @if diff $@ $@.new -q >/dev/null; then \ rm $@.new; \ else \ mv $@.new $@; echo "Modules updated"; \ fi BUILT_SOURCES = Modules ## Fortran module dependencies # Get module lists from other directories module_lists = \ ../basics/Modules \ ../testing/Modules \ ../system/Modules \ ../utilities/Modules $(module_lists): $(MAKE) -C `dirname $@` Modules Module_dependencies.sed: $(nodist_libcombinatorics_la_SOURCES) \ $(libcombinatorics_ut_la_SOURCES) Module_dependencies.sed: $(module_lists) @rm -f $@ echo 's/, *only:.*//' >> $@ echo 's/, *&//' >> $@ echo 's/, *.*=>.*//' >> $@ echo 's/$$/.lo/' >> $@ for list in $(module_lists); do \ dir="`dirname $$list`"; \ for mod in `cat $$list`; do \ echo 's!: '$$mod'.lo$$!': $$dir/$$mod'.lo!' >> $@; \ done \ done DISTCLEANFILES += Module_dependencies.sed # The following line just says # include Makefile.depend # but in a portable fashion (depending on automake's AM_MAKE_INCLUDE @am__include@ @am__quote@Makefile.depend@am__quote@ Makefile.depend: Module_dependencies.sed Makefile.depend: $(nodist_libcombinatorics_la_SOURCES) \ $(libcombinatorics_ut_la_SOURCES) @rm -f $@ for src in $^; do \ module="`basename $$src | sed 's/\.f[90][0358]//'`"; \ grep '^ *use ' $$src \ | grep -v '!NODEP!' \ | sed -e 's/^ *use */'$$module'.lo: /' \ -f Module_dependencies.sed; \ done > $@ DISTCLEANFILES += Makefile.depend # Fortran90 module files are generated at the same time as object files .lo.$(FCMOD): @: # touch $@ AM_FCFLAGS = -I../basics -I../testing -I../system -I../utilities ######################################################################## # For the moment, the submodule dependencies will be hard-coded bytes_sub.lo: bytes.lo hashes_sub.lo: hashes.lo md5_sub.lo: md5.lo permutations_sub.lo: permutations.lo sorting_sub.lo: sorting.lo solver_sub.lo: solver.lo grids_sub.lo: grids.lo ######################################################################## ## Default Fortran compiler options ## Profiling if FC_USE_PROFILING AM_FCFLAGS += $(FCFLAGS_PROFILING) endif ## OpenMP if FC_USE_OPENMP AM_FCFLAGS += $(FCFLAGS_OPENMP) endif # MPI if FC_USE_MPI AM_FCFLAGS += $(FCFLAGS_MPI) endif ######################################################################## ## Non-standard targets and dependencies ## (Re)create F90 sources from NOWEB source. if NOWEB_AVAILABLE FILTER = -filter "sed 's/defn MPI:/defn/'" PRELUDE = $(top_srcdir)/src/noweb-frame/whizard-prelude.nw POSTLUDE = $(top_srcdir)/src/noweb-frame/whizard-postlude.nw combinatorics.stamp: $(PRELUDE) $(srcdir)/combinatorics.nw $(POSTLUDE) @rm -f combinatorics.tmp @touch combinatorics.tmp for src in $(COMMON_F90) $(libcombinatorics_ut_la_SOURCES); do \ $(NOTANGLE) -R[[$$src]] $^ | $(CPIF) $$src; \ done for src in $(COMMON_SUBMODULES); do \ $(NOTANGLE) -R[[$$src]] $^ | $(CPIF) $$src; \ done for src in $(SERIAL_F90:.f90_serial=.f90); do \ $(NOTANGLE) -R[[$$src]] $^ | $(CPIF) $$src'_serial'; \ done for src in $(MPI_F90:.f90_mpi=.f90); do \ $(NOTANGLE) -R[[$$src]] $(FILTER) $^ | $(CPIF) $$src'_mpi'; \ done @mv -f combinatorics.tmp combinatorics.stamp $(COMMON_F90) $(COMMON_SUBMODULES) $(SERIAL_F90) $(MPI_F90) $(libcombinatorics_ut_la_SOURCES): combinatorics.stamp ## Recover from the removal of $@ @if test -f $@; then :; else \ rm -f combinatorics.stamp; \ $(MAKE) $(AM_MAKEFLAGS) combinatorics.stamp; \ fi endif ######################################################################## ## Non-standard cleanup tasks ## Remove sources that can be recreated using NOWEB if NOWEB_AVAILABLE maintainer-clean-noweb: -rm -f *.f90 *.f90_serial *.f90_mpi *.c endif .PHONY: maintainer-clean-noweb ## Remove those sources also if builddir and srcdir are different if NOWEB_AVAILABLE clean-noweb: test "$(srcdir)" != "." && rm -f *.f90 *.f90_serial *.f90_mpi *.c || true endif .PHONY: clean-noweb ## Remove F90 module files clean-local: clean-noweb -rm -f combinatorics.stamp combinatorics.tmp -rm -f *.$(FCMOD) if FC_SUBMODULES -rm -f *.smod *.sub endif ## Remove backup files maintainer-clean-backup: -rm -f *~ .PHONY: maintainer-clean-backup ## Register additional clean targets maintainer-clean-local: maintainer-clean-noweb maintainer-clean-backup Index: trunk/src/pdf_builtin/Makefile.am =================================================================== --- trunk/src/pdf_builtin/Makefile.am (revision 8773) +++ trunk/src/pdf_builtin/Makefile.am (revision 8774) @@ -1,111 +1,122 @@ ## Makefile.am -- Makefile for 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. # ######################################################################## ## The files in this directory end up in an auxiliary libtool library. noinst_LTLIBRARIES = libpdf_builtin.la +PDF_MODULES = \ + cteq6pdf.f90 mrst2004qed.f90 mstwpdf.f90 ct10pdf.f90 \ + CJpdf.f90 ct14pdf.f90 ct18pdf.f90 pdf_builtin.f90 + +PDF_SUBMODULES = \ + pdf_builtin_sub.f90 + libpdf_builtin_la_SOURCES = \ - cteq6pdf.f90 pdf_builtin.f90 mrst2004qed.f90 mstwpdf.f90 \ - ct10pdf.f90 CJpdf.f90 ct14pdf.f90 ct18pdf.f90 + $(PDF_MODULES) \ + $(PDF_SUBMODULES) $(libpdf_builtin_la_OBJECTS): \ ../basics/kinds.$(FCMOD) \ ../basics/iso_varying_string.$(FCMOD) \ ../basics/io_units.$(FCMOD) \ ../utilities/format_utils.$(FCMOD) \ ../system/diagnostics.$(FCMOD) # Modules and installation # Dump module names into file Modules execmoddir = $(fmoddir)/whizard nodist_execmod_HEADERS = \ pdf_builtin.$(FCMOD) \ cteq6pdf.$(FCMOD) \ cj_pdf.$(FCMOD) \ mrst2004qed.$(FCMOD) \ mstwpdf.$(FCMOD) \ ct10pdf.$(FCMOD) \ ct14pdf.$(FCMOD) \ ct18pdf.$(FCMOD) SUFFIXES = .lo .$(FCMOD) # Fortran90 module files are generated at the same time as object files .lo.$(FCMOD): @: # touch $@ AM_FCFLAGS = -I../basics -I../utilities -I../system ######################################################################## +# For the moment, the submodule dependencies will be hard-coded +pdf_builtin_sub.lo: pdf_builtin.lo + +######################################################################## ## Default Fortran compiler options ## Profiling if FC_USE_PROFILING AM_FCFLAGS += $(FCFLAGS_PROFILING) endif ## OpenMP if FC_USE_OPENMP AM_FCFLAGS += $(FCFLAGS_OPENMP) endif # MPI if FC_USE_MPI AM_FCFLAGS += $(FCFLAGS_MPI) endif ######################################################################## ## Explicit dependencies pdf_builtin.lo: cteq6pdf.lo mrst2004qed.lo mstwpdf.lo ct10pdf.lo CJpdf.lo ct14pdf.lo ct18pdf.lo pdf_builtin.$(FCMOD): pdf_builtin.lo cteq6pdf.$(FCMOD): cteq6pdf.lo cj_pdf.$(FCMOD): CJpdf.lo mrst2004qed.$(FCMOD): mrst2004qed.lo mstwpdf.$(FCMOD): mstwpdf.lo ct10pdf.$(FCMOD): ct10pdf.lo ct14pdf.$(FCMOD): ct14pdf.lo ct18pdf.$(FCMOD): ct18pdf.lo ######################################################################## ## Non-standard cleanup tasks ## Remove F90 module files clean-local: -rm -f *.$(FCMOD) if FC_SUBMODULES - -rm -f *.smod + -rm -f *.smod *.sub endif ## Remove backup files maintainer-clean-local: -rm -f *~ Index: trunk/src/pdf_builtin/pdf_builtin.f90 =================================================================== --- trunk/src/pdf_builtin/pdf_builtin.f90 (revision 8773) +++ trunk/src/pdf_builtin/pdf_builtin.f90 (revision 8774) @@ -1,1065 +1,92 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! 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. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Wrap a common interface around the different PDF sets. module pdf_builtin use kinds, only: default, double use iso_varying_string, string_t => varying_string - use string_utils, only: upper_case - use constants, only: PI - use diagnostics - use mrst2004qed - use cteq6pdf - use mstwpdf - use ct10pdf - use CJ_pdf - use ct14pdf - use ct18pdf implicit none save private -! The available sets - integer, parameter :: nsets = 24 - integer, parameter, public :: & - CTEQ6M = 1, CTEQ6D = 2, CTEQ6L = 3, CTEQ6L1 = 4, & - MRST2004QEDp = 5, MRST2004QEDn = 6, MSTW2008LO = 7, MSTW2008NLO = 8, & - MSTW2008NNLO = 9, CT10 = 10, CJ12_max = 11, CJ12_mid = 12, & - CJ12_min = 13, MMHT2014LO = 14, MMHT2014NLO = 15, MMHT2014NNLO = 16, & - CT14LL = 17, CT14L = 18, CT14N = 19, CT14NN = 20, CJ15_LO = 21, & - CJ15_NLO = 22, CT18N = 23, CT18NN = 24 - -! Limits - real(kind=default), parameter :: & - cteq6_q_min = 1.3, cteq6_q_max = 10.E3, & - cteq6_x_min = 1.E-6, cteq6_x_max = 1., & - mrst2004qed_q_min = sqrt (1.26_default), mrst2004qed_q_max = sqrt (0.99E7_default), & - mrst2004qed_x_min = 1.01E-5, mrst2004qed_x_max = 1., & - mstw2008_q_min = 1.01, mstw2008_q_max = sqrt (0.99E9_default), & - mstw2008_x_min = 1.01E-6, mstw2008_x_max = 1., & - ct10_q_min = 1.3, ct10_q_max = 1.E5, & - ct10_x_min = 1.E-8, ct10_x_max = 1., & - cj12_q_min = 1.3, cj12_q_max = 1.E5, & - cj12_x_min = 1.E-6, cj12_x_max = 1., & - mmht2014_q_min = mstw2008_q_min, mmht2014_q_max = mstw2008_q_max, & - mmht2014_x_min = mstw2008_x_min, mmht2014_x_max = mstw2008_x_max, & - ct14_q_min = 1.3, ct14_q_max = 1.E5, ct14_x_min = 1.E-9, & - ct14_x_max = 1., ct18_q_min = 1.3, ct18_q_max = 1.E5, & - ct18_x_min = 1.E-9, ct18_x_max = 1. - -! Lambda_QCD and quark masses - - real(kind=double), dimension(6), parameter, public :: alam_ct10 = & - [ 0.37219423568859566_double, 0.37219423568859566_double, & - 0.37219423568859566_double, 0.37219423568859566_double, & - 0.32560730624033124_double, 0.22600000000000001_double ] - - real(kind=double), dimension(6), parameter, public :: alam_cteq6m = & - [ 0.37248648555333957_double, 0.37248648555333957_double, & - 0.37248648555333957_double, 0.37248648555333957_double, & - 0.32590307925177625_double, 0.22623045868581182_double ] - - real(kind=double), dimension(6), parameter, public :: alam_cteq6l = & - [ 0.37218603304249098_double, 0.37218603304249098_double, & - 0.37218603304249098_double, 0.37218603304249098_double, & - 0.32559900534407232_double, 0.22599353258372407_double ] - - real(kind=double), dimension(6), parameter, public :: alam_cteq6ll = & - [ 0.24560434095679567_double, 0.24560434095679567_double, & - 0.24560434095679567_double, 0.24560434095679567_double, & - 0.21495099184302788_double, 0.16499885346541945_double ] - - real(kind=double), dimension(0:5), parameter :: amhat = & - [ 0.0000000000000000_double, 0.0000000000000000_double, & - 0.0000000000000000_double, 0.0000000000000000_double, & - 1.3000000000000000_double, 4.500000000000000_double ] - - !!! Global parameter for minimal value of evolution - real(kind=double), parameter, private :: amn = 0.375_double, & - tol = .0005_double - - integer, parameter, private :: nhq = 2 - - real(kind=double), parameter :: & - ca = 3.D0, cf = 4./3.D0, tr = 0.5D0 - real(kind=double), parameter :: & - b00 = 11./3.d0 * ca, b01 = -4./3.d0 * tr, b10 = 34./3.d0 * ca**2, & - b11 = -20./3.d0 * ca*tr - 4.* cf*tr - - ! Init flags - integer :: cteq6_initialized = -1 - integer :: mstw2008_initialized = -1 - integer :: mmht2014_initialized = -1 - logical :: ct10_initialized = .false. - integer :: cj12_initialized = -1 - integer :: ct14_initialized = -1 - integer :: cj15_initialized = -1 - integer :: ct18_initialized = -1 - logical :: & - mrst2004qedp_initialized = .false., & - mrst2004qedn_initialized = .false. - type(string_t) :: mrst2004qedp_prefix, mrst2004qedn_prefix, & - mstw2008_prefix, mmht2014_prefix - -! Public stuff public :: pdf_init public :: pdf_get_name public :: pdf_evolve public :: pdf_evolve_LHAPDF public :: pdf_provides_photon public :: pdf_get_id public :: pdf_alphas public :: pdf_alphas_LHAPDF public :: pdf_getmass -contains - -! Get PDF name - function pdf_get_name (pdftype) result (name) - integer, intent(in) :: pdftype - type(string_t) :: name - select case (pdftype) - case (CTEQ6M) - name = var_str ("CTEQ6M") - case (CTEQ6D) - name = var_str ("CTEQ6D") - case (CTEQ6L) - name = var_str ("CTEQ6L") - case (CTEQ6L1) - name = var_str ("CTEQ6L1") - case (MRST2004QEDp) - name = var_str ("MRST2004QEDp") - case (MRST2004QEDn) - name = var_str ("MRST2004QEDn") - case (MSTW2008LO) - name = var_str ("MSTW2008LO") - case (MSTW2008NLO) - name = var_str ("MSTW2008NLO") - case (MSTW2008NNLO) - name = var_str ("MSTW2008NNLO") - case (CT10) - name = var_str ("CT10") - case (CJ12_max) - name = var_str ("CJ12_max") - case (CJ12_mid) - name = var_str ("CJ12_mid") - case (CJ12_min) - name = var_str ("CJ12_min") - case (CJ15_LO) - name = var_str ("CJ15LO") - case (CJ15_NLO) - name = var_str ("CJ15NLO") - case (MMHT2014LO) - name = var_str ("MMHT2014LO") - case (MMHT2014NLO) - name = var_str ("MMHT2014NLO") - case (MMHT2014NNLO) - name = var_str ("MMHT2014NNLO") - case (CT14LL) - name = var_str ("CT14LL") - case (CT14L) - name = var_str ("CT14L") - case (CT14N) - name = var_str ("CT14N") - case (CT14NN) - name = var_str ("CT14NN") - case (CT18N) - name = var_str ("CT18N") - case (CT18NN) - name = var_str ("CT18NN") - case default - call msg_fatal ("pdf_builtin: internal: invalid PDF set!") - end select - end function pdf_get_name - -! Get the ID of a PDF set - function pdf_get_id (name) result (id) - type(string_t), intent(in) :: name - integer :: id - do id = 1, nsets - if (upper_case (pdf_get_name (id)) == upper_case (name)) return - end do - id = -1 - end function pdf_get_id - -! Query whether a PDF supplies a photon distribution - function pdf_provides_photon (pdftype) result (flag) - integer, intent(in) :: pdftype - logical :: flag - select case (pdftype) - case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) - flag = .false. - case (MRST2004QEDp, MRST2004QEDn) - flag = .true. - case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO) - flag = .false. - case (CT10) - flag = .false. - case (CJ12_max, CJ12_mid, CJ12_min) - flag = .false. - case (CJ15_LO, CJ15_NLO) - flag = .false. - case (MMHT2014LO, MMHT2014NLO, MMHT2014NNLO) - flag = .false. - case (CT14LL, CT14L, CT14N, CT14NN) - flag = .false. - case (CT18N, CT18NN) - flag = .false. - case default - call msg_fatal ("pdf_builtin: internal: invalid PDF set!") - end select - end function pdf_provides_photon - -! Initialize a PDF - subroutine pdf_init (pdftype, prefix, verbose) - integer, intent(in) :: pdftype - type(string_t), intent(in), optional :: prefix - type(string_t) :: mprefix - logical, intent(in), optional :: verbose - logical :: mverbose - if (present (prefix)) then - mprefix = prefix - else - mprefix = "" - end if - if (present (verbose)) then - mverbose = verbose - else - mverbose = .true. - end if - select case (pdftype) - case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) - if (cteq6_initialized == pdftype) return - call setctq6 (pdftype, char (mprefix)) - cteq6_initialized = pdftype - case (MRST2004QEDp) - if (mrst2004qedp_initialized) return - mrst2004qedp_initialized = .true. - mrst2004qedp_prefix = mprefix - case (MRST2004QEDn) - if (mrst2004qedn_initialized) return - mrst2004qedn_initialized = .true. - mrst2004qedn_prefix = mprefix - case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO) - if (mstw2008_initialized == pdftype) return - mstw2008_initialized = pdftype - mstw2008_prefix = mprefix - case (CT10) - if (ct10_initialized) return - call setct10 (char (mprefix), 100) - ct10_initialized = .true. - case (CJ12_max) - if (cj12_initialized == pdftype) return - call setCJ (char (mprefix), 300) - cj12_initialized = pdftype - case (CJ12_mid) - if (cj12_initialized == pdftype) return - call setCJ (char (mprefix), 200) - cj12_initialized = pdftype - case (CJ12_min) - if (cj12_initialized == pdftype) return - call setCJ (char (mprefix), 100) - cj12_initialized = pdftype - case (CJ15_LO) - if (cj15_initialized == pdftype) return - call setCJ (char (mprefix), 400) - cj15_initialized = pdftype - case (CJ15_NLO) - if (cj15_initialized == pdftype) return - call setCJ (char (mprefix), 500) - cj15_initialized = pdftype - case (MMHT2014LO, MMHT2014NLO, MMHT2014NNLO) - if (mmht2014_initialized == pdftype) return - mmht2014_initialized = pdftype - mmht2014_prefix = mprefix - case (CT14LL) - if (ct14_initialized == pdftype) return - ct14_initialized = pdftype - call setct14 (char (mprefix), 1) - case (CT14L) - if (ct14_initialized == pdftype) return - ct14_initialized = pdftype - call setct14 (char (mprefix), 2) - case (CT14N) - if (ct14_initialized == pdftype) return - ct14_initialized = pdftype - call setct14 (char (mprefix), 3) - case (CT14NN) - if (ct14_initialized == pdftype) return - ct14_initialized = pdftype - call setct14 (char (mprefix), 4) - case (CT18N) - if (ct18_initialized == pdftype) return - ct18_initialized = pdftype - call setct18 (char (prefix), 1) - case (CT18NN) - if (ct18_initialized == pdftype) return - ct18_initialized = pdftype - call setct18 (char (prefix), 2) - case default - call msg_fatal ("pdf_builtin: internal: invalid PDF set!") - end select - if (mverbose) call msg_message ("Initialized builtin PDF " // & - char (pdf_get_name (pdftype))) - !!! Up to now only proton - end subroutine pdf_init - -! Evolve PDF - subroutine pdf_evolve (pdftype, x, q, f, fphoton) - integer, intent(in) :: pdftype - real(kind=default), intent(in) :: x, q - real(kind=default), intent(out), optional :: f(-6:6), fphoton - real(kind=double) :: mx, mq - real(kind=double) :: upv, dnv, ups, dns, str, chm, bot, glu, phot, & - sbar, bbar, cbar - type(string_t) :: setname - select case (pdftype) - case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) - if (cteq6_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (cteq6_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (cteq6_initialized)) // & - " cannot be used simultaneously") - mx = max (min (x, cteq6_x_max), cteq6_x_min) - mq = max (min (q, cteq6_q_max), cteq6_q_min) - if (present (f)) f = [ 0._double, & - ctq6pdf (-5, mx, mq), ctq6pdf (-4, mx, mq), & - ctq6pdf (-3, mx, mq), ctq6pdf (-1, mx, mq), & - ctq6pdf (-2, mx, mq), ctq6pdf ( 0, mx, mq), & - ctq6pdf ( 2, mx, mq), ctq6pdf ( 1, mx, mq), & - ctq6pdf ( 3, mx, mq), ctq6pdf ( 4, mx, mq), & - ctq6pdf ( 5, mx, mq), 0._double ] - if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & - char (pdf_get_name (pdftype)) // " which does not provide it!") - case (MRST2004QEDp) - if (.not. mrst2004qedp_initialized) call msg_fatal ( & - "pdf_builtin: internal: PDF set MRST2004QEDp requested without " // & - "initialization!") - mx = max (min (x, mrst2004qed_x_max), mrst2004qed_x_min) - mq = max (min (q, mrst2004qed_q_max), mrst2004qed_q_min) - call mrstqed (mx, mq, 1, upv, dnv, ups, dns, str, chm, bot, glu, phot, & - char (mrst2004qedp_prefix)) - if (present (f)) f = & - [ 0._double, bot, chm, str, ups, dns, glu, dns + dnv, ups + upv, & - str, chm, bot, 0._double ] / mx - if (present (fphoton)) fphoton = phot / mx - case (MRST2004QEDn) - if (.not. mrst2004qedn_initialized) call msg_fatal ( & - "pdf_builtin: internal: PDF set MRST2004QEDn requested without " // & - "initialization!") - mx = max (min (x, mrst2004qed_x_max), mrst2004qed_x_min) - mq = max (min (q, mrst2004qed_q_max), mrst2004qed_q_min) - call mrstqed (mx, mq, 2, upv, dnv, ups, dns, str, chm, bot, glu, phot, & - char (mrst2004qedn_prefix)) - if (present (f)) f = & - [ 0._double, bot, chm, str, ups, dns, glu, dns + dnv, ups + upv, & - str, chm, bot, 0._double ] / mx - if (present (fphoton)) fphoton = phot / mx - case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO) - if (mstw2008_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (mstw2008_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (mstw2008_initialized)) // & - " cannot be used simultaneously") - select case (pdftype) - case (MSTW2008LO) - setname = var_str ("mstw2008lo") - case (MSTW2008NLO) - setname = var_str ("mstw2008nlo") - case (MSTW2008NNLO) - setname = var_str ("mstw2008nnlo") - end select - mx = max (min (x, mstw2008_x_max), mstw2008_x_min) - mq = max (min (q, mstw2008_q_max), mstw2008_q_min) - call getmstwpdf (char (mstw2008_prefix), char (setname), 0, mx, mq, & - upv, dnv, ups, dns, str, sbar, chm, cbar, bot, bbar, glu, phot) - if (present (f)) f = & - [ 0._double, bbar, cbar, sbar, ups, dns, glu, dns + dnv, ups + upv, & - str, chm, bot, 0._double ] / mx - if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & - char (pdf_get_name (pdftype)) // " which does not provide it!") - case (CT10) - if (.not. ct10_initialized) call msg_fatal ( & - "pdf_builtin: internal: PDF set CT10 requested without " // & - "initialization!") - mx = max (min (x, ct10_x_max), ct10_x_min) - mq = max (min (q, ct10_q_max), ct10_q_min) - if (present (f)) f = [ 0._double, & - getct10pdf (-5, mx, mq), getct10pdf (-4, mx, mq), & - getct10pdf (-3, mx, mq), getct10pdf (-1, mx, mq), & - getct10pdf (-2, mx, mq), getct10pdf ( 0, mx, mq), & - getct10pdf ( 2, mx, mq), getct10pdf ( 1, mx, mq), & - getct10pdf ( 3, mx, mq), getct10pdf ( 4, mx, mq), & - getct10pdf ( 5, mx, mq), 0._double ] - if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & - "CT10 which does not provide it!") - case (CJ12_max, CJ12_mid, CJ12_min) - if (cj12_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (cj12_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (cj12_initialized)) // & - " cannot be used simultaneously") - mx = max (min (x, cj12_x_max), cj12_x_min) - mq = max (min (q, cj12_q_max), cj12_q_min) - if (present (f)) f = [ 0._double, & - CJpdf (-5, mx, mq), CJpdf (-4, mx, mq), CJpdf (-3, mx, mq), & - CJpdf (-2, mx, mq), CJpdf (-1, mx, mq), CJpdf ( 0, mx, mq), & - CJpdf ( 1, mx, mq), CJpdf ( 2, mx, mq), CJpdf ( 3, mx, mq), & - CJpdf ( 4, mx, mq), CJpdf ( 5, mx, mq), 0._double ] - if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & - "CJ12 which does not provide it!") - case (CJ15_LO, CJ15_NLO) - if (cj15_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (cj15_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (cj15_initialized)) // & - " cannot be used simultaneously") - mx = max (min (x, cj12_x_max), cj12_x_min) - mq = max (min (q, cj12_q_max), cj12_q_min) - if (present (f)) f = [ 0._double, & - CJpdf (-5, mx, mq), CJpdf (-4, mx, mq), CJpdf (-3, mx, mq), & - CJpdf (-2, mx, mq), CJpdf (-1, mx, mq), CJpdf ( 0, mx, mq), & - CJpdf ( 1, mx, mq), CJpdf ( 2, mx, mq), CJpdf ( 3, mx, mq), & - CJpdf ( 4, mx, mq), CJpdf ( 5, mx, mq), 0._double ] - if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & - "CJ15 which does not provide it!") - case (MMHT2014LO, MMHT2014NLO, MMHT2014NNLO) - if (mmht2014_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (mmht2014_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (mmht2014_initialized)) // & - " cannot be used simultaneously") - select case (pdftype) - case (MMHT2014LO) - setname = var_str ("mmht2014lo") - case (MMHT2014NLO) - setname = var_str ("mmht2014nlo") - case (MMHT2014NNLO) - setname = var_str ("mmht2014nnlo") - end select - mx = max (min (x, mmht2014_x_max), mmht2014_x_min) - mq = max (min (q, mmht2014_q_max), mmht2014_q_min) - call getmstwpdf (char (mmht2014_prefix), char (setname), 0, mx, mq, & - upv, dnv, ups, dns, str, sbar, chm, cbar, bot, bbar, glu, phot) - if (present (f)) f = & - [ 0._double, bbar, cbar, sbar, ups, dns, glu, dns + dnv, ups + upv, & - str, chm, bot, 0._double ] / mx - if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & - char (pdf_get_name (pdftype)) // " which does not provide it!") - case (CT14LL, CT14L, CT14N, CT14NN) - if (ct14_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (ct14_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (ct14_initialized)) // & - " cannot be used simultaneously") - mx = max (min (x, ct14_x_max), ct14_x_min) - mq = max (min (q, ct14_q_max), ct14_q_min) - if (present (f)) f = [ 0._double, & - getct14pdf (-5, mx, mq), getct14pdf (-4, mx, mq), & - getct14pdf (-3, mx, mq), getct14pdf (-1, mx, mq), & - getct14pdf (-2, mx, mq), getct14pdf ( 0, mx, mq), & - getct14pdf ( 2, mx, mq), getct14pdf ( 1, mx, mq), & - getct14pdf ( 3, mx, mq), getct14pdf ( 4, mx, mq), & - getct14pdf ( 5, mx, mq), 0._double ] - case (CT18N, CT18NN) - if (ct18_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (ct18_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (ct18_initialized)) // & - " cannot be used simultaneously") - mx = max (min (x, ct18_x_max), ct18_x_min) - mq = max (min (q, ct18_q_max), ct18_q_min) - if (present (f)) f = [ 0._double, & - getct18pdf (-5, mx, mq), getct18pdf (-4, mx, mq), & - getct18pdf (-3, mx, mq), getct18pdf (-1, mx, mq), & - getct18pdf (-2, mx, mq), getct18pdf ( 0, mx, mq), & - getct18pdf ( 2, mx, mq), getct18pdf ( 1, mx, mq), & - getct18pdf ( 3, mx, mq), getct18pdf ( 4, mx, mq), & - getct18pdf ( 5, mx, mq), 0._double ] - if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & - char (pdf_get_name (pdftype)) // " which does not provide it!") - case default - call msg_fatal ("pdf_builtin: internal: invalid PDF set!") - end select - end subroutine pdf_evolve - -! included for compatibility with LHAPDF -! use a double precision array for the pdfs instead of a -! default precision - subroutine pdf_evolve_LHAPDF (set, x, q, ff) - integer, intent(in) :: set - double precision, intent(in) :: x, q - real(kind=default) :: dx, dq - double precision, dimension(-6:6), intent(out) :: ff - - real(kind=default) :: f(-6:6) - dx = x - dq = q - - call pdf_evolve (set, dx, dq, f) - ff = f - end subroutine pdf_evolve_LHAPDF - -! PDF-specific running alphas - function pdf_alphas (pdftype, q) result (alphas) - real(kind=double) :: as - real(kind=default), intent(in) :: q - real(kind=double) :: qdummy - real(kind=default) :: alphas - integer, intent(in) :: pdftype - integer, parameter :: nset_dummy = 1 - qdummy = q - select case (pdftype) - case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) - if (cteq6_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (cteq6_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (cteq6_initialized)) // & - " cannot be used simultaneously") - select case (pdftype) - case (CTEQ6M, CTEQ6D) - as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6m) - case (CTEQ6L) - as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6l) - case (CTEQ6L1) - as = PI*pdf_builtin_alpi (qdummy,1,5,alam_cteq6ll) - case default - call msg_fatal ("Unrecognized PDF setup in evolution of alphas.") - end select - case (CT10) - if (.not. ct10_initialized) call msg_fatal ( & - "pdf_builtin: internal: PDF set CT10 requested without " // & - "initialization!") - as = PI*pdf_builtin_alpi (qdummy,2,5,alam_ct10) - case (MRST2004QEDp) - if (.not. mrst2004qedp_initialized) call msg_fatal ( & - "pdf_builtin: internal: PDF set MRST2004QEDp requested without " // & - "initialization!") - call mrst_alphas (as,qdummy) - case (MRST2004QEDn) - if (.not. mrst2004qedn_initialized) call msg_fatal ( & - "pdf_builtin: internal: PDF set MRST2004QEDn requested without " // & - "initialization!") - call mrst_alphas (as,qdummy) - case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO) - if (mstw2008_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (mstw2008_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (mstw2008_initialized)) // & - " cannot be used simultaneously") - select case (pdftype) - case (MSTW2008lo) - call pdf_builtin_alfa (as,qdummy,0) - case (MSTW2008nlo) - call pdf_builtin_alfa (as,qdummy,1) - case (MSTW2008nnlo) - call pdf_builtin_alfa (as,qdummy,2) - end select - case (CJ12_max, CJ12_mid, CJ12_min) - if (cj12_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (cj12_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (cj12_initialized)) // & - " cannot be used simultaneously") - as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6m) - case (CJ15_LO, CJ15_NLO) - if (cj15_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (cj15_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (cj15_initialized)) // & - " cannot be used simultaneously") - select case (pdftype) - case (CJ15_NLO) - as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6m) - case (CJ15_LO) - as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6l) - end select - case (MMHT2014LO, MMHT2014NLO, MMHT2014NNLO) - if (mmht2014_initialized < 0) & - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (mmht2014_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (mmht2014_initialized)) // & - " cannot be used simultaneously") - select case (pdftype) - case (MMHT2014lo) - call pdf_builtin_alfa (as,qdummy,0) - case (MMHT2014nlo) - call pdf_builtin_alfa (as,qdummy,1) - case (MMHT2014nnlo) - call pdf_builtin_alfa (as,qdummy,2) - end select - case (CT14LL, CT14L, CT14N, CT14NN) - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (ct14_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (ct14_initialized)) // & - " cannot be used simultaneously") - as = CT14Alphas(qdummy) - case (CT18N, CT18NN) - call msg_fatal ("pdf_builtin: internal: PDF set " // & - char (pdf_get_name (pdftype)) // " requested without initialization!") - if (ct18_initialized /= pdftype) & - call msg_fatal ( & - "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & - char (pdf_get_name (ct18_initialized)) // & - " cannot be used simultaneously") - as = CT18Alphas(qdummy) - case default - call msg_fatal ("pdf_builtin: internal: invalid PDF set!") - end select - alphas = as - end function pdf_alphas - - function pdf_alphas_LHAPDF (pdftype, q) result (alphas) - integer, intent(in) :: pdftype - real(kind=double), intent(in) :: q - real(kind=double) :: alphas - real(kind=default) :: q_def - q_def = q - alphas = pdf_alphas (pdftype, q_def) - end function pdf_alphas_LHAPDF - - function pdf_getmass (nf) result (mass) - integer, intent(in) :: nf - real(kind=double) :: mass - select case (abs (nf)) - case (1:3) - mass = 0._double - case (4) - mass = 1.3_double - case (5) - mass = 4.5_double - case default - call msg_fatal ("PDF builtin: invalid PDG code for quark mass.") - end select - end function pdf_getmass - - function pdf_builtin_alpi (AMU,NORDER,NF,alam) - integer, intent(in) :: nf, norder - real(kind=double), intent(in) :: amu - real(kind=double), dimension(6), intent(in) :: alam - real(kind=double) :: alm - real(kind=double) :: pdf_builtin_alpi - integer :: irt, n_eff - n_eff = num_flavor (AMU,NF) - alm = alam (n_eff+1) - pdf_builtin_alpi = pdf_builtin_alpqcd (norder, n_eff, amu/alm, irt) - select case (irt) - case (1) - call msg_warning ("AMU < ALAM in pdf_builtin_alpi") - case (2) - call msg_warning ("pdf_builtin_alpi > 3; Be aware!") - end select - end function pdf_builtin_alpi - - double precision function pdf_builtin_ALPQCD (IRDR, NF, RML, IRT) - real(kind=double) :: b0, b1 - integer, intent(in) :: irdr, nf - integer :: irt - real(kind=double), intent(in) :: rml - real(kind=double) :: al, aln - real(kind=double) :: rm2 - IRT = 0 - if (IRDR .LT. 1 .OR. IRDR .GT. 2) THEN - print *, & - & 'Order out of range in pdf_builtin_ALPQCD: IRDR = ', IRDR - STOP - end if - b0 = (11.d0*CA - 2.* NF) / 3.d0 - b1 = (34.d0*CA**2 - 10.d0*CA*NF - 6.d0*CF*NF) / 3.d0 - rm2 = rml**2 - if (RM2 .LE. 1.) THEN - IRT = 1 - pdf_builtin_ALPQCD = 99. - RETURN - end if - aln = log (RM2) - al = 4.d0/ B0 / aln - IF (IRDR .GE. 2) AL = AL * (1.d0-B1*log(ALN) / ALN / B0**2) - if (AL .GE. 3.) THEN - IRT = 2 - end if - pdf_builtin_ALPQCD = AL - end function pdf_builtin_ALPQCD - - function num_flavor (amu,nf) - real(kind=double), intent(in) :: amu - integer, intent(in) :: nf - integer :: num_flavor, i - num_flavor = nf - nhq - IF ((num_flavor .EQ. nf) .OR. (amu .LE. amn)) GOTO 20 - do 10 I = nf - NHQ + 1, nf - if (amu .GE. amhat(I)) THEN - num_flavor = I - else - goto 20 - end if -10 continue -20 return - end function num_flavor - - subroutine pdf_builtin_alfa (alfas,Qalfa,norder) - integer, intent(in) :: norder - real(kind=double), intent(in) :: Qalfa - real(kind=double), intent(out) :: alfas - real(kind=double) :: alphasq0 - real(kind=double) :: mCharmSave, mBottomSave, mTopSave - mTopSave = 10000000000.0_double - mBottomSave = 4.75_double - mCharmSave = 1.4_double - select case (norder) - case (0) - alphasq0 = 0.68183_double - case (1) - alphasq0 = 0.491278_double - case (2) - alphasq0 = 0.45077_double - case default - call msg_fatal ("Wrong MSTW order!") - end select - if (Qalfa.GT.mTopSave) call msg_fatal & - ("No MSTW alphas for such high scales.") - alfas = pdf_builtin_mstw_alphas (Qalfa,norder,mTopSave**2) - end subroutine pdf_builtin_alfa - - function pdf_builtin_mstw_alphas (mur, norder, m2t) result (alfas) - real(kind=double), intent(in) :: mur - real(kind=double) :: alfas - integer, intent(in) :: norder - integer :: nf - real(kind=double), intent(in) :: m2t - real(kind=double) :: as0, asc, asb - real(kind=double) :: r2, m2, asi, asf, r20, r2b, r2c, r2t - real(kind=double), parameter :: m20 = 1.0_double, m2c = 1.96_double, & - m2b = 22.5625_double - real(kind=double) :: logfr - integer :: nff - real(kind=double) :: ast - nff = 4 - select case (norder) - case (0) - as0 = 5.4258307424173556E-002_double - asc = 4.0838232991033577E-002_double - asb = 2.2297506503539639E-002_double - case (1) - as0 = 3.9094820221093209E-002_double - asc = 3.0202751103489092E-002_double - asb = 1.7751676069301441E-002_double - case (2) - as0 = 3.5871136848766867E-002_double - asc = 2.8092349518054758E-002_double - asb = 1.6936586710596203E-002_double - case default - call msg_fatal ("Wrong evolution order in MSTW alphas evolution.") - end select - r2 = mur**2 - m2 = r2 * exp(+logfr) - if (m2 .gt. m2t) then - nf = 6 - r2t = m2t * r2/m2 - asi = ast - asf = mstw_alphas (r2,r2t,ast,nf,norder) - else - if (m2 .gt. m2b) then - nf = 5 - r2b = m2b * r2/m2 - asi = asb - asf = mstw_alphas (r2,r2b,asb,nf,norder) - else - if (m2 .gt. m2c) then - nf = 4 - r2c = m2c * r2/m2 - asi = asc - asf = mstw_alphas (r2,r2c,asc,nf,norder) - else - nf = 3 - r20 = m20 * r2/m2 - asi = as0 - asf = mstw_alphas (r2,r20,as0,nf,norder) - end if - end if - end if - alfas = 4.D0*PI*ASF - end function pdf_builtin_mstw_alphas - - function mstw_alphas (r2, r20, as0, nf, naord) result (as) - real(kind=double) :: as - real(kind=double), intent(in) :: r2, r20, as0 - integer, intent(in) :: nf, naord - integer, parameter :: nastps = 20 - integer :: k1 - real(kind=double), parameter :: sxth = 1./6. - real(kind=double) :: lrrat, dlr, xk0, xk1, xk2, xk3 - as = as0 - lrrat = log (r2/r20) - dlr = lrrat / nastps - ! ..Solution of the evolution equation depending on NAORD - ! (fourth-order Runge-Kutta beyond the leading order) - select case (naord) - case (0) - as = as0 / (1.+ beta0(nf) * as0 * lrrat) - case (1) - do k1 = 1, nastps - xk0 = dlr * fbeta1 (as, nf) - xk1 = dlr * fbeta1 (as + 0.5 * xk0, nf) - xk2 = dlr * fbeta1 (as + 0.5 * xk1, nf) - xk3 = dlr * fbeta1 (as + xk2, nf) - as = as + sxth * (xk0 + 2.* xk1 + 2.* xk2 + xk3) - end do - case (2) - do k1 = 1, nastps - xk0 = dlr * fbeta2 (as, nf) - xk1 = dlr * fbeta2 (as + 0.5 * xk0, nf) - xk2 = dlr * fbeta2 (as + 0.5 * xk1, nf) - xk3 = dlr * fbeta2 (as + xk2, nf) - as = as + sxth * (xk0 + 2.* xk1 + 2.* xk2 + xk3) - end do - case (3) - do k1 = 1, nastps - xk0 = dlr * fbeta3 (as, nf) - xk1 = dlr * fbeta3 (as + 0.5 * xk0, nf) - xk2 = dlr * fbeta3 (as + 0.5 * xk1, nf) - xk3 = dlr * fbeta3 (as + xk2, nf) - as = as + sxth * (xk0 + 2.* xk1 + 2.* xk2 + xk3) - end do - end select - end function mstw_alphas - - double precision function beta0 (nf) - integer, intent(in) :: nf - beta0 = b00 + b01 * nf - end function beta0 - - double precision function beta1 (nf) - integer, intent(in) :: nf - beta1 = b10 + b11 * nf - end function beta1 - - double precision function beta2 (nf) - integer, intent(in) :: nf - beta2 = 1428.50 - 279.611 * nf + 6.01852 * nf**2 - end function beta2 - - double precision function beta3 (nf) - integer, intent(in) :: nf - beta3 = 29243.0 - 6946.30 * nf + 405.089 * nf**2 + & - 1.49931 * nf**3 - end function beta3 - - ! ..The beta functions FBETAn at N^nLO for n = 1, 2, and 3 - double precision function fbeta1 (a, nf) - real(kind=double), intent(in) :: a - integer, intent(in) :: nf - fbeta1 = - a**2 * ( beta0(nf) + a * beta1(nf) ) - end function fbeta1 - - double precision function fbeta2 (a, nf) - real(kind=double), intent(in) :: a - integer, intent(in) :: nf - fbeta2 = - a**2 * ( beta0(nf) + a * ( beta1(nf) & - + a * beta2(nf) ) ) - end function fbeta2 - - double precision function fbeta3 (a, nf) - real(kind=double), intent(in) :: a - integer, intent(in) :: nf - fbeta3 = - a**2 * ( beta0(nf) + a * ( beta1(nf) & - + a * ( beta2(nf) + a * beta3(nf)) ) ) - end function fbeta3 - - subroutine mrst_alphas (alpha, q) - real(kind=double), intent(in) :: q - real(kind=double), intent(out) :: alpha - integer :: idir - integer, parameter :: nflav = 4 - real(kind=double) :: alphas, qz2, qz, alambda, astep, & - tol2 - alphas = 0.1205_double - qz2=8315._double - qz = sqrt(qz2) - alambda = 0.3000_double - astep = 0.010_double - tol2 = 0.0000001_double - idir = +1 -100 continue - alambda = alambda + idir*astep - call mrst_lambda (nflav,alpha,qz,alambda) - if (idir*(alphas-alpha).gt.0.0) then - goto 200 - else - astep = 0.5*astep - idir = -1*idir - goto 200 - end if -200 continue - if(abs(alpha-alphas).gt.tol2) goto 100 - ! alambda found -- save it !!!! - ! next call mrst_lambda to get alphas at q with the correct alambda - call mrst_lambda (nflav,alpha,q,alambda) - RETURN - end subroutine mrst_alphas - - subroutine mrst_lambda (nflav,alpha,Q,alambda) - integer, intent(in) :: nflav - integer :: ith - real(kind=double), intent(in) :: q, alambda - real(kind=double), intent(out) :: alpha - real(kind=double) :: al, al2, alfinv, alfqc3, alfqc4, alfqc5, & - alfqs3, alfqs5 - real(kind=double) :: x1, x2, t, tt - real(kind=double) :: qsdt, qsdtt, qsct, qsctt, qs, q2 - real(kind=double) :: b0, b1, del, f, fp, as, as2 - real(kind=double) :: flav - ! The value of Lambda required corresponds to nflav=4 - qsdt=8.18 !! This is the value of 4m_c^2 - qsct=74.0 !! This is the value of 4m_b^2 - al2=alambda*alambda - q2=q*q - t = log (q2/al2) - ! CHECK: explicitly initialising ALFQC{3,4,5} (by AB) - alfqc3 = 0 - alfqc4 = 0 - alfqc5 = 0 - ith=0 - tt=t - qsdtt=qsdt/4. - qsctt=qsct/4. - al = alambda - al2=al*al - flav=4. - qs=al2*exp(T) - if (qs.lt.0.5d0) then !! running stops below 0.5 - qs=0.5d0 - t = log (qs/al2) - tt= t - end if - if (qs.lt.qsdtt .and. nflav.gt.3) GOTO 312 -11 continue - b0 = 11._double - 2./3._double * flav - x1 = 4.*PI/b0 - b1=102.-38.*FLAV/3. - X2=B1/B0**2 - AS2=X1/T*(1.-X2* log(T)/T) -5 AS=AS2 - F=-T+X1/AS-X2*log(X1/AS+X2) - FP=-X1/AS**2*(1.-X2/(X1/AS+X2)) - AS2=AS-F/FP - DEL=ABS(F/FP/AS) - IF((DEL-tol).GT.0.) goto 5 - ALPHA=AS2 -51 continue - IF(ITH.EQ.0) RETURN - GOTO (13,14,15) ITH - ! GOTO 5 -12 ITH=1 - T=log(QSCTT/AL2) - GOTO 11 -13 ALFQC4=ALPHA - FLAV=5. - ITH=2 - GOTO 11 -14 ALFQC5=ALPHA - ITH=3 - T=TT - GOTO 11 -15 ALFQS5=ALPHA - ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5 - ALPHA=1./ALFINV - RETURN -311 CONTINUE - B0=11-2.*FLAV/3. - X1=4.*PI/B0 - B1=102.-38.*FLAV/3. - X2=B1/B0**2 - AS2=X1/T*(1.-X2*log(T)/T) -35 AS=AS2 - F=-T+X1/AS-X2*log(X1/AS+X2) - FP=-X1/AS**2*(1.-X2/(X1/AS+X2)) - AS2=AS-F/FP - DEL=ABS(F/FP/AS) - IF((DEL-tol).GT.0.) goto 35 - ALPHA=AS2 -351 continue - IF(ITH.EQ.0) RETURN - GOTO (313,314,315) ITH -312 ITH=1 - T=log(QSDTT/AL2) - GOTO 311 -313 ALFQC4=ALPHA - FLAV=3. - ITH=2 - GOTO 311 -314 ALFQC3=ALPHA - ITH=3 - T=TT - GOTO 311 -315 ALFQS3=ALPHA - ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3 - ALPHA=1./ALFINV - end subroutine mrst_lambda + interface + module subroutine pdf_init (pdftype, prefix, verbose) + integer, intent(in) :: pdftype + type(string_t), intent(in), optional :: prefix + logical, intent(in), optional :: verbose + end subroutine pdf_init + module subroutine pdf_evolve (pdftype, x, q, f, fphoton) + integer, intent(in) :: pdftype + real(kind=default), intent(in) :: x, q + real(kind=default), intent(out), optional :: f(-6:6), fphoton + end subroutine pdf_evolve + module subroutine pdf_evolve_LHAPDF (set, x, q, ff) + integer, intent(in) :: set + double precision, intent(in) :: x, q + double precision, dimension(-6:6), intent(out) :: ff + end subroutine pdf_evolve_LHAPDF + module function pdf_get_name (pdftype) result (name) + integer, intent(in) :: pdftype + type(string_t) :: name + end function pdf_get_name + module function pdf_get_id (name) result (id) + type(string_t), intent(in) :: name + integer :: id + end function pdf_get_id + module function pdf_provides_photon (pdftype) result (flag) + integer, intent(in) :: pdftype + logical :: flag + end function pdf_provides_photon + module function pdf_getmass (nf) result (mass) + integer, intent(in) :: nf + real(kind=double) :: mass + end function pdf_getmass + module function pdf_alphas (pdftype, q) result (alphas) + real(kind=default), intent(in) :: q + real(kind=default) :: alphas + integer, intent(in) :: pdftype + end function pdf_alphas + module function pdf_alphas_LHAPDF (pdftype, q) result (alphas) + integer, intent(in) :: pdftype + real(kind=double), intent(in) :: q + real(kind=double) :: alphas + end function pdf_alphas_LHAPDF + end interface end module pdf_builtin Index: trunk/src/pdf_builtin/pdf_builtin_sub.f90 =================================================================== --- trunk/src/pdf_builtin/pdf_builtin_sub.f90 (revision 0) +++ trunk/src/pdf_builtin/pdf_builtin_sub.f90 (revision 8774) @@ -0,0 +1,1051 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! +! 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. +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! Wrap a common interface around the different PDF sets. + +submodule (pdf_builtin) pdf_builtin_s + + use string_utils, only: upper_case + use constants, only: PI + use diagnostics + use mrst2004qed + use cteq6pdf + use mstwpdf + use ct10pdf + use CJ_pdf + use ct14pdf + use ct18pdf + + implicit none + save + +! The available sets + integer, parameter :: nsets = 24 + integer, parameter :: & + CTEQ6M = 1, CTEQ6D = 2, CTEQ6L = 3, CTEQ6L1 = 4, & + MRST2004QEDp = 5, MRST2004QEDn = 6, MSTW2008LO = 7, MSTW2008NLO = 8, & + MSTW2008NNLO = 9, CT10 = 10, CJ12_max = 11, CJ12_mid = 12, & + CJ12_min = 13, MMHT2014LO = 14, MMHT2014NLO = 15, MMHT2014NNLO = 16, & + CT14LL = 17, CT14L = 18, CT14N = 19, CT14NN = 20, CJ15_LO = 21, & + CJ15_NLO = 22, CT18N = 23, CT18NN = 24 + +! Limits + real(kind=default), parameter :: & + cteq6_q_min = 1.3, cteq6_q_max = 10.E3, & + cteq6_x_min = 1.E-6, cteq6_x_max = 1., & + mrst2004qed_q_min = sqrt (1.26_default), mrst2004qed_q_max = sqrt (0.99E7_default), & + mrst2004qed_x_min = 1.01E-5, mrst2004qed_x_max = 1., & + mstw2008_q_min = 1.01, mstw2008_q_max = sqrt (0.99E9_default), & + mstw2008_x_min = 1.01E-6, mstw2008_x_max = 1., & + ct10_q_min = 1.3, ct10_q_max = 1.E5, & + ct10_x_min = 1.E-8, ct10_x_max = 1., & + cj12_q_min = 1.3, cj12_q_max = 1.E5, & + cj12_x_min = 1.E-6, cj12_x_max = 1., & + mmht2014_q_min = mstw2008_q_min, mmht2014_q_max = mstw2008_q_max, & + mmht2014_x_min = mstw2008_x_min, mmht2014_x_max = mstw2008_x_max, & + ct14_q_min = 1.3, ct14_q_max = 1.E5, ct14_x_min = 1.E-9, & + ct14_x_max = 1., ct18_q_min = 1.3, ct18_q_max = 1.E5, & + ct18_x_min = 1.E-9, ct18_x_max = 1. + +! Lambda_QCD and quark masses + + real(kind=double), dimension(6), parameter :: alam_ct10 = & + [ 0.37219423568859566_double, 0.37219423568859566_double, & + 0.37219423568859566_double, 0.37219423568859566_double, & + 0.32560730624033124_double, 0.22600000000000001_double ] + + real(kind=double), dimension(6), parameter :: alam_cteq6m = & + [ 0.37248648555333957_double, 0.37248648555333957_double, & + 0.37248648555333957_double, 0.37248648555333957_double, & + 0.32590307925177625_double, 0.22623045868581182_double ] + + real(kind=double), dimension(6), parameter :: alam_cteq6l = & + [ 0.37218603304249098_double, 0.37218603304249098_double, & + 0.37218603304249098_double, 0.37218603304249098_double, & + 0.32559900534407232_double, 0.22599353258372407_double ] + + real(kind=double), dimension(6), parameter :: alam_cteq6ll = & + [ 0.24560434095679567_double, 0.24560434095679567_double, & + 0.24560434095679567_double, 0.24560434095679567_double, & + 0.21495099184302788_double, 0.16499885346541945_double ] + + real(kind=double), dimension(0:5), parameter :: amhat = & + [ 0.0000000000000000_double, 0.0000000000000000_double, & + 0.0000000000000000_double, 0.0000000000000000_double, & + 1.3000000000000000_double, 4.500000000000000_double ] + + !!! Global parameter for minimal value of evolution + real(kind=double), parameter :: amn = 0.375_double, & + tol = .0005_double + + integer, parameter :: nhq = 2 + + real(kind=double), parameter :: & + ca = 3.D0, cf = 4./3.D0, tr = 0.5D0 + real(kind=double), parameter :: & + b00 = 11./3.d0 * ca, b01 = -4./3.d0 * tr, b10 = 34./3.d0 * ca**2, & + b11 = -20./3.d0 * ca*tr - 4.* cf*tr + + ! Init flags + integer :: cteq6_initialized = -1 + integer :: mstw2008_initialized = -1 + integer :: mmht2014_initialized = -1 + logical :: ct10_initialized = .false. + integer :: cj12_initialized = -1 + integer :: ct14_initialized = -1 + integer :: cj15_initialized = -1 + integer :: ct18_initialized = -1 + logical :: & + mrst2004qedp_initialized = .false., & + mrst2004qedn_initialized = .false. + type(string_t) :: mrst2004qedp_prefix, mrst2004qedn_prefix, & + mstw2008_prefix, mmht2014_prefix + +contains + +! Get PDF name + module function pdf_get_name (pdftype) result (name) + integer, intent(in) :: pdftype + type(string_t) :: name + select case (pdftype) + case (CTEQ6M) + name = var_str ("CTEQ6M") + case (CTEQ6D) + name = var_str ("CTEQ6D") + case (CTEQ6L) + name = var_str ("CTEQ6L") + case (CTEQ6L1) + name = var_str ("CTEQ6L1") + case (MRST2004QEDp) + name = var_str ("MRST2004QEDp") + case (MRST2004QEDn) + name = var_str ("MRST2004QEDn") + case (MSTW2008LO) + name = var_str ("MSTW2008LO") + case (MSTW2008NLO) + name = var_str ("MSTW2008NLO") + case (MSTW2008NNLO) + name = var_str ("MSTW2008NNLO") + case (CT10) + name = var_str ("CT10") + case (CJ12_max) + name = var_str ("CJ12_max") + case (CJ12_mid) + name = var_str ("CJ12_mid") + case (CJ12_min) + name = var_str ("CJ12_min") + case (CJ15_LO) + name = var_str ("CJ15LO") + case (CJ15_NLO) + name = var_str ("CJ15NLO") + case (MMHT2014LO) + name = var_str ("MMHT2014LO") + case (MMHT2014NLO) + name = var_str ("MMHT2014NLO") + case (MMHT2014NNLO) + name = var_str ("MMHT2014NNLO") + case (CT14LL) + name = var_str ("CT14LL") + case (CT14L) + name = var_str ("CT14L") + case (CT14N) + name = var_str ("CT14N") + case (CT14NN) + name = var_str ("CT14NN") + case (CT18N) + name = var_str ("CT18N") + case (CT18NN) + name = var_str ("CT18NN") + case default + call msg_fatal ("pdf_builtin: internal: invalid PDF set!") + end select + end function pdf_get_name + +! Get the ID of a PDF set + module function pdf_get_id (name) result (id) + type(string_t), intent(in) :: name + integer :: id + do id = 1, nsets + if (upper_case (pdf_get_name (id)) == upper_case (name)) return + end do + id = -1 + end function pdf_get_id + +! Query whether a PDF supplies a photon distribution + module function pdf_provides_photon (pdftype) result (flag) + integer, intent(in) :: pdftype + logical :: flag + select case (pdftype) + case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) + flag = .false. + case (MRST2004QEDp, MRST2004QEDn) + flag = .true. + case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO) + flag = .false. + case (CT10) + flag = .false. + case (CJ12_max, CJ12_mid, CJ12_min) + flag = .false. + case (CJ15_LO, CJ15_NLO) + flag = .false. + case (MMHT2014LO, MMHT2014NLO, MMHT2014NNLO) + flag = .false. + case (CT14LL, CT14L, CT14N, CT14NN) + flag = .false. + case (CT18N, CT18NN) + flag = .false. + case default + call msg_fatal ("pdf_builtin: internal: invalid PDF set!") + end select + end function pdf_provides_photon + +! Initialize a PDF + module subroutine pdf_init (pdftype, prefix, verbose) + integer, intent(in) :: pdftype + type(string_t), intent(in), optional :: prefix + type(string_t) :: mprefix + logical, intent(in), optional :: verbose + logical :: mverbose + if (present (prefix)) then + mprefix = prefix + else + mprefix = "" + end if + if (present (verbose)) then + mverbose = verbose + else + mverbose = .true. + end if + select case (pdftype) + case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) + if (cteq6_initialized == pdftype) return + call setctq6 (pdftype, char (mprefix)) + cteq6_initialized = pdftype + case (MRST2004QEDp) + if (mrst2004qedp_initialized) return + mrst2004qedp_initialized = .true. + mrst2004qedp_prefix = mprefix + case (MRST2004QEDn) + if (mrst2004qedn_initialized) return + mrst2004qedn_initialized = .true. + mrst2004qedn_prefix = mprefix + case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO) + if (mstw2008_initialized == pdftype) return + mstw2008_initialized = pdftype + mstw2008_prefix = mprefix + case (CT10) + if (ct10_initialized) return + call setct10 (char (mprefix), 100) + ct10_initialized = .true. + case (CJ12_max) + if (cj12_initialized == pdftype) return + call setCJ (char (mprefix), 300) + cj12_initialized = pdftype + case (CJ12_mid) + if (cj12_initialized == pdftype) return + call setCJ (char (mprefix), 200) + cj12_initialized = pdftype + case (CJ12_min) + if (cj12_initialized == pdftype) return + call setCJ (char (mprefix), 100) + cj12_initialized = pdftype + case (CJ15_LO) + if (cj15_initialized == pdftype) return + call setCJ (char (mprefix), 400) + cj15_initialized = pdftype + case (CJ15_NLO) + if (cj15_initialized == pdftype) return + call setCJ (char (mprefix), 500) + cj15_initialized = pdftype + case (MMHT2014LO, MMHT2014NLO, MMHT2014NNLO) + if (mmht2014_initialized == pdftype) return + mmht2014_initialized = pdftype + mmht2014_prefix = mprefix + case (CT14LL) + if (ct14_initialized == pdftype) return + ct14_initialized = pdftype + call setct14 (char (mprefix), 1) + case (CT14L) + if (ct14_initialized == pdftype) return + ct14_initialized = pdftype + call setct14 (char (mprefix), 2) + case (CT14N) + if (ct14_initialized == pdftype) return + ct14_initialized = pdftype + call setct14 (char (mprefix), 3) + case (CT14NN) + if (ct14_initialized == pdftype) return + ct14_initialized = pdftype + call setct14 (char (mprefix), 4) + case (CT18N) + if (ct18_initialized == pdftype) return + ct18_initialized = pdftype + call setct18 (char (prefix), 1) + case (CT18NN) + if (ct18_initialized == pdftype) return + ct18_initialized = pdftype + call setct18 (char (prefix), 2) + case default + call msg_fatal ("pdf_builtin: internal: invalid PDF set!") + end select + if (mverbose) call msg_message ("Initialized builtin PDF " // & + char (pdf_get_name (pdftype))) + !!! Up to now only proton + end subroutine pdf_init + +! Evolve PDF + module subroutine pdf_evolve (pdftype, x, q, f, fphoton) + integer, intent(in) :: pdftype + real(kind=default), intent(in) :: x, q + real(kind=default), intent(out), optional :: f(-6:6), fphoton + real(kind=double) :: mx, mq + real(kind=double) :: upv, dnv, ups, dns, str, chm, bot, glu, phot, & + sbar, bbar, cbar + type(string_t) :: setname + select case (pdftype) + case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) + if (cteq6_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (cteq6_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (cteq6_initialized)) // & + " cannot be used simultaneously") + mx = max (min (x, cteq6_x_max), cteq6_x_min) + mq = max (min (q, cteq6_q_max), cteq6_q_min) + if (present (f)) f = [ 0._double, & + ctq6pdf (-5, mx, mq), ctq6pdf (-4, mx, mq), & + ctq6pdf (-3, mx, mq), ctq6pdf (-1, mx, mq), & + ctq6pdf (-2, mx, mq), ctq6pdf ( 0, mx, mq), & + ctq6pdf ( 2, mx, mq), ctq6pdf ( 1, mx, mq), & + ctq6pdf ( 3, mx, mq), ctq6pdf ( 4, mx, mq), & + ctq6pdf ( 5, mx, mq), 0._double ] + if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & + char (pdf_get_name (pdftype)) // " which does not provide it!") + case (MRST2004QEDp) + if (.not. mrst2004qedp_initialized) call msg_fatal ( & + "pdf_builtin: internal: PDF set MRST2004QEDp requested without " // & + "initialization!") + mx = max (min (x, mrst2004qed_x_max), mrst2004qed_x_min) + mq = max (min (q, mrst2004qed_q_max), mrst2004qed_q_min) + call mrstqed (mx, mq, 1, upv, dnv, ups, dns, str, chm, bot, glu, phot, & + char (mrst2004qedp_prefix)) + if (present (f)) f = & + [ 0._double, bot, chm, str, ups, dns, glu, dns + dnv, ups + upv, & + str, chm, bot, 0._double ] / mx + if (present (fphoton)) fphoton = phot / mx + case (MRST2004QEDn) + if (.not. mrst2004qedn_initialized) call msg_fatal ( & + "pdf_builtin: internal: PDF set MRST2004QEDn requested without " // & + "initialization!") + mx = max (min (x, mrst2004qed_x_max), mrst2004qed_x_min) + mq = max (min (q, mrst2004qed_q_max), mrst2004qed_q_min) + call mrstqed (mx, mq, 2, upv, dnv, ups, dns, str, chm, bot, glu, phot, & + char (mrst2004qedn_prefix)) + if (present (f)) f = & + [ 0._double, bot, chm, str, ups, dns, glu, dns + dnv, ups + upv, & + str, chm, bot, 0._double ] / mx + if (present (fphoton)) fphoton = phot / mx + case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO) + if (mstw2008_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (mstw2008_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (mstw2008_initialized)) // & + " cannot be used simultaneously") + select case (pdftype) + case (MSTW2008LO) + setname = var_str ("mstw2008lo") + case (MSTW2008NLO) + setname = var_str ("mstw2008nlo") + case (MSTW2008NNLO) + setname = var_str ("mstw2008nnlo") + end select + mx = max (min (x, mstw2008_x_max), mstw2008_x_min) + mq = max (min (q, mstw2008_q_max), mstw2008_q_min) + call getmstwpdf (char (mstw2008_prefix), char (setname), 0, mx, mq, & + upv, dnv, ups, dns, str, sbar, chm, cbar, bot, bbar, glu, phot) + if (present (f)) f = & + [ 0._double, bbar, cbar, sbar, ups, dns, glu, dns + dnv, ups + upv, & + str, chm, bot, 0._double ] / mx + if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & + char (pdf_get_name (pdftype)) // " which does not provide it!") + case (CT10) + if (.not. ct10_initialized) call msg_fatal ( & + "pdf_builtin: internal: PDF set CT10 requested without " // & + "initialization!") + mx = max (min (x, ct10_x_max), ct10_x_min) + mq = max (min (q, ct10_q_max), ct10_q_min) + if (present (f)) f = [ 0._double, & + getct10pdf (-5, mx, mq), getct10pdf (-4, mx, mq), & + getct10pdf (-3, mx, mq), getct10pdf (-1, mx, mq), & + getct10pdf (-2, mx, mq), getct10pdf ( 0, mx, mq), & + getct10pdf ( 2, mx, mq), getct10pdf ( 1, mx, mq), & + getct10pdf ( 3, mx, mq), getct10pdf ( 4, mx, mq), & + getct10pdf ( 5, mx, mq), 0._double ] + if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & + "CT10 which does not provide it!") + case (CJ12_max, CJ12_mid, CJ12_min) + if (cj12_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (cj12_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (cj12_initialized)) // & + " cannot be used simultaneously") + mx = max (min (x, cj12_x_max), cj12_x_min) + mq = max (min (q, cj12_q_max), cj12_q_min) + if (present (f)) f = [ 0._double, & + CJpdf (-5, mx, mq), CJpdf (-4, mx, mq), CJpdf (-3, mx, mq), & + CJpdf (-2, mx, mq), CJpdf (-1, mx, mq), CJpdf ( 0, mx, mq), & + CJpdf ( 1, mx, mq), CJpdf ( 2, mx, mq), CJpdf ( 3, mx, mq), & + CJpdf ( 4, mx, mq), CJpdf ( 5, mx, mq), 0._double ] + if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & + "CJ12 which does not provide it!") + case (CJ15_LO, CJ15_NLO) + if (cj15_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (cj15_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (cj15_initialized)) // & + " cannot be used simultaneously") + mx = max (min (x, cj12_x_max), cj12_x_min) + mq = max (min (q, cj12_q_max), cj12_q_min) + if (present (f)) f = [ 0._double, & + CJpdf (-5, mx, mq), CJpdf (-4, mx, mq), CJpdf (-3, mx, mq), & + CJpdf (-2, mx, mq), CJpdf (-1, mx, mq), CJpdf ( 0, mx, mq), & + CJpdf ( 1, mx, mq), CJpdf ( 2, mx, mq), CJpdf ( 3, mx, mq), & + CJpdf ( 4, mx, mq), CJpdf ( 5, mx, mq), 0._double ] + if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & + "CJ15 which does not provide it!") + case (MMHT2014LO, MMHT2014NLO, MMHT2014NNLO) + if (mmht2014_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (mmht2014_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (mmht2014_initialized)) // & + " cannot be used simultaneously") + select case (pdftype) + case (MMHT2014LO) + setname = var_str ("mmht2014lo") + case (MMHT2014NLO) + setname = var_str ("mmht2014nlo") + case (MMHT2014NNLO) + setname = var_str ("mmht2014nnlo") + end select + mx = max (min (x, mmht2014_x_max), mmht2014_x_min) + mq = max (min (q, mmht2014_q_max), mmht2014_q_min) + call getmstwpdf (char (mmht2014_prefix), char (setname), 0, mx, mq, & + upv, dnv, ups, dns, str, sbar, chm, cbar, bot, bbar, glu, phot) + if (present (f)) f = & + [ 0._double, bbar, cbar, sbar, ups, dns, glu, dns + dnv, ups + upv, & + str, chm, bot, 0._double ] / mx + if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & + char (pdf_get_name (pdftype)) // " which does not provide it!") + case (CT14LL, CT14L, CT14N, CT14NN) + if (ct14_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (ct14_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (ct14_initialized)) // & + " cannot be used simultaneously") + mx = max (min (x, ct14_x_max), ct14_x_min) + mq = max (min (q, ct14_q_max), ct14_q_min) + if (present (f)) f = [ 0._double, & + getct14pdf (-5, mx, mq), getct14pdf (-4, mx, mq), & + getct14pdf (-3, mx, mq), getct14pdf (-1, mx, mq), & + getct14pdf (-2, mx, mq), getct14pdf ( 0, mx, mq), & + getct14pdf ( 2, mx, mq), getct14pdf ( 1, mx, mq), & + getct14pdf ( 3, mx, mq), getct14pdf ( 4, mx, mq), & + getct14pdf ( 5, mx, mq), 0._double ] + case (CT18N, CT18NN) + if (ct18_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (ct18_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (ct18_initialized)) // & + " cannot be used simultaneously") + mx = max (min (x, ct18_x_max), ct18_x_min) + mq = max (min (q, ct18_q_max), ct18_q_min) + if (present (f)) f = [ 0._double, & + getct18pdf (-5, mx, mq), getct18pdf (-4, mx, mq), & + getct18pdf (-3, mx, mq), getct18pdf (-1, mx, mq), & + getct18pdf (-2, mx, mq), getct18pdf ( 0, mx, mq), & + getct18pdf ( 2, mx, mq), getct18pdf ( 1, mx, mq), & + getct18pdf ( 3, mx, mq), getct18pdf ( 4, mx, mq), & + getct18pdf ( 5, mx, mq), 0._double ] + if (present (fphoton)) call msg_fatal ("photon pdf requested for " // & + char (pdf_get_name (pdftype)) // " which does not provide it!") + case default + call msg_fatal ("pdf_builtin: internal: invalid PDF set!") + end select + end subroutine pdf_evolve + +! included for compatibility with LHAPDF +! use a double precision array for the pdfs instead of a +! default precision + module subroutine pdf_evolve_LHAPDF (set, x, q, ff) + integer, intent(in) :: set + double precision, intent(in) :: x, q + real(kind=default) :: dx, dq + double precision, dimension(-6:6), intent(out) :: ff + + real(kind=default) :: f(-6:6) + dx = x + dq = q + + call pdf_evolve (set, dx, dq, f) + ff = f + end subroutine pdf_evolve_LHAPDF + +! PDF-specific running alphas + module function pdf_alphas (pdftype, q) result (alphas) + real(kind=double) :: as + real(kind=default), intent(in) :: q + real(kind=double) :: qdummy + real(kind=default) :: alphas + integer, intent(in) :: pdftype + integer, parameter :: nset_dummy = 1 + qdummy = q + select case (pdftype) + case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) + if (cteq6_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (cteq6_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (cteq6_initialized)) // & + " cannot be used simultaneously") + select case (pdftype) + case (CTEQ6M, CTEQ6D) + as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6m) + case (CTEQ6L) + as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6l) + case (CTEQ6L1) + as = PI*pdf_builtin_alpi (qdummy,1,5,alam_cteq6ll) + case default + call msg_fatal ("Unrecognized PDF setup in evolution of alphas.") + end select + case (CT10) + if (.not. ct10_initialized) call msg_fatal ( & + "pdf_builtin: internal: PDF set CT10 requested without " // & + "initialization!") + as = PI*pdf_builtin_alpi (qdummy,2,5,alam_ct10) + case (MRST2004QEDp) + if (.not. mrst2004qedp_initialized) call msg_fatal ( & + "pdf_builtin: internal: PDF set MRST2004QEDp requested without " // & + "initialization!") + call mrst_alphas (as,qdummy) + case (MRST2004QEDn) + if (.not. mrst2004qedn_initialized) call msg_fatal ( & + "pdf_builtin: internal: PDF set MRST2004QEDn requested without " // & + "initialization!") + call mrst_alphas (as,qdummy) + case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO) + if (mstw2008_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (mstw2008_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (mstw2008_initialized)) // & + " cannot be used simultaneously") + select case (pdftype) + case (MSTW2008lo) + call pdf_builtin_alfa (as,qdummy,0) + case (MSTW2008nlo) + call pdf_builtin_alfa (as,qdummy,1) + case (MSTW2008nnlo) + call pdf_builtin_alfa (as,qdummy,2) + end select + case (CJ12_max, CJ12_mid, CJ12_min) + if (cj12_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (cj12_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (cj12_initialized)) // & + " cannot be used simultaneously") + as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6m) + case (CJ15_LO, CJ15_NLO) + if (cj15_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (cj15_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (cj15_initialized)) // & + " cannot be used simultaneously") + select case (pdftype) + case (CJ15_NLO) + as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6m) + case (CJ15_LO) + as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6l) + end select + case (MMHT2014LO, MMHT2014NLO, MMHT2014NNLO) + if (mmht2014_initialized < 0) & + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (mmht2014_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (mmht2014_initialized)) // & + " cannot be used simultaneously") + select case (pdftype) + case (MMHT2014lo) + call pdf_builtin_alfa (as,qdummy,0) + case (MMHT2014nlo) + call pdf_builtin_alfa (as,qdummy,1) + case (MMHT2014nnlo) + call pdf_builtin_alfa (as,qdummy,2) + end select + case (CT14LL, CT14L, CT14N, CT14NN) + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (ct14_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (ct14_initialized)) // & + " cannot be used simultaneously") + as = CT14Alphas(qdummy) + case (CT18N, CT18NN) + call msg_fatal ("pdf_builtin: internal: PDF set " // & + char (pdf_get_name (pdftype)) // " requested without initialization!") + if (ct18_initialized /= pdftype) & + call msg_fatal ( & + "PDF sets " // char (pdf_get_name (pdftype)) // " and " // & + char (pdf_get_name (ct18_initialized)) // & + " cannot be used simultaneously") + as = CT18Alphas(qdummy) + case default + call msg_fatal ("pdf_builtin: internal: invalid PDF set!") + end select + alphas = as + end function pdf_alphas + + module function pdf_alphas_LHAPDF (pdftype, q) result (alphas) + integer, intent(in) :: pdftype + real(kind=double), intent(in) :: q + real(kind=double) :: alphas + real(kind=default) :: q_def + q_def = q + alphas = pdf_alphas (pdftype, q_def) + end function pdf_alphas_LHAPDF + + module function pdf_getmass (nf) result (mass) + integer, intent(in) :: nf + real(kind=double) :: mass + select case (abs (nf)) + case (1:3) + mass = 0._double + case (4) + mass = 1.3_double + case (5) + mass = 4.5_double + case default + call msg_fatal ("PDF builtin: invalid PDG code for quark mass.") + end select + end function pdf_getmass + + function pdf_builtin_alpi (AMU,NORDER,NF,alam) + integer, intent(in) :: nf, norder + real(kind=double), intent(in) :: amu + real(kind=double), dimension(6), intent(in) :: alam + real(kind=double) :: alm + real(kind=double) :: pdf_builtin_alpi + integer :: irt, n_eff + n_eff = num_flavor (AMU,NF) + alm = alam (n_eff+1) + pdf_builtin_alpi = pdf_builtin_alpqcd (norder, n_eff, amu/alm, irt) + select case (irt) + case (1) + call msg_warning ("AMU < ALAM in pdf_builtin_alpi") + case (2) + call msg_warning ("pdf_builtin_alpi > 3; Be aware!") + end select + end function pdf_builtin_alpi + + double precision function pdf_builtin_ALPQCD (IRDR, NF, RML, IRT) + real(kind=double) :: b0, b1 + integer, intent(in) :: irdr, nf + integer :: irt + real(kind=double), intent(in) :: rml + real(kind=double) :: al, aln + real(kind=double) :: rm2 + IRT = 0 + if (IRDR .LT. 1 .OR. IRDR .GT. 2) THEN + print *, & + & 'Order out of range in pdf_builtin_ALPQCD: IRDR = ', IRDR + STOP + end if + b0 = (11.d0*CA - 2.* NF) / 3.d0 + b1 = (34.d0*CA**2 - 10.d0*CA*NF - 6.d0*CF*NF) / 3.d0 + rm2 = rml**2 + if (RM2 .LE. 1.) THEN + IRT = 1 + pdf_builtin_ALPQCD = 99. + RETURN + end if + aln = log (RM2) + al = 4.d0/ B0 / aln + IF (IRDR .GE. 2) AL = AL * (1.d0-B1*log(ALN) / ALN / B0**2) + if (AL .GE. 3.) THEN + IRT = 2 + end if + pdf_builtin_ALPQCD = AL + end function pdf_builtin_ALPQCD + + function num_flavor (amu,nf) + real(kind=double), intent(in) :: amu + integer, intent(in) :: nf + integer :: num_flavor, i + num_flavor = nf - nhq + IF ((num_flavor .EQ. nf) .OR. (amu .LE. amn)) GOTO 20 + do 10 I = nf - NHQ + 1, nf + if (amu .GE. amhat(I)) THEN + num_flavor = I + else + goto 20 + end if +10 continue +20 return + end function num_flavor + + subroutine pdf_builtin_alfa (alfas,Qalfa,norder) + integer, intent(in) :: norder + real(kind=double), intent(in) :: Qalfa + real(kind=double), intent(out) :: alfas + real(kind=double) :: alphasq0 + real(kind=double) :: mCharmSave, mBottomSave, mTopSave + mTopSave = 10000000000.0_double + mBottomSave = 4.75_double + mCharmSave = 1.4_double + select case (norder) + case (0) + alphasq0 = 0.68183_double + case (1) + alphasq0 = 0.491278_double + case (2) + alphasq0 = 0.45077_double + case default + call msg_fatal ("Wrong MSTW order!") + end select + if (Qalfa.GT.mTopSave) call msg_fatal & + ("No MSTW alphas for such high scales.") + alfas = pdf_builtin_mstw_alphas (Qalfa,norder,mTopSave**2) + end subroutine pdf_builtin_alfa + + function pdf_builtin_mstw_alphas (mur, norder, m2t) result (alfas) + real(kind=double), intent(in) :: mur + real(kind=double) :: alfas + integer, intent(in) :: norder + integer :: nf + real(kind=double), intent(in) :: m2t + real(kind=double) :: as0, asc, asb + real(kind=double) :: r2, m2, asi, asf, r20, r2b, r2c, r2t + real(kind=double), parameter :: m20 = 1.0_double, m2c = 1.96_double, & + m2b = 22.5625_double + real(kind=double) :: logfr + integer :: nff + real(kind=double) :: ast + nff = 4 + select case (norder) + case (0) + as0 = 5.4258307424173556E-002_double + asc = 4.0838232991033577E-002_double + asb = 2.2297506503539639E-002_double + case (1) + as0 = 3.9094820221093209E-002_double + asc = 3.0202751103489092E-002_double + asb = 1.7751676069301441E-002_double + case (2) + as0 = 3.5871136848766867E-002_double + asc = 2.8092349518054758E-002_double + asb = 1.6936586710596203E-002_double + case default + call msg_fatal ("Wrong evolution order in MSTW alphas evolution.") + end select + r2 = mur**2 + m2 = r2 * exp(+logfr) + if (m2 .gt. m2t) then + nf = 6 + r2t = m2t * r2/m2 + asi = ast + asf = mstw_alphas (r2,r2t,ast,nf,norder) + else + if (m2 .gt. m2b) then + nf = 5 + r2b = m2b * r2/m2 + asi = asb + asf = mstw_alphas (r2,r2b,asb,nf,norder) + else + if (m2 .gt. m2c) then + nf = 4 + r2c = m2c * r2/m2 + asi = asc + asf = mstw_alphas (r2,r2c,asc,nf,norder) + else + nf = 3 + r20 = m20 * r2/m2 + asi = as0 + asf = mstw_alphas (r2,r20,as0,nf,norder) + end if + end if + end if + alfas = 4.D0*PI*ASF + end function pdf_builtin_mstw_alphas + + function mstw_alphas (r2, r20, as0, nf, naord) result (as) + real(kind=double) :: as + real(kind=double), intent(in) :: r2, r20, as0 + integer, intent(in) :: nf, naord + integer, parameter :: nastps = 20 + integer :: k1 + real(kind=double), parameter :: sxth = 1./6. + real(kind=double) :: lrrat, dlr, xk0, xk1, xk2, xk3 + as = as0 + lrrat = log (r2/r20) + dlr = lrrat / nastps + ! ..Solution of the evolution equation depending on NAORD + ! (fourth-order Runge-Kutta beyond the leading order) + select case (naord) + case (0) + as = as0 / (1.+ beta0(nf) * as0 * lrrat) + case (1) + do k1 = 1, nastps + xk0 = dlr * fbeta1 (as, nf) + xk1 = dlr * fbeta1 (as + 0.5 * xk0, nf) + xk2 = dlr * fbeta1 (as + 0.5 * xk1, nf) + xk3 = dlr * fbeta1 (as + xk2, nf) + as = as + sxth * (xk0 + 2.* xk1 + 2.* xk2 + xk3) + end do + case (2) + do k1 = 1, nastps + xk0 = dlr * fbeta2 (as, nf) + xk1 = dlr * fbeta2 (as + 0.5 * xk0, nf) + xk2 = dlr * fbeta2 (as + 0.5 * xk1, nf) + xk3 = dlr * fbeta2 (as + xk2, nf) + as = as + sxth * (xk0 + 2.* xk1 + 2.* xk2 + xk3) + end do + case (3) + do k1 = 1, nastps + xk0 = dlr * fbeta3 (as, nf) + xk1 = dlr * fbeta3 (as + 0.5 * xk0, nf) + xk2 = dlr * fbeta3 (as + 0.5 * xk1, nf) + xk3 = dlr * fbeta3 (as + xk2, nf) + as = as + sxth * (xk0 + 2.* xk1 + 2.* xk2 + xk3) + end do + end select + end function mstw_alphas + + double precision function beta0 (nf) + integer, intent(in) :: nf + beta0 = b00 + b01 * nf + end function beta0 + + double precision function beta1 (nf) + integer, intent(in) :: nf + beta1 = b10 + b11 * nf + end function beta1 + + double precision function beta2 (nf) + integer, intent(in) :: nf + beta2 = 1428.50 - 279.611 * nf + 6.01852 * nf**2 + end function beta2 + + double precision function beta3 (nf) + integer, intent(in) :: nf + beta3 = 29243.0 - 6946.30 * nf + 405.089 * nf**2 + & + 1.49931 * nf**3 + end function beta3 + + ! ..The beta functions FBETAn at N^nLO for n = 1, 2, and 3 + double precision function fbeta1 (a, nf) + real(kind=double), intent(in) :: a + integer, intent(in) :: nf + fbeta1 = - a**2 * ( beta0(nf) + a * beta1(nf) ) + end function fbeta1 + + double precision function fbeta2 (a, nf) + real(kind=double), intent(in) :: a + integer, intent(in) :: nf + fbeta2 = - a**2 * ( beta0(nf) + a * ( beta1(nf) & + + a * beta2(nf) ) ) + end function fbeta2 + + double precision function fbeta3 (a, nf) + real(kind=double), intent(in) :: a + integer, intent(in) :: nf + fbeta3 = - a**2 * ( beta0(nf) + a * ( beta1(nf) & + + a * ( beta2(nf) + a * beta3(nf)) ) ) + end function fbeta3 + + subroutine mrst_alphas (alpha, q) + real(kind=double), intent(in) :: q + real(kind=double), intent(out) :: alpha + integer :: idir + integer, parameter :: nflav = 4 + real(kind=double) :: alphas, qz2, qz, alambda, astep, & + tol2 + alphas = 0.1205_double + qz2=8315._double + qz = sqrt(qz2) + alambda = 0.3000_double + astep = 0.010_double + tol2 = 0.0000001_double + idir = +1 +100 continue + alambda = alambda + idir*astep + call mrst_lambda (nflav,alpha,qz,alambda) + if (idir*(alphas-alpha).gt.0.0) then + goto 200 + else + astep = 0.5*astep + idir = -1*idir + goto 200 + end if +200 continue + if(abs(alpha-alphas).gt.tol2) goto 100 + ! alambda found -- save it !!!! + ! next call mrst_lambda to get alphas at q with the correct alambda + call mrst_lambda (nflav,alpha,q,alambda) + RETURN + end subroutine mrst_alphas + + subroutine mrst_lambda (nflav,alpha,Q,alambda) + integer, intent(in) :: nflav + integer :: ith + real(kind=double), intent(in) :: q, alambda + real(kind=double), intent(out) :: alpha + real(kind=double) :: al, al2, alfinv, alfqc3, alfqc4, alfqc5, & + alfqs3, alfqs5 + real(kind=double) :: x1, x2, t, tt + real(kind=double) :: qsdt, qsdtt, qsct, qsctt, qs, q2 + real(kind=double) :: b0, b1, del, f, fp, as, as2 + real(kind=double) :: flav + ! The value of Lambda required corresponds to nflav=4 + qsdt=8.18 !! This is the value of 4m_c^2 + qsct=74.0 !! This is the value of 4m_b^2 + al2=alambda*alambda + q2=q*q + t = log (q2/al2) + ! CHECK: explicitly initialising ALFQC{3,4,5} (by AB) + alfqc3 = 0 + alfqc4 = 0 + alfqc5 = 0 + ith=0 + tt=t + qsdtt=qsdt/4. + qsctt=qsct/4. + al = alambda + al2=al*al + flav=4. + qs=al2*exp(T) + if (qs.lt.0.5d0) then !! running stops below 0.5 + qs=0.5d0 + t = log (qs/al2) + tt= t + end if + if (qs.lt.qsdtt .and. nflav.gt.3) GOTO 312 +11 continue + b0 = 11._double - 2./3._double * flav + x1 = 4.*PI/b0 + b1=102.-38.*FLAV/3. + X2=B1/B0**2 + AS2=X1/T*(1.-X2* log(T)/T) +5 AS=AS2 + F=-T+X1/AS-X2*log(X1/AS+X2) + FP=-X1/AS**2*(1.-X2/(X1/AS+X2)) + AS2=AS-F/FP + DEL=ABS(F/FP/AS) + IF((DEL-tol).GT.0.) goto 5 + ALPHA=AS2 +51 continue + IF(ITH.EQ.0) RETURN + GOTO (13,14,15) ITH + ! GOTO 5 +12 ITH=1 + T=log(QSCTT/AL2) + GOTO 11 +13 ALFQC4=ALPHA + FLAV=5. + ITH=2 + GOTO 11 +14 ALFQC5=ALPHA + ITH=3 + T=TT + GOTO 11 +15 ALFQS5=ALPHA + ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5 + ALPHA=1./ALFINV + RETURN +311 CONTINUE + B0=11-2.*FLAV/3. + X1=4.*PI/B0 + B1=102.-38.*FLAV/3. + X2=B1/B0**2 + AS2=X1/T*(1.-X2*log(T)/T) +35 AS=AS2 + F=-T+X1/AS-X2*log(X1/AS+X2) + FP=-X1/AS**2*(1.-X2/(X1/AS+X2)) + AS2=AS-F/FP + DEL=ABS(F/FP/AS) + IF((DEL-tol).GT.0.) goto 35 + ALPHA=AS2 +351 continue + IF(ITH.EQ.0) RETURN + GOTO (313,314,315) ITH +312 ITH=1 + T=log(QSDTT/AL2) + GOTO 311 +313 ALFQC4=ALPHA + FLAV=3. + ITH=2 + GOTO 311 +314 ALFQC3=ALPHA + ITH=3 + T=TT + GOTO 311 +315 ALFQS3=ALPHA + ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3 + ALPHA=1./ALFINV + end subroutine mrst_lambda + +end submodule pdf_builtin_s Index: trunk/src/utilities/Makefile.am =================================================================== --- trunk/src/utilities/Makefile.am (revision 8773) +++ trunk/src/utilities/Makefile.am (revision 8774) @@ -1,226 +1,226 @@ ## Makefile.am -- Makefile for 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. # ######################################################################## ## The files in this directory are simple utilities used by WHIZARD ## We create a library which is still to be combined with auxiliary libs. noinst_LTLIBRARIES = libutilities.la check_LTLIBRARIES = libutilities_ut.la libutilities_la_SOURCES = \ $(UTILITIES_MODULES) \ $(UTILITIES_SUBMODULES) UTILITIES_MODULES = \ file_utils.f90 \ file_registries.f90 \ string_utils.f90 \ format_utils.f90 \ format_defs.f90 \ numeric_utils.f90 \ binary_tree.f90 \ array_list.f90 \ queue.f90 \ iterator.f90 UTILITIES_SUBMODULES = \ file_utils_sub.f90 \ file_registries_sub.f90 \ string_utils_sub.f90 \ format_utils_sub.f90 \ numeric_utils_sub.f90 \ binary_tree_sub.f90 \ array_list_sub.f90 \ queue_sub.f90 \ iterator_sub.f90 libutilities_ut_la_SOURCES = \ binary_tree_ut.f90 binary_tree_uti.f90 \ array_list_ut.f90 array_list_uti.f90 \ iterator_ut.f90 iterator_uti.f90 ## Omitting this would exclude it from the distribution dist_noinst_DATA = utilities.nw # Modules and installation # Dump module names into file Modules execmoddir = $(fmoddir)/whizard nodist_execmod_HEADERS = \ ${UTILITIES_MODULES:.f90=.$(FCMOD)} -#Submodules must not be included here +# Submodules must not be included here # Dump module names into file Modules libutilities_Modules = ${UTILITIES_MODULES:.f90=} ${libutilities_ut_la_SOURCES:.f90=} Modules: Makefile @for module in $(libutilities_Modules); do \ echo $$module >> $@.new; \ done @if diff $@ $@.new -q >/dev/null; then \ rm $@.new; \ else \ mv $@.new $@; echo "Modules updated"; \ fi BUILT_SOURCES = Modules ## Fortran module dependencies # Get module lists from other directories module_lists = \ ../basics/Modules \ ../testing/Modules $(module_lists): $(MAKE) -C `dirname $@` Modules Module_dependencies.sed: $(libutilities_la_SOURCES) \ $(libutilities_ut_la_SOURCES) Module_dependencies.sed: $(module_lists) @rm -f $@ echo 's/, *only:.*//' >> $@ echo 's/, *&//' >> $@ echo 's/, *.*=>.*//' >> $@ echo 's/$$/.lo/' >> $@ for list in $(module_lists); do \ dir="`dirname $$list`"; \ for mod in `cat $$list`; do \ echo 's!: '$$mod'.lo$$!': $$dir/$$mod'.lo!' >> $@; \ done \ done DISTCLEANFILES = Module_dependencies.sed # The following line just says # include Makefile.depend # but in a portable fashion (depending on automake's AM_MAKE_INCLUDE @am__include@ @am__quote@Makefile.depend@am__quote@ Makefile.depend: Module_dependencies.sed Makefile.depend: $(libutilities_la_SOURCES) $(libutilities_ut_la_SOURCES) @rm -f $@ for src in $^; do \ module="`basename $$src | sed 's/\.f[90][0358]//'`"; \ grep '^ *use ' $$src \ | grep -v '!NODEP!' \ | sed -e 's/^ *use */'$$module'.lo: /' \ -f Module_dependencies.sed; \ done > $@ DISTCLEANFILES += Makefile.depend # Fortran90 module files are generated at the same time as object files .lo.$(FCMOD): @: # touch $@ AM_FCFLAGS = -I../basics -I../testing ######################################################################## # For the moment, the submodule dependencies will be hard-coded file_utils_sub.lo: file_utils.lo file_registries_sub.lo: file_registries.lo string_utils_sub.lo: string_utils.lo format_utils_sub.lo: format_utils.lo numeric_utils_sub.lo: numeric_utils.lo binary_tree_sub.lo: binary_tree.lo array_list_sub.lo: array_list.lo queue_sub.lo: queue.lo iterator_sub.lo: iterator.lo ######################################################################## ## Default Fortran compiler options ## Profiling if FC_USE_PROFILING AM_FCFLAGS += $(FCFLAGS_PROFILING) endif ## OpenMP if FC_USE_OPENMP AM_FCFLAGS += $(FCFLAGS_OPENMP) endif ## MPI if FC_USE_MPI AM_FCFLAGS += $(FCFLAGS_MPI) endif ######################################################################## ## Non-standard targets and dependencies ## (Re)create F90 sources from NOWEB source. if NOWEB_AVAILABLE PRELUDE = $(top_srcdir)/src/noweb-frame/whizard-prelude.nw POSTLUDE = $(top_srcdir)/src/noweb-frame/whizard-postlude.nw utilities.stamp: $(PRELUDE) $(srcdir)/utilities.nw $(POSTLUDE) @rm -f utilities.tmp @touch utilities.tmp for src in $(libutilities_la_SOURCES) $(libutilities_ut_la_SOURCES); do \ $(NOTANGLE) -R[[$$src]] $^ | $(CPIF) $$src; \ done @mv -f utilities.tmp utilities.stamp $(libutilities_la_SOURCES) $(libutilities_ut_la_SOURCES): utilities.stamp ## Recover from the removal of $@ @if test -f $@; then :; else \ rm -f utilities.stamp; \ $(MAKE) $(AM_MAKEFLAGS) utilities.stamp; \ fi endif ######################################################################## ## Non-standard cleanup tasks ## Remove sources that can be recreated using NOWEB if NOWEB_AVAILABLE maintainer-clean-noweb: -rm -f *.f90 *.c endif .PHONY: maintainer-clean-noweb ## Remove those sources also if builddir and srcdir are different if NOWEB_AVAILABLE clean-noweb: test "$(srcdir)" != "." && rm -f *.f90 *.c || true endif .PHONY: clean-noweb ## Remove F90 module files clean-local: clean-noweb -rm -f utilities.stamp utilities.tmp -rm -f *.$(FCMOD) if FC_SUBMODULES -rm -f *.smod *.sub endif ## Remove backup files maintainer-clean-backup: -rm -f *~ .PHONY: maintainer-clean-backup ## Register additional clean targets maintainer-clean-local: maintainer-clean-noweb maintainer-clean-backup Index: trunk/src/testing/Makefile.am =================================================================== --- trunk/src/testing/Makefile.am (revision 8773) +++ trunk/src/testing/Makefile.am (revision 8774) @@ -1,183 +1,193 @@ ## Makefile.am -- Makefile for 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. # ######################################################################## ## The files in this directory handle unit tests in Fortran. ## We create a library which is still to be combined with auxiliary libs. noinst_LTLIBRARIES = libtesting.la libtesting_la_SOURCES = \ + $(TESTING_MODULES) \ + $(TESTING_SUBMODULES) + +TESTING_MODULES = \ unit_tests.f90 +TESTING_SUBMODULES = \ + unit_tests_sub.f90 ## Omitting this would exclude it from the distribution dist_noinst_DATA = testing.nw # Modules and installation # Dump module names into file Modules execmoddir = $(fmoddir)/whizard nodist_execmod_HEADERS = \ - ${libtesting_la_SOURCES:.f90=.$(FCMOD)} + ${TESTING_MODULES:.f90=.$(FCMOD)} +# Submodules must not be included here # Dump module names into file Modules -libtesting_Modules = ${libtesting_la_SOURCES:.f90=} +libtesting_Modules = ${TESTING_MODULES:.f90=} Modules: Makefile @for module in $(libtesting_Modules); do \ echo $$module >> $@.new; \ done @if diff $@ $@.new -q >/dev/null; then \ rm $@.new; \ else \ mv $@.new $@; echo "Modules updated"; \ fi BUILT_SOURCES = Modules ## Fortran module dependencies # Get module lists from other directories module_lists = \ ../basics/Modules $(module_lists): $(MAKE) -C `dirname $@` Modules Module_dependencies.sed: $(libtesting_la_SOURCES) Module_dependencies.sed: $(module_lists) @rm -f $@ echo 's/, *only:.*//' >> $@ echo 's/, *&//' >> $@ echo 's/, *.*=>.*//' >> $@ echo 's/$$/.lo/' >> $@ for list in $(module_lists); do \ dir="`dirname $$list`"; \ for mod in `cat $$list`; do \ echo 's!: '$$mod'.lo$$!': $$dir/$$mod'.lo!' >> $@; \ done \ done DISTCLEANFILES = Module_dependencies.sed # The following line just says # include Makefile.depend # but in a portable fashion (depending on automake's AM_MAKE_INCLUDE @am__include@ @am__quote@Makefile.depend@am__quote@ Makefile.depend: Module_dependencies.sed Makefile.depend: $(libtesting_la_SOURCES) @rm -f $@ for src in $^; do \ module="`basename $$src | sed 's/\.f[90][0358]//'`"; \ grep '^ *use ' $$src \ | grep -v '!NODEP!' \ | sed -e 's/^ *use */'$$module'.lo: /' \ -f Module_dependencies.sed; \ done > $@ DISTCLEANFILES += Makefile.depend # Fortran90 module files are generated at the same time as object files .lo.$(FCMOD): @: # touch $@ AM_FCFLAGS = -I../basics -I../utilities -I../system +######################################################################## +# For the moment, the submodule dependencies will be hard-coded +unit_tests_sub.lo: unit_tests.lo ######################################################################## ## Default Fortran compiler options ## Profiling if FC_USE_PROFILING AM_FCFLAGS += $(FCFLAGS_PROFILING) endif ## OpenMP if FC_USE_OPENMP AM_FCFLAGS += $(FCFLAGS_OPENMP) endif ## MPI if FC_USE_MPI AM_FCFLAGS += $(FCFLAGS_MPI) endif ######################################################################## ## Non-standard targets and dependencies ## (Re)create F90 sources from NOWEB source. if NOWEB_AVAILABLE PRELUDE = $(top_srcdir)/src/noweb-frame/whizard-prelude.nw POSTLUDE = $(top_srcdir)/src/noweb-frame/whizard-postlude.nw testing.stamp: $(PRELUDE) $(srcdir)/testing.nw $(POSTLUDE) @rm -f testing.tmp @touch testing.tmp for src in $(libtesting_la_SOURCES); do \ $(NOTANGLE) -R[[$$src]] $^ | $(CPIF) $$src; \ done @mv -f testing.tmp testing.stamp $(libtesting_la_SOURCES): testing.stamp ## Recover from the removal of $@ @if test -f $@; then :; else \ rm -f testing.stamp; \ $(MAKE) $(AM_MAKEFLAGS) testing.stamp; \ fi endif ######################################################################## ## Non-standard cleanup tasks ## Remove sources that can be recreated using NOWEB if NOWEB_AVAILABLE maintainer-clean-noweb: -rm -f *.f90 *.c endif .PHONY: maintainer-clean-noweb ## Remove those sources also if builddir and srcdir are different if NOWEB_AVAILABLE clean-noweb: test "$(srcdir)" != "." && rm -f *.f90 *.c || true endif .PHONY: clean-noweb ## Remove F90 module files clean-local: clean-noweb -rm -f testing.stamp testing.tmp -rm -f *.$(FCMOD) if FC_SUBMODULES - -rm -f *.smod + -rm -f *.smod *.sub endif ## Remove backup files maintainer-clean-backup: -rm -f *~ .PHONY: maintainer-clean-backup ## Register additional clean targets maintainer-clean-local: maintainer-clean-noweb maintainer-clean-backup Index: trunk/src/testing/testing.nw =================================================================== --- trunk/src/testing/testing.nw (revision 8773) +++ trunk/src/testing/testing.nw (revision 8774) @@ -1,342 +1,418 @@ % -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*- % WHIZARD code as NOWEB source \chapter{Testing} \includemodulegraph{testing} This part contains tools for automatic testing. \begin{description} \item[unit\_tests] A handler that executes test procedures and compares and collects results. \end{description} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Unit tests} We provide functionality for automated unit tests. Each test is required to produce output which is compared against a reference file. If the two are identical, we signal success. Otherwise, we signal failure and write the output to a file. <<[[unit_tests.f90]]>>= <> module unit_tests <> - use io_units <> <> -<> - <> <> + interface +<> + end interface + +end module unit_tests +@ %def unit_tests +@ +<<[[unit_tests_sub.f90]]>>= +<> + +submodule (unit_tests) unit_tests_s + + use io_units + + implicit none + +<> + contains <> -end module unit_tests -@ %def unit_tests +end submodule unit_tests_s + +@ %def unit_tests_s @ \subsection{Parameters} Building blocks of file names. The directory names and suffixes are hard-coded here, and they must reflect actual Makefile targets where applicable. <>= character(*), parameter :: ref_prefix = "ref-output/" character(*), parameter :: ref = ".ref" character(*), parameter :: err_prefix = "err-output/" character(*), parameter :: err = ".out" @ %def prefix ref err @ \subsection{Type for storing test results} We store the results of the individual unit tests in a linked list. Here is the entry: <>= public :: test_results_t <>= type :: test_result_t logical :: success = .false. type(string_t) :: name type(string_t) :: description type(test_result_t), pointer :: next => null () end type test_result_t type :: test_results_t private type(test_result_t), pointer :: first => null () type(test_result_t), pointer :: last => null () integer :: n_success = 0 integer :: n_failure = 0 contains <> end type test_results_t @ %def test_result_t test_results_t @ Add a test result. <>= procedure :: add => test_results_add +<>= + module subroutine test_results_add (list, name, description, success) + class(test_results_t), intent(inout) :: list + character(len=*), intent(in) :: name + character(len=*), intent(in) :: description + logical, intent(in) :: success + end subroutine test_results_add <>= - subroutine test_results_add (list, name, description, success) + module subroutine test_results_add (list, name, description, success) class(test_results_t), intent(inout) :: list character(len=*), intent(in) :: name character(len=*), intent(in) :: description logical, intent(in) :: success type(test_result_t), pointer :: result allocate (result) result%success = success result%name = name result%description = description if (associated (list%first)) then list%last%next => result else list%first => result end if list%last => result if (success) then list%n_success = list%n_success + 1 else list%n_failure = list%n_failure + 1 end if end subroutine test_results_add @ %def test_results_add @ Display the current state. <>= procedure, private :: write => test_results_write +<>= + module subroutine test_results_write (list, u) + class(test_results_t), intent(in) :: list + integer, intent(in) :: u + end subroutine test_results_write <>= - subroutine test_results_write (list, u) + module subroutine test_results_write (list, u) class(test_results_t), intent(in) :: list integer, intent(in) :: u type(test_result_t), pointer :: result write (u, "(A)") "*** Test Summary ***" if (list%n_success > 0) then write (u, "(2x,A)") "Success:" result => list%first do while (associated (result)) if (result%success) write (u, "(4x,A,': ',A)") & char (result%name), char (result%description) result => result%next end do end if if (list%n_failure > 0) then write (u, "(2x,A)") "Failure:" result => list%first do while (associated (result)) if (.not. result%success) write (u, "(4x,A,': ',A)") & char (result%name), char (result%description) result => result%next end do end if write (u, "(A,I0)") "Total = ", list%n_success + list%n_failure write (u, "(A,I0)") "Success = ", list%n_success write (u, "(A,I0)") "Failure = ", list%n_failure write (u, "(A)") "*** End of test Summary ***" end subroutine test_results_write @ %def test_results_write @ Return true if all tests were successful (or no test). <>= procedure, private :: report => test_results_report +<>= + module subroutine test_results_report (list, success) + class(test_results_t), intent(in) :: list + logical, intent(out) :: success + end subroutine test_results_report <>= - subroutine test_results_report (list, success) + module subroutine test_results_report (list, success) class(test_results_t), intent(in) :: list logical, intent(out) :: success success = list%n_failure == 0 end subroutine test_results_report @ %def test_results_report @ Return the number of successful/failed/total tests.. <>= procedure :: get_n_pass => test_results_get_n_pass procedure :: get_n_fail => test_results_get_n_fail procedure :: get_n_total => test_results_get_n_total +<>= + module function test_results_get_n_pass (list) result (n) + class(test_results_t), intent(in) :: list + integer :: n + end function test_results_get_n_pass + module function test_results_get_n_fail (list) result (n) + class(test_results_t), intent(in) :: list + integer :: n + end function test_results_get_n_fail + module function test_results_get_n_total (list) result (n) + class(test_results_t), intent(in) :: list + integer :: n + end function test_results_get_n_total <>= - function test_results_get_n_pass (list) result (n) + module function test_results_get_n_pass (list) result (n) class(test_results_t), intent(in) :: list integer :: n n = list%n_success end function test_results_get_n_pass - function test_results_get_n_fail (list) result (n) + module function test_results_get_n_fail (list) result (n) class(test_results_t), intent(in) :: list integer :: n n = list%n_failure end function test_results_get_n_fail - function test_results_get_n_total (list) result (n) + module function test_results_get_n_total (list) result (n) class(test_results_t), intent(in) :: list integer :: n n = list%n_success + list%n_failure end function test_results_get_n_total @ %def test_results_get_n_pass @ %def test_results_get_n_fail @ %def test_results_get_n_total @ Delete the list. <>= procedure, private :: final => test_results_final +<>= + module subroutine test_results_final (list) + class(test_results_t), intent(inout) :: list + end subroutine test_results_final <>= - subroutine test_results_final (list) + module subroutine test_results_final (list) class(test_results_t), intent(inout) :: list type(test_result_t), pointer :: result do while (associated (list%first)) result => list%first list%first => result%next deallocate (result) end do list%last => null () list%n_success = 0 list%n_failure = 0 end subroutine test_results_final @ %def test_results_final @ \subsection{Wrapup} This will write results, report status, and finalize. This is the only method which we need to access from outside. <>= procedure :: wrapup => test_results_wrapup +<>= + module subroutine test_results_wrapup (list, u, success) + class(test_results_t), intent(inout) :: list + integer, intent(in) :: u + logical, intent(out), optional :: success + end subroutine test_results_wrapup <>= - subroutine test_results_wrapup (list, u, success) + module subroutine test_results_wrapup (list, u, success) class(test_results_t), intent(inout) :: list integer, intent(in) :: u logical, intent(out), optional :: success call list%write (u) if (present (success)) call list%report (success) call list%final () end subroutine test_results_wrapup @ %def test_results_wrapup @ \subsection{Tool for Unit Tests} This procedure takes a test routine as an argument. It runs the test, output directed to a temporary file. Then, it compares the file against a reference file. The test routine must take the output unit as argument. We export this abstract interface, so the test drivers can reference it for declaring the actual test routines. <>= public :: unit_test <>= abstract interface subroutine unit_test (u) integer, intent(in) :: u end subroutine unit_test end interface @ %def unit_test @ The test routine can print to screen and, optionally, to a logging unit. This driver runs the test given as an argument, directing its output to a scratch file which is then checked against a ref file. Only if it differs, it is copied to an err file. NB: We call [[fatal_force_crash]] to produce a deliberate crash with backtrace on any fatal error, if the runtime system does allow it. This is not normal behavior, but should be useful if something goes wrong. <>= public :: test +<>= + module subroutine test (test_proc, name, description, u_log, results) + procedure(unit_test) :: test_proc + character(*), intent(in) :: name + character(*), intent(in) :: description + integer, intent(in) :: u_log + type(test_results_t), intent(inout) :: results + end subroutine test <>= - subroutine test (test_proc, name, description, u_log, results) + module subroutine test (test_proc, name, description, u_log, results) procedure(unit_test) :: test_proc character(*), intent(in) :: name character(*), intent(in) :: description integer, intent(in) :: u_log type(test_results_t), intent(inout) :: results integer :: u_test logical :: success call start_test (u_log, name) u_test = free_unit () open (unit = u_test, status="scratch", action="readwrite") call fatal_force_crash () call test_proc (u_test) rewind (u_test) call compare_test_results (u_test, u_log, name, success) call results%add (name, description, success) close (u_test) end subroutine test @ %def test @ Startup message. We expose this, as well as the comparison routine, for re-use in external test programs. <>= public :: start_test +<>= + module subroutine start_test (u_log, name) + integer, intent(in) :: u_log + character(*), intent(in) :: name + end subroutine start_test <>= - subroutine start_test (u_log, name) + module subroutine start_test (u_log, name) integer, intent(in) :: u_log character(*), intent(in) :: name write (*, "(A)", advance="no") "Running test: " // name write (u_log, "(A)") "Test: " // name end subroutine start_test @ %def start_test @ Do the actual comparison, given an open I/O unit [[u_test]]. In case of failure, copy the test results to an error file. <>= public :: compare_test_results +<>= + module subroutine compare_test_results (u_test, u_log, name, success) + integer, intent(in) :: u_test + integer, intent(in) :: u_log + character(*), intent(in) :: name + logical, intent(out) :: success + end subroutine compare_test_results <>= - subroutine compare_test_results (u_test, u_log, name, success) + module subroutine compare_test_results (u_test, u_log, name, success) integer, intent(in) :: u_test integer, intent(in) :: u_log character(*), intent(in) :: name logical, intent(out) :: success integer :: u_ref, u_err logical :: exist character(256) :: buffer1, buffer2 integer :: iostat1, iostat2 inquire (file=ref_prefix//name//ref, exist=exist) if (exist) then open (newunit = u_ref, file=ref_prefix//name//ref, & status="old", action="read") COMPARE_FILES: do read (u_test, "(A)", iostat=iostat1) buffer1 read (u_ref, "(A)", iostat=iostat2) buffer2 if (iostat1 /= iostat2) then success = .false. exit COMPARE_FILES else if (iostat1 < 0) then success = .true. exit COMPARE_FILES else if (buffer1 /= buffer2) then success = .false. exit COMPARE_FILES end if end do COMPARE_FILES close (u_ref) else write (*, "(A)", advance="no") " ... no reference output available" write (u_log, "(A)") " No reference output available." success = .false. end if if (success) then write (*, "(A)") " ... success." write (u_log, "(A)") " Success." else write (*, "(A)") " ... failure. See: " // err_prefix//name//err write (u_log, "(A)") " Failure." rewind (u_test) open (newunit = u_err, file=err_prefix//name//err, & action="write", status="replace") WRITE_OUTPUT: do read (u_test, "(A)", end=1) buffer1 write (u_err, "(A)") trim (buffer1) end do WRITE_OUTPUT 1 close (u_err) end if end subroutine compare_test_results @ %def compare_test_results