Index: trunk/omega/tests/parameters_SM.f90 =================================================================== --- trunk/omega/tests/parameters_SM.f90 (revision 8756) +++ trunk/omega/tests/parameters_SM.f90 (revision 8757) @@ -1,152 +1,152 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! 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. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module parameters_sm use kinds use constants implicit none private real(default), dimension(27), public :: mass, width real(default), public :: as complex(default), public :: gs, igs, ig, unit real(default), public :: e, g, e_em real(default), public :: sinthw, costhw, sin2thw, tanthw real(default), public :: qelep, qeup, qedwn real(default), public :: ttop, tbot, tch, ttau, tw real(default), public :: ltop, lbot, lc, ltau, lw complex(default), public :: qlep, qup, qdwn, gcc, qw, & gzww, gwww, ghww, ghhww, ghzz, ghhzz, & ghbb, ghtt, ghcc, ghtautau, gh3, gh4, ghmm, & ! ghgaga, ghgaz, ghgg, ghmm, & iqw, igzww, igwww, gw4, gzzww, gazww, gaaww real(default), public :: vev complex(default), dimension(2), public :: & gncneu, gnclep, gncup, gncdwn public :: init_parameters, model_update_alpha_s real(default), parameter :: & GF = 1.16639E-5_default ! Fermi constant !!! This corresponds to 1/alpha = 137.03598949333 real(default), parameter :: & alpha = 1.0_default/137.03598949333_default complex(default), parameter :: & alphas = 0.1178_default ! Strong coupling constant (Z point) contains subroutine init_parameters mass(1:27) = 0 width(1:27) = 0 mass(3) = 0.095_default ! s-quark mass mass(4) = 1.2_default ! c-quark mass - ! mass(5) = 4.2_default ! b-quark mass - ! mass(6) = 173.1_default ! t-quark mass + mass(5) = 4.2_default ! b-quark mass + mass(6) = 173.1_default ! t-quark mass ! width(6) = 1.523_default ! t-quark width ! mass(11) = 0.000510997_default ! electron mass mass(13) = 0.105658389_default ! muon mass mass(15) = 1.77705_default ! tau-lepton mas mass(23) = 91.1882_default ! Z-boson mass width(23) = 2.443_default ! Z-boson width mass(24) = 80.419_default ! W-boson mass - width(24) = 2.049_default ! W-boson width + ! width(24) = 2.049_default ! W-boson width mass(25) = 200._default ! Higgs mass width(25) = 1.419_default ! Higgs width ttop = 4.0_default * mass(6)**2 / mass(25)**2 tbot = 4.0_default * mass(5)**2 / mass(25)**2 tch = 4.0_default * mass(4)**2 / mass(25)**2 ttau = 4.0_default * mass(15)**2 / mass(25)**2 tw = 4.0_default * mass(24)**2 / mass(25)**2 !ltop = 4.0_default * mass(6)**2 / mass(23)**2 !lbot = 4.0_default * mass(5)**2 / mass(23)**2 !lc = 4.0_default * mass(4)**2 / mass(23)**2 !ltau = 4.0_default * mass(15)**2 / mass(23)**2 !lw = 4.0_default * mass(24)**2 / mass(23)**2 e_em = sqrt(4.0_default * PI * alpha) vev = 1 / sqrt (sqrt (2.0_default) * GF) ! v (Higgs vev) costhw = mass(24) / mass(23) ! cos(theta-W) sinthw = sqrt (1.0_default-costhw**2) ! sin(theta-W) sin2thw = sinthw**2 tanthw = sinthw/costhw e = 2.0_default * sinthw * mass(24) / vev ! em-coupling (GF scheme) qelep = - 1 qeup = 2.0_default / 3.0_default qedwn = - 1.0_default / 3.0_default g = e / sinthw ig = cmplx (0.0_default, 1.0_default, kind=default) * g gcc = - g / 2 / sqrt (2.0_default) gncneu(1) = - g / 2 / costhw * ( + 0.5_default) gnclep(1) = - g / 2 / costhw * ( - 0.5_default - 2 * qelep * sin2thw) gncup(1) = - g / 2 / costhw * ( + 0.5_default - 2 * qeup * sin2thw) gncdwn(1) = - g / 2 / costhw * ( - 0.5_default - 2 * qedwn * sin2thw) gncneu(2) = - g / 2 / costhw * ( + 0.5_default) gnclep(2) = - g / 2 / costhw * ( - 0.5_default) gncup(2) = - g / 2 / costhw * ( + 0.5_default) gncdwn(2) = - g / 2 / costhw * ( - 0.5_default) qlep = - e * qelep qup = - e * qeup qdwn = - e * qedwn qw = e iqw = (0,1)*qw gzww = g * costhw igzww = (0,1)*gzww gwww = g igwww = (0,1)*gwww gw4 = gwww**2 gzzww = gzww**2 gazww = gzww * qw gaaww = qw**2 ghww = mass(24) * g ghhww = g**2 / 2.0_default ghzz = mass(23) * g / costhw ghhzz = g**2 / 2.0_default / costhw**2 ghtt = - mass(6) / vev ghbb = - mass(5) / vev ghcc = - mass(4) / vev ghtautau = - mass(15) / vev ghmm = - mass(13) / vev gh3 = - 3 * mass(25)**2 / vev gh4 = - 3 * mass(25)**2 / vev**2 !!! Color flow basis, divide by sqrt(2) gs = sqrt(2.0_default*PI*alphas) igs = cmplx (0.0_default, 1.0_default, kind=default) * gs unit = 1.0_default end subroutine init_parameters subroutine model_update_alpha_s (alpha_s) real(default), intent(in) :: alpha_s gs = sqrt(2.0_default*PI*alpha_s) igs = cmplx (0.0_default, 1.0_default, kind=default) * gs end subroutine model_update_alpha_s end module parameters_sm Index: trunk/omega/tests/ward_lib.f90 =================================================================== --- trunk/omega/tests/ward_lib.f90 (revision 8756) +++ trunk/omega/tests/ward_lib.f90 (revision 8757) @@ -1,325 +1,448 @@ ! ward_lib.f90 -- check On Shell Ward Identities in O'Mega !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Copyright (C) 1999-2021 by ! Wolfgang Kilian ! Thorsten Ohl ! Juergen Reuter -! Christian Speckner +! 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. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module ward_lib ! use ieee_arithmetic use kinds use constants use tao_random_numbers use omega95 use omega_interface use omega_testtools implicit none private public :: check + + interface boost + module procedure boost_one, boost_many + end interface + contains elemental function ieee_is_nan (x) result (yorn) logical :: yorn real (kind=default), intent(in) :: x yorn = (x /= x) end function ieee_is_nan - subroutine check (physical, unphysical, roots, threshold, n, failures, attempts, seed) + subroutine check (physical, unphysical, roots, m1, m2, m3, m4, & + threshold, n, failures, attempts, seed) type(omega_procedures), intent(in) :: physical, unphysical real(kind=default), intent(in) :: roots, threshold + real(kind=default), intent(in), optional :: m1, m2, m3, m4 integer, intent(in) :: n integer, intent(out) :: failures, attempts integer, intent(in), optional :: seed logical :: match, passed integer :: n_out, n_flv, n_hel, n_col integer :: i, i_flv, i_hel, i_col ! integer :: i_prt integer, dimension(:,:), allocatable :: spin_states_phys, spin_states_unphys real(kind=default), dimension(:,:), allocatable :: p + real(kind=default), dimension(:), allocatable :: mass complex(kind=default), dimension(:), allocatable :: a character(len=80) :: msg complex(kind=default) :: wi - real(kind=default) :: a_avg + real(kind=default) :: a_avg, mass_thr failures = 0 attempts = 0 + !!! We take a numerical threshold which is smaller than the electron mass + mass_thr = 0.0001_default call quantum_numbers (physical, unphysical, n_out, n_flv, n_hel, n_col, match) if (.not.match) then failures = 1 return end if if (n_out <= 0) then print *, "no outgoing particles" failures = 1 return end if if (n_flv <= 0) then print *, "no allowed flavor combinations" failures = 1 return end if if (n_hel <= 0) then print *, "no allowed helicity combinations" failures = 1 return end if if (n_col <= 0) then print *, "no allowed color flows" failures = 1 return end if if (present (seed)) then call tao_random_seed (seed) end if call physical%reset_helicity_selection (-1.0_default, -1) call unphysical%reset_helicity_selection (-1.0_default, -1) allocate (p(0:3,2+n_out)) allocate (a(n_hel)) allocate (spin_states_phys(2+n_out,n_hel)) allocate (spin_states_unphys(2+n_out,unphysical%number_spin_states())) call physical%spin_states(spin_states_phys) call unphysical%spin_states(spin_states_unphys) - call beams (ROOTS, 0.0_default, 0.0_default, p(:,1), p(:,2)) + call beams (ROOTS, m1, m2, p(:,1), p(:,2)) do i = 1, N - call massless_isotropic_decay (ROOTS, p(:,3:)) + if (m3 > mass_thr .or. m4 > mass_thr) then + allocate (mass (n_out)) + mass = 0 + mass(1) = m3 + mass(2) = m4 + call massive_decay (ROOTS, mass, p(:,3:)) + deallocate (mass) + else + call massless_isotropic_decay (ROOTS, p(:,3:)) + end if call physical%new_event (p) call unphysical%new_event (p) do i_flv = 1, n_flv do i_col = 1, n_col do i_hel = 1, n_hel a(i_hel) = physical%get_amplitude (i_flv, i_hel, i_col) ! do i_prt = 1, (2+n_out) ! if (spin_states_phys(i_prt,i_hel).eq.0) then ! a(i_hel) = 0.0_default ! exit ! end if ! end do end do a_avg = sum (abs (a)) / n_hel !write (*, "(1X,'a_avg=',E15.5)") a_avg if (.not. ieee_is_nan (a_avg)) then if (a_avg > 0) then do i_hel = 1, n_hel / 2 ! do i_hel = 1, size(spin_states_unphys,dim=2) wi = unphysical%get_amplitude (i_flv, i_hel, i_col) ! do i_prt = 1, (2+n_out) ! if (spin_states_unphys(i_prt,i_hel).eq.0) then ! wi = 0.0_default ! exit ! end if ! end do attempts = attempts + 1 write (msg, "(1X,'evt=',I5,', flv=',I3,', col=',I3,', hel=',I3)") & i, i_flv, i_col, i_hel passed = .true. call expect_zero (wi, a_avg, trim(msg), passed, & quiet=.true., threshold = threshold) if (.not.passed) then failures = failures + 1 end if end do else ! write (*, "(1X,'evt=',I5,', flv=',I3,', col=',I3,': ', A)") & ! i, i_flv, i_col, "skipped: physical amplitude vanishes" end if else write (*, "(1X,'evt=',I5,', flv=',I3,', col=',I3,': ', A)") & i, i_flv, i_col, "physical amplitude NaN" attempts = attempts + 1 failures = failures + 1 end if end do end do end do deallocate (p) deallocate (a) deallocate (spin_states_phys) deallocate (spin_states_unphys) end subroutine check subroutine quantum_numbers (physical, unphysical, n_out, n_flv, n_hel, n_col, match) type(omega_procedures), intent(in) :: physical, unphysical integer, intent(out) :: n_out, n_flv, n_hel, n_col logical, intent(out) :: match integer, dimension(:,:), allocatable :: & physical_flavor_states, unphysical_flavor_states, & physical_spin_states, unphysical_spin_states integer, dimension(:,:,:), allocatable :: & physical_color_flows, unphysical_color_flows logical, dimension(:,:), allocatable :: & physical_ghost_flags, unphysical_ghost_flags type(omega_color_factor), dimension(:), allocatable :: & physical_color_factors, unphysical_color_factors integer :: n_in, n_prt, n_cix, n_cfs n_in = physical%number_particles_in () n_out = physical%number_particles_out () n_prt = n_in + n_out n_flv = physical%number_flavor_states () n_hel = physical%number_spin_states () n_cix = physical%number_color_indices () n_col = physical%number_color_flows () n_cfs = physical%number_color_factors () match = .true. if (unphysical%number_particles_in () .ne. n_in) then print *, "#particles_in don't match!" match = .false. end if if (unphysical%number_particles_out () .ne. n_out) then print *, "#particles_out don't match!" match = .false. end if if (unphysical%number_flavor_states () .ne. n_flv) then print *, "#flavor_states don't match!" match = .false. end if if (unphysical%number_spin_states () .ne. n_hel/2) then print *, "#spin_states don't match!" ! match = .false. end if if (unphysical%number_color_indices () .ne. n_cix) then print *, "#color_indices don't match!" match = .false. end if if (unphysical%number_color_flows () .ne. n_col) then print *, "#color_flows don't match!" match = .false. end if if (unphysical%number_color_factors () .ne. n_cfs) then print *, "#color_factors don't match!" match = .false. end if if (match) then allocate (physical_flavor_states(n_prt,n_flv), unphysical_flavor_states(n_prt,n_flv)) allocate (physical_spin_states(n_prt,n_hel), unphysical_spin_states(n_prt,n_hel/2)) allocate (physical_color_flows(n_cix,n_prt,n_col), & unphysical_color_flows(n_cix,n_prt,n_col)) allocate (physical_ghost_flags(n_prt,n_col), unphysical_ghost_flags(n_prt,n_col)) allocate (physical_color_factors(n_cfs), unphysical_color_factors(n_cfs)) call physical%flavor_states (physical_flavor_states) call unphysical%flavor_states (unphysical_flavor_states) call physical%spin_states (physical_spin_states) call unphysical%spin_states (unphysical_spin_states) call physical%color_flows (physical_color_flows, physical_ghost_flags) call unphysical%color_flows (unphysical_color_flows, unphysical_ghost_flags) call physical%color_factors (physical_color_factors) call unphysical%color_factors (unphysical_color_factors) if (any (physical_flavor_states .ne. unphysical_flavor_states)) then print *, "flavor states don't match!" print *, "CAVEAT: this might be due to simple reordering!" match = .false. end if ! if (any (physical_spin_states .ne. unphysical_spin_states)) then ! print *, "spin states don't match!" ! print *, "CAVEAT: this might be due to simple reordering!" ! match = .false. ! end if if (any (physical_color_flows .ne. unphysical_color_flows)) then print *, "color flows don't match!" print *, "CAVEAT: this might be due to simple reordering!" match = .false. end if if (any (physical_ghost_flags .neqv. unphysical_ghost_flags)) then print *, "ghost flags don't match!" print *, "CAVEAT: this might be due to simple reordering!" match = .false. end if if (any (.not. color_factors_equal (physical_color_factors, & unphysical_color_factors))) then print *, "color_factors don't match!" print *, "CAVEAT: this might be due to simple reordering!" match = .false. end if deallocate (physical_flavor_states, unphysical_flavor_states) deallocate (physical_spin_states, unphysical_spin_states) deallocate (physical_color_flows, unphysical_color_flows) deallocate (physical_ghost_flags, unphysical_ghost_flags) deallocate (physical_color_factors, unphysical_color_factors) end if end subroutine quantum_numbers elemental function color_factors_equal (cf1, cf2) result (eq) logical :: eq type(omega_color_factor), intent(in) :: cf1, cf2 eq = (cf1%i1 .eq. cf2%i1) .and. (cf1%i2 .eq. cf2%i2) .and. (cf1%factor .eq. cf2%factor) end function color_factors_equal pure function dot (p, q) result (pq) real(kind=default), dimension(0:), intent(in) :: p, q real(kind=default) :: pq pq = p(0)*q(0) - dot_product (p(1:), q(1:)) end function dot pure subroutine beams (roots, m1, m2, p1, p2) real(kind=default), intent(in) :: roots, m1, m2 real(kind=default), dimension(0:), intent(out) :: p1, p2 real(kind=default) :: m12, m22 m12 = m1**2 m22 = m2**2 p1(0) = (roots**2 + m12 - m22) / (2*roots) p1(1:2) = 0 p1(3) = sqrt (p1(0)**2 - m12) p2(0) = roots - p1(0) p2(1:3) = - p1(1:3) end subroutine beams ! The massless RAMBO algorithm subroutine massless_isotropic_decay (roots, p) real(kind=default), intent(in) :: roots real(kind=default), dimension(0:,:), intent(out) :: p real(kind=default), dimension(0:3,size(p,dim=2)) :: q real(kind=default), dimension(0:3) :: qsum real(kind=double), dimension(4) :: ran_double real(kind=default), dimension(4) :: ran real(kind=default) :: c, s, f, qabs, x, r, z integer :: k ! Generate isotropic null vectors do k = 1, size (p, dim = 2) ! if default is not double or single, we can't use ! tao_random_number directly ... call tao_random_number (ran_double) ran = ran_double ! generate a x*exp(-x) distribution for q(0,k) q(0,k)= -log(ran(1)*ran(2)) c = 2*ran(3)-1 f = 2*PI*ran(4) s = sqrt(1-c*c) q(2,k) = q(0,k)*s*sin(f) q(3,k) = q(0,k)*s*cos(f) q(1,k) = q(0,k)*c enddo ! Boost and rescale the vectors qsum = sum (q, dim = 2) qabs = sqrt (dot (qsum, qsum)) x = roots/qabs do k = 1, size (p, dim = 2) r = dot (q(0:,k), qsum) / qabs z = (q(0,k)+r)/(qsum(0)+qabs) p(1:3,k) = x*(q(1:3,k)-qsum(1:3)*z) p(0,k) = x*r enddo end subroutine massless_isotropic_decay + pure function mass2 (p) result (m2) + real (kind=default), dimension(0:), intent(in) :: p + real (kind=default) :: m2 + m2 = p(0)*p(0) - p(1)*p(1) - p(2)*p(2) - p(3)*p(3) + end function mass2 + + pure subroutine boost_one (v, p, q) + real (kind=default), dimension(0:), intent(in) :: v, p + real (kind=default), dimension(0:), intent(out) :: q + q(0) = dot_product (p, v) + q(1:3) = p(1:3) & + + v(1:3) * (p(0) + dot_product (p(1:3), v(1:3)) / (1 + v(0))) + end subroutine boost_one + + pure subroutine boost_many (v, p, q) + real (kind=default), dimension(0:), intent(in) :: v + real (kind=default), dimension(0:,:), intent(in) :: p + real (kind=default), dimension(0:,:), intent(out) :: q + integer :: k + do k = 1, size (p, dim = 2) + call boost_one (v, p(:,k), q(:,k)) + enddo + end subroutine boost_many + + !!! The massive RAMBO algorithm (not reweighted, therefore not isotropic) + subroutine massive_decay (roots, m, p) + real(kind=default), intent(in) :: roots + real(kind=default), dimension(:), intent(in) :: m + real(kind=default), dimension(0:,:), intent(out) :: p + real(kind=default), dimension(0:3,size(p,dim=2)) :: q + real(kind=default), dimension(size(p,dim=2)) :: p2, m2, p0 + real(kind=default), dimension(0:3) :: qsum + real(kind=double), dimension(2) :: ran_double + real(kind=default), dimension(2) :: ran + real(kind=default) :: c, s, f, qq + real(kind=default) :: w, a, xu, u, umax, xv, v, vmax, x + real(kind=default) :: xi, delta + integer :: k, i + if (sum(m) > roots) then + print *, "no solution: sum(m) > roots" + p = 0 + return + end if + m2 = m*m + ! Generate isotropic massive vectors + w = 1 + do k = 1, size (p, dim = 2) + ! Kinderman/Monahan (a la Kleiss/Sterling) + a = 2 * m(k) / w + xu = 0.5 * (1 - a + sqrt (1 + a*a)) + xv = 0.5 * (3 - a + sqrt (9 + 4*a + a*a)) + umax = exp (-0.5*xu) * sqrt (sqrt (xu*xu + a*xu)) + vmax = xv * exp (-0.5*xv) * sqrt (sqrt (xv*xv + a*xv)) + rejection: do + call tao_random_number (ran_double) + ran = ran_double + u = ran(1) * umax + v = ran(2) * vmax + x = v / u + if (u*u < exp(-x) * sqrt (x*x + a*x)) then + qq = m(k) + w*x + exit rejection + end if + end do rejection + call tao_random_number (ran_double) + ran = ran_double + c = 2*ran(1) - 1 + !!! select case (k) + !!! case (1,3) + !!! c = 1 - 0.0000002*ran(1) + !!! case (2,4) + !!! c = 0.0000002*ran(1) - 1 + !!! end select + f = 2*PI*ran(2) + s = sqrt (1 - c*c) + q(0,k) = sqrt (qq*qq + m2(k)) + q(1,k) = qq * s * sin(f) + q(2,k) = qq * s * cos(f) + q(3,k) = qq * c + enddo + ! Boost the vectors to the common rest frame + qsum = sum (q, dim = 2) + call boost ((/ qsum(0), - qsum(1:3) /) / sqrt (mass2 (qsum)), q, p) + ! rescale momenta + do k = 1, size (p, dim = 2) + p2(k) = dot_product (p(1:3,k), p(1:3,k)) + end do + i = 1 + xi = 1 + find_xi: do + p0 = sqrt (xi*xi*p2 + m2) + delta = sum (p0) - roots + if ((i > 100) .or. (abs (delta) <= 10 * epsilon (roots))) then + exit find_xi + end if + ! Newton / Ralphson iteration + xi = xi - delta / (xi * sum (p2 / p0)) + i = i + 1 + end do find_xi + p(0,:) = p0 + p(1:3,:) = xi * p(1:3,:) + end subroutine massive_decay + end module ward_lib Index: trunk/omega/tests/ward_driver_UFO.sh =================================================================== --- trunk/omega/tests/ward_driver_UFO.sh (revision 8756) +++ trunk/omega/tests/ward_driver_UFO.sh (revision 8757) @@ -1,159 +1,161 @@ #! /bin/sh # ward_driver_UFO.sh -- ######################################################################## omega="$1" shift models="sm_ufo" modules="" ######################################################################## while read module threshold n roots model unphysical mode process; do case $module in '#'*) # skip comments ;; '') # skip empty lines ;; '!'*) break ;; *) ######################################################################## modules="$modules $module" eval threshold_$module=$threshold eval n_$module=$n eval roots_$module=$roots eval process_$module="'$process'" ######################################################################## ######################################################################## # echo "running $omega_bin -$mode '$process'" 1>&2 $omega "$@" -model:exec \ -target:parameter_module parameters_sm_ufo \ -target:module amplitude_ward_ufo_physical_$module \ -$mode "$process" 2>/dev/null $omega "$@" -model:exec \ -target:parameter_module parameters_sm_ufo \ -target:module amplitude_ward_ufo_unphysical_$module \ -$mode "$process" -unphysical $unphysical 2>/dev/null ;; esac done ######################################################################## for module in $modules; do for mode in physical unphysical; do cat < number_particles_in p%number_particles_out => number_particles_out p%number_spin_states => number_spin_states p%spin_states => spin_states p%number_flavor_states => number_flavor_states p%flavor_states => flavor_states p%number_color_indices => number_color_indices p%number_color_flows => number_color_flows p%color_flows => color_flows p%number_color_factors => number_color_factors p%color_factors => color_factors p%color_sum => color_sum p%new_event => new_event p%reset_helicity_selection => reset_helicity_selection p%is_allowed => is_allowed p%get_amplitude => get_amplitude end function load end module interface_ward_ufo_${mode}_${module} EOF done done ######################################################################## cat < load EOF done done for model in $models; do cat < setup_parameters EOF done cat < 0) then print *, failures, " failures in ", attempts, " attempts" failed_processes = failed_processes + 1 end if EOF done cat < 0) then print *, failed_processes, " failed processes in ", attempted_processes, " attempts" stop 1 end if end program ward_ufo_driver EOF exit 0 Index: trunk/omega/tests/ward_identities_long.list =================================================================== --- trunk/omega/tests/ward_identities_long.list (revision 8756) +++ trunk/omega/tests/ward_identities_long.list (revision 8757) @@ -1,10 +1,10 @@ # ward_identities_long.list -- # ---------------------------------------------------------------------- -# thr n roots model i process ... +# thr n roots m1 m2 model i process ... # ---------------------------------------------------------------------- -uuggg 0.7 50 1000 QCD 3 scatter u ubar -> gl gl gl -uggggg 0.7 10 1000 QCD 3 scatter u ubar -> gl gl gl gl -ggggg 0.7 50 1000 QCD 3 scatter gl gl -> gl gl gl -gggggg 0.7 10 1000 QCD 3 scatter gl gl -> gl gl gl gl +uuggg 0.7 50 1000 0.000 0.000 QCD 3 scatter u ubar -> gl gl gl +uggggg 0.7 10 1000 0.000 0.000 QCD 3 scatter u ubar -> gl gl gl gl +ggggg 0.7 50 1000 0.000 0.000 QCD 3 scatter gl gl -> gl gl gl +gggggg 0.7 10 1000 0.000 0.000 QCD 3 scatter gl gl -> gl gl gl gl Index: trunk/omega/tests/ward_identities_fail.list =================================================================== --- trunk/omega/tests/ward_identities_fail.list (revision 8756) +++ trunk/omega/tests/ward_identities_fail.list (revision 8757) @@ -1,21 +1,21 @@ # ward_identities_fail.list -- # ---------------------------------------------------------------------- -# thr n roots model i process ... +# thr n roots m1 m2 model i process ... # ---------------------------------------------------------------------- # # Needs a mass for the Higgs: -#gg1 0.7 1000 1000 SYM 3 scatter phi g1 -> g1 +#gg1 0.7 1000 1000 0.000 0.000 SYM 3 scatter phi g1 -> g1 # # Works without a Hgg coupling (the required Hggg coupling is missing) -#gggg1 0.7 1 1000 SYM 3 scatter g1 g1 -> g1 g1 g1 -#ggggg1 0.7 1 1000 SYM 3 scatter g1 g1 -> g1 g1 g1 g1 +#gggg1 0.7 1 1000 0.000 0.000 SYM 3 scatter g1 g1 -> g1 g1 g1 +#ggggg1 0.7 1 1000 0.000 0.000 SYM 3 scatter g1 g1 -> g1 g1 g1 g1 # # The required Hggg coupling is missing -#ggp1 0.7 1 1000 SYM 3 scatter g1 g1 -> g1 phi +#ggp1 0.7 1 1000 0.000 0.000 SYM 3 scatter g1 g1 -> g1 phi # # Without Hgg, SYM should work just as well as QCD -qqggg 0.7 500 1000 QCD 3 scatter u ubar -> gl gl gl -qqggg1 0.7 500 1000 SYM 3 scatter q1 Q1 -> g1 g1 g1 -ssggg1 0.7 500 1000 SYM 3 scatter sq1 sQ1 -> g1 g1 g1 -ssgss1 0.7 500 1000 SYM 3 scatter sq1 sQ1 -> g1 sq1 sQ1 +qqggg 0.7 500 1000 0.000 0.000 QCD 3 scatter u ubar -> gl gl gl +qqggg1 0.7 500 1000 0.000 0.000 SYM 3 scatter q1 Q1 -> g1 g1 g1 +ssggg1 0.7 500 1000 0.000 0.000 SYM 3 scatter sq1 sQ1 -> g1 g1 g1 +ssgss1 0.7 500 1000 0.000 0.000 SYM 3 scatter sq1 sQ1 -> g1 sq1 sQ1 Index: trunk/omega/tests/ward_identities.list =================================================================== --- trunk/omega/tests/ward_identities.list (revision 8756) +++ trunk/omega/tests/ward_identities.list (revision 8757) @@ -1,92 +1,102 @@ # ward_identities.list -- +# note: massive states need more relaxed thresholds (more numerical noise) +# in general avg. thresholds have to be lowered from 0.35 to 0.15 # ---------------------------------------------------------------------- -# thr n roots model i process ... +# thr n roots m1 m2 m3 m4 model i process ... # ---------------------------------------------------------------------- -eeaa 0.75 100 1000 QED 3 scatter e+ e- -> A A -eeaaa 0.75 100 1000 QED 3 scatter e+ e- -> A A A -eeaee 0.35 100 1000 QED 3 scatter e+ e- -> A e+ e- -aeae 0.35 100 1000 QED 3 scatter A e+ -> A e+ -aeaea 0.35 100 1000 QED 3 scatter A e+ -> A e+ A -# -uugg 0.70 100 1000 QCD 3 scatter u ubar -> gl gl -gggg 0.70 100 1000 QCD 3 scatter gl gl -> gl gl -uuguu 0.35 50 1000 QCD 3 scatter u ubar -> gl u ubar -gugu 0.30 100 1000 QCD 3 scatter gl u -> gl u -# numerical noise improves w/quad: -uugguu 0.25 50 1000 QCD 3 scatter u ubar -> gl gl u ubar -# -qqgg1 0.70 100 1000 SYM 3 scatter q1 Q1 -> g1 g1 -ssgg1 0.70 100 1000 SYM 3 scatter sq1 sQ1 -> g1 g1 -ppgg1 0.70 100 1000 SYM 3 scatter phi phi -> g1 g1 -gggg1 0.70 100 1000 SYM 3 scatter g1 g1 -> g1 g1 -# -ggtt 0.35 100 1000 SM 1 scatter gl gl -> t tbar -gtgt 0.35 100 1000 SM 3 scatter gl t -> gl t -ggtt_a 0.35 100 1000 SM_top_anom 1 scatter gl gl -> t tbar -gtgt_a 0.35 100 1000 SM_top_anom 3 scatter gl t -> gl t -# -aatt 0.30 100 1000 SM 3 scatter t tbar -> A A -ttaa 0.30 100 1000 SM 1 scatter A A -> t tbar -tata 0.30 100 1000 SM 3 scatter A t -> A t -atat 0.30 100 1000 SM 3 scatter A tbar -> A tbar -aatt_a 0.30 100 1000 SM_top_anom 3 scatter t tbar -> A A -ttaa_a 0.30 100 1000 SM_top_anom 1 scatter A A -> t tbar -tata_a 0.30 100 1000 SM_top_anom 3 scatter A t -> A t -atat_a 0.30 100 1000 SM_top_anom 3 scatter A tbar -> A tbar -# -#ttaz 0.30 100 1000 SM 3 scatter t tbar -> A Z -#aztt 0.30 100 1000 SM 1 scatter A Z -> t tbar -#tzta 0.30 100 1000 SM 3 scatter Z t -> A t -#ztat 0.30 100 1000 SM 3 scatter Z tbar -> A tbar -#ttaz_a 0.30 100 1000 SM_top_anom 3 scatter t tbar -> A Z -#aztt_a 0.30 100 1000 SM_top_anom 1 scatter A Z -> t tbar -#tzta_a 0.30 100 1000 SM_top_anom 3 scatter Z t -> A t -#ztat_a 0.30 100 1000 SM_top_anom 3 scatter Z tbar -> A tbar -# -#ttza 0.30 100 1000 SM 3 scatter t tbar -> Z A -#zatt 0.30 100 1000 SM 1 scatter Z A -> t tbar -#tatz 0.30 100 1000 SM 3 scatter A t -> Z t -#atzt 0.30 100 1000 SM 3 scatter A tbar -> Z tbar -#ttza_a 0.30 100 1000 SM_top_anom 3 scatter t tbar -> Z A -#zatt_a 0.30 100 1000 SM_top_anom 1 scatter Z A -> t tbar -#tatz_a 0.30 100 1000 SM_top_anom 3 scatter A t -> Z t -#atzt_a 0.30 100 1000 SM_top_anom 3 scatter A tbar -> Z tbar -# -#btwz 0.35 100 1000 SM 4 scatter t bbar -> W+ Z -#wbzt 0.35 100 1000 SM 3 scatter W+ b -> Z t -#btwz_a 0.35 100 1000 SM_top_anom 4 scatter t bbar -> W+ Z -#wbzt_a 0.35 100 1000 SM_top_anom 3 scatter W+ b -> Z t -# -#tbwz 0.35 100 1000 SM 4 scatter b tbar -> W- Z -#wtzb 0.35 100 1000 SM 3 scatter W- t -> Z b -#tbwz_a 0.35 100 1000 SM_top_anom 4 scatter b tbar -> W- Z -#wtzb_a 0.35 100 1000 SM_top_anom 3 scatter W- t -> Z b -# -#btwa 0.35 100 1000 SM 4 scatter t bbar -> W+ A -#wbat 0.35 100 1000 SM 3 scatter W+ b -> A t -#btwa_a 0.35 100 1000 SM_top_anom 4 scatter t bbar -> W+ A -#wbat_a 0.35 100 1000 SM_top_anom 3 scatter W+ b -> A t -# -#tbaw 0.35 100 1000 SM 3 scatter b tbar -> A W- -#wtab 0.35 100 1000 SM 3 scatter W- t -> A b -#tbaw_a 0.35 100 1000 SM_top_anom 3 scatter b tbar -> A W- -#wtab_a 0.35 100 1000 SM_top_anom 3 scatter W- t -> A b -# -#ttzz 0.30 100 1000 SM 3 scatter t tbar -> Z Z -#zztt 0.30 100 1000 SM 1 scatter Z Z -> t tbar -#tztz 0.30 100 1000 SM 3 scatter Z t -> Z t -#ztzt 0.30 100 1000 SM 3 scatter Z tbar -> Z tbar -#ttzz_a 0.30 100 1000 SM_top_anom 3 scatter t tbar -> Z Z -#zztt_a 0.30 100 1000 SM_top_anom 1 scatter Z Z -> t tbar -#tztz_a 0.30 100 1000 SM_top_anom 3 scatter Z t -> Z t -#ztzt_a 0.30 100 1000 SM_top_anom 3 scatter Z tbar -> Z tbar -# -#wwtt 0.35 100 1000 SM 3 scatter t tbar -> W+ W- -#wtwt 0.35 100 1000 SM 3 scatter W+ t -> W+ t -#wwtt_a 0.35 100 1000 SM_top_anom 3 scatter t tbar -> W+ W- -#wtwt_a 0.35 100 1000 SM_top_anom 3 scatter W+ t -> W+ t -# -#wwbb 0.35 100 1000 SM 3 scatter b bbar -> W+ W- -#wbwb 0.35 100 1000 SM 3 scatter W+ b -> W+ b -#wwbb_a 0.35 100 1000 SM_top_anom 3 scatter b bbar -> W+ W- -#wbwb_a 0.35 100 1000 SM_top_anom 3 scatter W+ b -> W+ b +eeaa 0.75 100 1000 0.000 0.000 0.000 0.000 QED 3 scatter e+ e- -> A A +eeaaa 0.75 100 1000 0.000 0.000 0.000 0.000 QED 3 scatter e+ e- -> A A A +eeaee 0.35 100 1000 0.000 0.000 0.000 0.000 QED 3 scatter e+ e- -> A e+ e- +aeae 0.35 100 1000 0.000 0.000 0.000 0.000 QED 3 scatter A e+ -> A e+ +aeaea 0.35 100 1000 0.000 0.000 0.000 0.000 QED 3 scatter A e+ -> A e+ A +# +uugg 0.70 100 1000 0.000 0.000 0.000 0.000 QCD 3 scatter u ubar -> gl gl +gggg 0.70 100 1000 0.000 0.000 0.000 0.000 QCD 3 scatter gl gl -> gl gl +uuguu 0.35 50 1000 0.000 0.000 0.000 0.000 QCD 3 scatter u ubar -> gl u ubar +gugu 0.30 100 1000 0.000 0.000 0.000 0.000 QCD 3 scatter gl u -> gl u +# numerical noise improves w/quad: +uugguu 0.25 50 1000 0.000 0.000 0.000 0.000 QCD 3 scatter u ubar -> gl gl u ubar +# +qqgg1 0.70 100 1000 0.000 0.000 0.000 0.000 SYM 3 scatter q1 Q1 -> g1 g1 +ssgg1 0.70 100 1000 0.000 0.000 0.000 0.000 SYM 3 scatter sq1 sQ1 -> g1 g1 +ppgg1 0.70 100 1000 0.000 0.000 0.000 0.000 SYM 3 scatter phi phi -> g1 g1 +gggg1 0.70 100 1000 0.000 0.000 0.000 0.000 SYM 3 scatter g1 g1 -> g1 g1 +# +ggtt 0.20 100 1000 0.000 0.000 173.100 173.100 SM 1 scatter gl gl -> t tbar +gtgt 0.15 100 1000 0.000 173.100 0.000 173.100 SM 3 scatter gl t -> gl t +ggtt_a 0.25 100 1000 0.000 0.000 173.100 173.100 SM_top_anom 1 scatter gl gl -> t tbar +gtgt_a 0.25 100 1000 0.000 173.100 0.000 173.100 SM_top_anom 3 scatter gl t -> gl t +# +ttaa 0.20 100 1000 173.100 173.100 0.000 0.000 SM 3 scatter t tbar -> A A +aatt 0.20 100 1000 0.000 0.000 173.100 173.100 SM 1 scatter A A -> t tbar +tata 0.20 100 1000 0.000 173.100 0.000 173.100 SM 3 scatter A t -> A t +atat 0.20 100 1000 0.000 173.100 0.000 173.100 SM 3 scatter A tbar -> A tbar +ttaa_a 0.30 100 1000 173.100 173.100 0.000 0.000 SM_top_anom 3 scatter t tbar -> A A +aatt_a 0.30 100 1000 0.000 0.000 173.100 173.100 SM_top_anom 1 scatter A A -> t tbar +tata_a 0.30 100 1000 0.000 173.100 0.000 173.100 SM_top_anom 3 scatter A t -> A t +atat_a 0.30 100 1000 0.000 173.100 0.000 173.100 SM_top_anom 3 scatter A tbar -> A tbar +# +ttaz 0.20 100 1000 173.100 173.100 0.000 91.1882 SM 3 scatter t tbar -> A Z +aztt 0.20 100 1000 0.000 91.1882 173.100 173.100 SM 1 scatter A Z -> t tbar +tzta 0.15 100 1000 91.1882 173.100 0.000 173.100 SM 3 scatter Z t -> A t +ztat 0.15 100 1000 91.1882 173.100 0.000 173.100 SM 3 scatter Z tbar -> A tbar +ttaz_a 0.30 100 1000 173.100 173.100 0.000 91.1882 SM_top_anom 3 scatter t tbar -> A Z +aztt_a 0.30 100 1000 0.000 91.1882 173.10 173.100 SM_top_anom 1 scatter A Z -> t tbar +tzta_a 0.30 100 1000 91.1882 173.100 0.000 173.100 SM_top_anom 3 scatter Z t -> A t +ztat_a 0.30 100 1000 91.1882 173.100 0.000 173.100 SM_top_anom 3 scatter Z tbar -> A tbar +# Ward identities for the Z boson +#ttza 0.30 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter t tbar -> Z A +#zatt 0.30 100 1000 0.000 0.000 0.000 0.000 SM 1 scatter Z A -> t tbar +#tatz 0.30 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter A t -> Z t +#atzt 0.30 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter A tbar -> Z tbar +#ttza_a 0.30 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter t tbar -> Z A +#zatt_a 0.30 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 1 scatter Z A -> t tbar +#tatz_a 0.30 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter A t -> Z t +#atzt_a 0.30 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter A tbar -> Z tbar +# Ward identities for the Z boson +#btwz 0.35 100 1000 0.000 0.000 0.000 0.000 SM 4 scatter t bbar -> W+ Z +#wbzt 0.35 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter W+ b -> Z t +#btwz_a 0.35 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 4 scatter t bbar -> W+ Z +#wbzt_a 0.35 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter W+ b -> Z t +# Ward identities for the Z boson +#tbwz 0.35 100 1000 0.000 0.000 0.000 0.000 SM 4 scatter b tbar -> W- Z +#wtzb 0.35 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter W- t -> Z b +#tbwz_a 0.35 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 4 scatter b tbar -> W- Z +#wtzb_a 0.35 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter W- t -> Z b +# asymmetric processes need lower thresholds +btwa 0.10 100 1000 173.100 4.200 80.419 0.000 SM 4 scatter t bbar -> W+ A +wbat 0.20 100 1000 80.419 4.200 0.000 173.100 SM 3 scatter W+ b -> A t +btwa_a 0.10 100 1000 173.100 4.200 80.419 0.000 SM_top_anom 4 scatter t bbar -> W+ A +wbat_a 0.30 100 1000 80.419 4.200 0.000 173.100 SM_top_anom 3 scatter W+ b -> A t +# +tbaw 0.10 100 1000 4.200 173.100 0.000 80.419 SM 3 scatter b tbar -> A W- +wtab 0.20 100 1000 80.419 173.100 0.000 4.200 SM 3 scatter W- t -> A b +tbaw_a 0.10 100 1000 4.200 173.100 0.000 80.419 SM_top_anom 3 scatter b tbar -> A W- +wtab_a 0.30 100 1000 80.419 173.100 0.000 4.200 SM_top_anom 3 scatter W- t -> A b +# Ward identities for Z bosons +#ttzz 0.30 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter t tbar -> Z Z +#zztt 0.30 100 1000 0.000 0.000 0.000 0.000 SM 1 scatter Z Z -> t tbar +#tztz 0.30 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter Z t -> Z t +#ztzt 0.30 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter Z tbar -> Z tbar +#ttzz_a 0.30 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter t tbar -> Z Z +#zztt_a 0.30 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 1 scatter Z Z -> t tbar +#tztz_a 0.30 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter Z t -> Z t +#ztzt_a 0.30 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter Z tbar -> Z tbar +# Ward identities for W bosons +#wwtt 0.35 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter t tbar -> W+ W- +#wtwt 0.35 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter W+ t -> W+ t +#wwtt_a 0.35 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter t tbar -> W+ W- +#wtwt_a 0.35 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter W+ t -> W+ t +# Ward identities for W bosons +#wwbb 0.35 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter b bbar -> W+ W- +#wbwb 0.35 100 1000 0.000 0.000 0.000 0.000 SM 3 scatter W+ b -> W+ b +#wwbb_a 0.35 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter b bbar -> W+ W- +#wbwb_a 0.35 100 1000 0.000 0.000 0.000 0.000 SM_top_anom 3 scatter W+ b -> W+ b +# +wwaz 0.15 100 1000 80.419 80.419 0.000 91.1882 SM 3 scatter W+ W- -> A Z +wwza 0.15 100 1000 80.419 80.419 91.1882 0.000 SM 4 scatter W+ W- -> Z A +azww 0.15 100 1000 0.000 91.1882 80.419 80.419 SM 1 scatter A Z -> W+ W- +wawz 0.15 100 1000 80.419 0.000 80.419 91.1882 SM 2 scatter W+ A -> W+ Z +wwaa 0.20 100 1000 80.419 80.419 0.000 0.000 SM 3 scatter W+ W- -> A A +aaww 0.20 100 1000 0.000 0.000 80.419 80.419 SM 1 scatter A A -> W+ W- +awaw 0.20 100 1000 0.000 80.419 0.000 80.419 SM 1 scatter A W+ -> A W+ Index: trunk/omega/tests/parameters_SM_top_anom.f90 =================================================================== --- trunk/omega/tests/parameters_SM_top_anom.f90 (revision 8756) +++ trunk/omega/tests/parameters_SM_top_anom.f90 (revision 8757) @@ -1,530 +1,548 @@ ! parameters.SM_top_anom.f90 -- ! ! 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. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module parameters_sm_top_anom use kinds use constants implicit none private real(default), dimension(27), public :: mass, width real(default), public :: as complex(default), public :: gs, igs, ig, unit, half real(default), public :: e, g, e_em real(default), public :: sinthw, costhw, sin2thw, tanthw real(default), public :: qelep, qeup, qedwn real(default), public :: ttop, tbot, tch, ttau, tw real(default), public :: ltop, lbot, lc, ltau, lw complex(default), public :: qlep, qup, qdwn, gcc, qw, & gzww, gwww, ghww, ghhww, ghzz, ghhzz, & ghbb, ghtt, ghcc, ghtautau, gh3, gh4, ghmm, & iqw, igzww, igwww, gw4, gzzww, gazww, gaaww, & gvl_qbub_n, gvl_qw, gvl_qw_u, gvl_qw_d, & gsl_dttr, gsr_dttr, gsl_dttl, gsl_dbtl, & c_quqd1_1, c_quqd1_2, c_quqd8_1, c_quqd8_2 real(default), public :: vev, lambda, gi_flag, norm_flag, norm, & n_tvaa, n_vlrz, n_tvaz, n_vlrw, n_tlrw, n_tvag, n_sph complex(default), dimension(2), public :: & gncneu, gnclep, gncup, gncdwn, & - tvaa, tvaabb, vlrz, vlrcz, tvaz, tcvaz, tvazbb, tcvaa, tuvaa, & - vlrw, tlrw, tvag, tcvag, tuvag, sph, & + tvaa, tvaabb, vlrz, vlrcz, vlruz, tvaz, tcvaz, tvazbb, tcvaa, tuvaa, & + tuvaz, vlrw, tlrw, tvag, tcvag, tuvag, sph, & gvlr_qbub, gvlr_qbub_u, gvlr_qbub_d, gvlr_qbub_e, & gvlr_qgug, gslr_dbtr integer, public :: fun_flag logical, public :: bz=.false., bw=.false., ba=.false. public :: init_parameters, model_update_alpha_s, & - gmom, gtva_tta, gtva_tca, gtva_tua, & + gmom, gtva_tta, gtva_tca, gtva_tua, gtva_tuz, gvlr_tuz, & gvlr_ttz, gvlr_tcz, gtva_ttz, gtva_tcz, gvlr_btw, gvlr_tbw, & gtlr_btw, gtrl_tbw, gtlr_btwz, gtrl_tbwz, gtlr_btwa, gtrl_tbwa, & gtva_ttww, gtva_bba, gtva_bbz, gtva_bbww, & gtva_ttg, gtva_ttgg, gtva_tcg, gtva_tug, gtva_tcgg, gtva_tugg, gsp_tth real(default), parameter :: & GF = 1.16639E-5_default ! Fermi constant !!! This corresponds to 1/alpha = 137.03598949333 real(default), parameter :: & alpha = 1.0_default/137.03598949333_default complex(default), parameter :: & alphas = 0.1178_default ! Strong coupling constant (Z point) contains subroutine init_parameters tvaa(1) = 1 tvaa(2) = 1 * imago tcvaa(1) = 1 tcvaa(2) = 1 * imago tuvaa(1) = 1 tuvaa(2) = 1 * imago + tuvaz(1) = 1 + tuvaz(2) = 1 * imago vlrz(1) = 1 vlrz(2) = 1 vlrcz(1) = 1 vlrcz(2) = 1 * imago + vlruz(1) = 1 + vlruz(2) = 1 * imago tvaz(1) = 1 tvaz(2) = 1 * imago tcvaz(1) = 1 tcvaz(2) = 1 * imago vlrw(1) = 1 + 1 * imago vlrw(2) = 1 + 1 * imago tlrw(1) = 1 + 1 * imago tlrw(2) = 1 + 1 * imago tvag(1) = 1 tvag(2) = 1 * imago tcvag(1) = 1 tcvag(2) = 1 * imago tuvag(1) = 1 tuvag(2) = 1 * imago sph(1) = 0 sph(2) = 0 lambda = 2000 fun_flag = 0 ! norm_flag = 1 gi_flag = 1 mass(1:27) = 0 width(1:27) = 0 ! mass(3) = 0.095_default ! s-quark mass ! mass(4) = 1.2_default ! c-quark mass -! mass(5) = 4.2_default ! b-quark mass -! mass(6) = 173.1_default ! t-quark mass + mass(5) = 4.2_default ! b-quark mass + mass(6) = 173.1_default ! t-quark mass ! width(6) = 1.523_default ! t-quark width ! mass(11) = 0.000510997_default ! electron mass ! mass(13) = 0.105658389_default ! muon mass ! mass(15) = 1.77705_default ! tau-lepton mas mass(23) = 91.1882_default ! Z-boson mass width(23) = 2.443_default ! Z-boson width mass(24) = 80.419_default ! W-boson mass width(24) = 2.049_default ! W-boson width mass(25) = 200._default ! Higgs mass width(25) = 1.419_default ! Higgs width ttop = 4.0_default * mass(6)**2 / mass(25)**2 tbot = 4.0_default * mass(5)**2 / mass(25)**2 tch = 4.0_default * mass(4)**2 / mass(25)**2 ttau = 4.0_default * mass(15)**2 / mass(25)**2 tw = 4.0_default * mass(24)**2 / mass(25)**2 ! ltop = 4.0_default * mass(6)**2 / mass(23)**2 ! lbot = 4.0_default * mass(5)**2 / mass(23)**2 ! lc = 4.0_default * mass(4)**2 / mass(23)**2 ! ltau = 4.0_default * mass(15)**2 / mass(23)**2 ! lw = 4.0_default * mass(24)**2 / mass(23)**2 e_em = sqrt(4.0_default * PI * alpha) vev = 1 / sqrt (sqrt (2.0_default) * GF) ! v (Higgs vev) ! costhw = mass(24) / mass(23) ! cos(theta-W) costhw = 0.881901_default sinthw = sqrt (1.0_default-costhw**2) ! sin(theta-W) sin2thw = sinthw**2 tanthw = sinthw/costhw ! e = 2.0_default * sinthw * mass(24) / vev ! em-coupling (GF scheme) e = 0.349196_default qelep = - 1 qeup = 2.0_default / 3.0_default qedwn = - 1.0_default / 3.0_default g = e / sinthw ig = cmplx (0.0_default, 1.0_default, kind=default) * g unit = 1.0_default half = 0.5_default gcc = - g / 2 / sqrt (2.0_default) gncneu(1) = - g / 2 / costhw * ( + 0.5_default) gnclep(1) = - g / 2 / costhw * ( - 0.5_default - 2 * qelep * sin2thw) gncup(1) = - g / 2 / costhw * ( + 0.5_default - 2 * qeup * sin2thw) gncdwn(1) = - g / 2 / costhw * ( - 0.5_default - 2 * qedwn * sin2thw) gncneu(2) = - g / 2 / costhw * ( + 0.5_default) gnclep(2) = - g / 2 / costhw * ( - 0.5_default) gncup(2) = - g / 2 / costhw * ( + 0.5_default) gncdwn(2) = - g / 2 / costhw * ( - 0.5_default) qlep = - e * qelep qup = - e * qeup qdwn = - e * qedwn qw = e iqw = imago*qw gzww = g * costhw igzww = imago*gzww gwww = g igwww = imago*gwww gw4 = gwww**2 gzzww = gzww**2 gazww = gzww * qw gaaww = qw**2 ghww = mass(24) * g ghhww = g**2 / 2.0_default ghzz = mass(23) * g / costhw ghhzz = g**2 / 2.0_default / costhw**2 ghtt = - mass(6) / vev ghbb = - mass(5) / vev ghcc = - mass(4) / vev ghtautau = - mass(15) / vev ghmm = - mass(13) / vev gh3 = - 3 * mass(25)**2 / vev gh4 = - 3 * mass(25)**2 / vev**2 !!! Color flow basis, divide by sqrt(2) gs = sqrt(2.0_default*PI*alphas) igs = cmplx (0.0_default, 1.0_default, kind=default) * gs norm = vev / lambda**2 n_tvaa = e * norm n_vlrz = mass(23) * norm / 2 n_tvaz = norm n_vlrw = sqrt(2.0_default) * mass(24) * norm / 2 n_tlrw = sqrt(2.0_default) * norm / 2 ! Assuming that imag(gs).eq.0 ... n_tvag = real (gs, kind=default) * norm n_sph = sqrt(0.5_default) * vev * norm tvaabb(1) = 0.0_default tvaabb(2) = 0.0_default tvazbb(1) = 0.0_default tvazbb(2) = 0.0_default if ( gi_flag > 0. ) then if ( abs(real(vlrw(1))) > 0. ) then print *, "WARNING: gauge invariance and vanishing anomalous" print *, " bbZ vector coupling implies a relation" print *, " vl_ttZ ~ vl_tbW_Re ." print *, " Inferring vl_ttZ from vl_tbW_Re and IGNORING any" print *, " inconsistent values set!" vlrz(1) = real(vlrw(1)) * sqrt(2.0_default) / costhw print *, " vl_ttZ = ", real(vlrz(1))/n_vlrz end if if ( ( abs(real(tlrw(1))) > 0. ).or.( abs(aimag(tlrw(1))) > 0. ) ) then print *, "WARNING: anomalous tbW tensor couplings are related to anomalous" print *, " bbZ and bbA tensor couplings by gauge invariance:" print *, " Inferring bottom couplings from tl_tbW" tvaabb(1) = real(tlrw(1)) * sinthw / sqrt(2.0_default) tvaabb(2) = aimag(tlrw(1))*imago * sinthw / sqrt(2.0_default) tvazbb(1) = real(tlrw(1)) * costhw / sqrt(2.0_default) tvazbb(2) = aimag(tlrw(1))*imago * costhw / sqrt(2.0_default) print *, " tv_bbA = ", real(tvaabb(1))/n_tvaa print *, " ta_bbA = ", aimag(tvaabb(2))/n_tvaa print *, " tv_bbZ = ", real(tvazbb(1))/n_tvaz print *, " ta_bbZ = ", aimag(tvazbb(2))/n_tvaz end if if ( ( abs(real(tvaz(1))) > 0. ).or.( abs(aimag(tvaz(2))) > 0. ) ) then bz = .true. end if if ( ( abs(real(tlrw(2))) > 0. ).or.( abs(real(tlrw(2))) > 0. ) ) then bw = .true. end if if ( ( abs(real(tvaa(1))) > 0. ).or.( abs(aimag(tvaa(2)))> 0. ) ) then ba = .true. end if if ( bz.or.bw.or.ba ) then print *, "WARNING: top anomalous tensor couplings to W, A and Z" print *, " are interrelated by gauge invariance:" print *, " Inferring Z couplings from W/A couplings according to" print *, " the relation in the model file and IGNORING any inconsistent" print *, " values set manually! (Exception: only tX_ttZ != 0:" print *, " tr_tbW ~ tv_ttZ + i*ta_ttZ and tX_ttA = 0)" if ( ( bw.and.bz ).and..not.ba ) then tvaa(1) = ( real(tlrw(2)) - sqrt(2.0_default)*costhw*tvaz(1) ) / ( 2.0_default*sinthw ) tvaa(2) = ( aimag(tlrw(2))*imago - sqrt(2.0_default)*costhw*tvaz(2) ) / ( 2.0_default*sinthw ) else if ( bz.and..not.bw ) then tlrw(2) = sqrt(2.0_default)*costhw*( tvaz(1) + tvaz(2) ) + 2.0_default*sinthw*( tvaa(1) + tvaa(2) ) else tvaz(1) = ( real(tlrw(2)) - 2.0_default*sinthw*tvaa(1) ) / ( sqrt(2.0_default)*costhw ) tvaz(2) = ( aimag(tlrw(2))*imago - 2.0_default*sinthw*tvaa(2) ) / ( sqrt(2.0_default)*costhw ) end if print *, " tv_ttA = ", real(tvaa(1))/n_tvaa print *, " ta_ttA = ", aimag(tvaa(2))/n_tvaa print *, " tv_ttZ = ", real(tvaz(1))/n_tvaz print *, " ta_ttZ = ", aimag(tvaz(2))/n_tvaz print *, " tr_tbW = ", real(tlrw(2))/n_tlrw, " + ", aimag(tlrw(2))/n_tlrw, "I" end if end if gvlr_qgug(1) = 0 gvlr_qgug(2) = 0 gvlr_qbub(1) = 0 gvlr_qbub(2) = 0 gvlr_qbub_u(1) = 1.0_default / 3.0_default / 2 gvlr_qbub_u(2) = 4.0_default / 3.0_default / 2 gvlr_qbub_d(1) = 1.0_default / 3.0_default / 2 gvlr_qbub_d(2) = - 2.0_default / 3.0_default / 2 gvlr_qbub_e(1) = - 1.0_default / 2 gvlr_qbub_e(2) = - 2.0_default / 2 gvl_qbub_n = - 1.0_default / 2 gvl_qw = 0 gvl_qw_u = 0.5_default / 2 gvl_qw_d = - 0.5_default / 2 gsl_dttr = 0 gsr_dttr = 0 gsl_dttl = 0 gslr_dbtr(1) = 0 gslr_dbtr(2) = 0 gsl_dbtl = 0 c_quqd1_1 = 0 c_quqd1_2 = 0 c_quqd8_1 = 0 c_quqd8_2 = 0 end subroutine init_parameters subroutine model_update_alpha_s (alpha_s) real(default), intent(in) :: alpha_s gs = sqrt(2.0_default*PI*alpha_s) igs = cmplx (0.0_default, 1.0_default, kind=default) * gs end subroutine model_update_alpha_s pure function gmom (k2, i, coeff, lam) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i complex(default), dimension(2), intent(in) :: coeff real(default), intent(in) :: lam select case (fun_flag) case (0) c = coeff(i) case (1) c = coeff(i) / (1.0_default + k2/lam**2) case (2) c = coeff(i) / (1.0_default + k2/lam**2)**2 case (3) c = coeff(i) * exp(-k2/lam**2) case default c = 0 ! hoping we never get here ... end select end function gmom pure function gtva_tta (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tvaa, lambda) end function gtva_tta pure function gtva_tca (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tcvaa, lambda) end function gtva_tca pure function gtva_tua (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tuvaa, lambda) end function gtva_tua pure function gtva_bba (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tvaabb, lambda) end function gtva_bba pure function gvlr_ttz (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, vlrz, lambda) end function gvlr_ttz pure function gvlr_tcz (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, vlrcz, lambda) end function gvlr_tcz + pure function gvlr_tuz (k2, i) result (c) + complex(default) :: c + real(default), intent(in) :: k2 + integer, intent(in) :: i + c = - gmom (k2, i, vlruz, lambda) + end function gvlr_tuz + pure function gtva_ttz (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tvaz, lambda) end function gtva_ttz pure function gtva_tcz (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tcvaz, lambda) end function gtva_tcz + pure function gtva_tuz (k2, i) result (c) + complex(default) :: c + real(default), intent(in) :: k2 + integer, intent(in) :: i + c = - gmom (k2, i, tuvaz, lambda) + end function gtva_tuz + pure function gtva_bbz (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tvazbb, lambda) end function gtva_bbz pure function gvlr_btw (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, vlrw, lambda) end function gvlr_btw pure function gvlr_tbw (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, conjg(vlrw), lambda) end function gvlr_tbw pure function gtlr_btw (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tlrw, lambda) end function gtlr_btw pure function gtrl_tbw (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, conjg(tlrw), lambda) end function gtrl_tbw pure function gtlr_btwa (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i !!! don't touch this relative factor: fixed by ward identity! c = - sinthw * gtlr_btw(k2, i) end function gtlr_btwa pure function gtrl_tbwa (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i !!! don't touch this relative factor: fixed by ward identity! c = - sinthw * gtrl_tbw(k2, i) end function gtrl_tbwa pure function gtlr_btwz (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i !!! don't touch this relative factor: fixed by ward identity! c = - costhw * gtlr_btw(k2, i) end function gtlr_btwz pure function gtrl_tbwz (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i !!! don't touch this relative factor: fixed by ward identity! c = - costhw * gtrl_tbw(k2, i) end function gtrl_tbwz pure function gtva_ttww (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i !!! additional factor 2 wrt. ward relation to cancel explicit 1/2 in "gtlr" select case (i) case (1) c = - sqrt(2.0_default) * real(gtlr_btw(k2, 2)) case (2) c = - sqrt(2.0_default) * aimag(gtlr_btw(k2, 2)) * imago case default c = 0 ! hoping we never get here ... end select end function gtva_ttww pure function gtva_bbww (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i !!! additional factor 2 wrt. ward relation to cancel explicit 1/2 in "gtlr" select case (i) case (1) c = - sqrt(2.0_default) * real(gtlr_btw(k2, 1)) case (2) c = sqrt(2.0_default) * aimag(gtlr_btw(k2, 1)) * imago case default c = 0 ! hoping we never get here ... end select end function gtva_bbww pure function gtva_ttg (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tvag, lambda) end function gtva_ttg pure function gtva_ttgg (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i !!! don't touch this relative factor: fixed by ward identity! c = - gtva_ttg(k2, i) end function gtva_ttgg pure function gtva_tcg (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tcvag, lambda) end function gtva_tcg pure function gtva_tug (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, tuvag, lambda) end function gtva_tug pure function gtva_tcgg (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i !!! don't touch this relative factor: fixed by ward identity! c = - gtva_tcg(k2, i) end function gtva_tcgg pure function gtva_tugg (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i !!! don't touch this relative factor: fixed by ward identity! c = - gtva_tug(k2, i) end function gtva_tugg pure function gsp_tth (k2, i) result (c) complex(default) :: c real(default), intent(in) :: k2 integer, intent(in) :: i c = - gmom (k2, i, sph, lambda) end function gsp_tth end module parameters_sm_top_anom Index: trunk/omega/tests/ward_driver.sh =================================================================== --- trunk/omega/tests/ward_driver.sh (revision 8756) +++ trunk/omega/tests/ward_driver.sh (revision 8757) @@ -1,160 +1,175 @@ #! /bin/sh # ward_driver.sh -- ######################################################################## omega="$1" shift models="qed qcd sym sm sm_top_anom" modules="" ######################################################################## -while read module threshold n roots model unphysical mode process; do +# m1 m2 are the masses of the incoming particles +# m3 m4 are the masses of the first two outgoing particles, all further +# outgoing particles are assumed to be massless +######################################################################## +while read module threshold n roots m1 m2 m3 m4 model unphysical mode process; do case $module in '#'*) # skip comments ;; '') # skip empty lines ;; '!'*) break ;; *) ######################################################################## modules="$modules $module" eval threshold_$module=$threshold eval n_$module=$n eval roots_$module=$roots + eval m1_$module=$m1 + eval m2_$module=$m2 + eval m3_$module=$m3 + eval m4_$module=$m4 eval process_$module="'$process'" ######################################################################## ######################################################################## omega_bin="`echo $omega | sed s/%%%/$model/g`" # echo "running $omega_bin -$mode '$process'" 1>&2 $omega_bin "$@" \ -target:parameter_module parameters_$model \ -target:module amplitude_ward_physical_$module \ -$mode "$process" 2>/dev/null $omega_bin "$@" \ -target:parameter_module parameters_$model \ -target:module amplitude_ward_unphysical_$module \ -$mode "$process" -unphysical $unphysical 2>/dev/null ;; esac done ######################################################################## for module in $modules; do for mode in physical unphysical; do cat < number_particles_in p%number_particles_out => number_particles_out p%number_spin_states => number_spin_states p%spin_states => spin_states p%number_flavor_states => number_flavor_states p%flavor_states => flavor_states p%number_color_indices => number_color_indices p%number_color_flows => number_color_flows p%color_flows => color_flows p%number_color_factors => number_color_factors p%color_factors => color_factors p%color_sum => color_sum p%new_event => new_event p%reset_helicity_selection => reset_helicity_selection p%is_allowed => is_allowed p%get_amplitude => get_amplitude end function load end module interface_ward_${mode}_${module} EOF done done ######################################################################## cat < load EOF done done for model in $models; do cat < init_parameters EOF done cat < 0) then print *, failures, " failures in ", attempts, " attempts" failed_processes = failed_processes + 1 end if EOF done cat < 0) then print *, failed_processes, " failed processes in ", attempted_processes, " attempts" stop 1 end if end program ward_driver EOF exit 0