Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8723567
write_module.f
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
53 KB
Subscribers
None
write_module.f
View Options
!-----------------------
! madgraph - a Feynman Diagram package by Tim Stelzer and Bill Long
! (c) 1993
!
! Filename: write_module.f
!-----------------------
!***************************************************************************
!***************************************************************************
! The Fortran 90 interface for Madgraph to run with WHIZARD
! adapted by W. Kilian 2000-2004
!***************************************************************************
!***************************************************************************
subroutine
writeamp_f90
&
(
ngraphs
,
nexternal
,
nproc
,
sflows
,
&
unit_par_def
,
unit_par_set
)
!***************************************************************************
! writes out the subroutine for Fortran 90
!***************************************************************************
implicit none
! Constants
include
'params.inc'
! integer maxlines, maxfactors, maxflows
! parameter (maxlines=9,maxfactors=10,maxflows=200)
! integer maxgraphs
! parameter (maxgraphs=9999)
! Arguments
integer
ngraphs
,
nexternal
,
nproc
,
ncolor
double precision
sflows
(
maxflows
,
maxflows
)
integer
::
unit_par_def
,
unit_par_set
! Local Variables
character
(
120
)
buff
character
(
80
)
name
character
(
3
)
temp
character
(
3
)
ifmt1
,
ifmt2
,
ifmt3
character
(
80
)
fmt1
,
fmt2
character
*
4
str1
integer
i
,
j
,
k
,
ichar
,
jchar
,
ishift
,
nlength
,
nlen
integer
isym
,
jhel
,
initc
,
iden
(
maxcross
)
integer
ncross
,
ipart
(
maxlines
,
maxcross
),
ic
(
maxlines
,
maxcross
)
integer
jline
(
0
:
4
),
imatch
(
0
:
maxcross
)
integer
icross
,
iboson
integer
nhel
(
maxlines
,
3
**
maxlines
),
step
(
maxlines
),
ncomb
,
ncomb_in
integer
cflow
(
maxlines
,
maxflows
)
integer
aflow
(
maxlines
,
maxflows
)
integer
flow_c
(
maxflows
,
maxflows
,
2
)
integer
ncolors
,
nflows
integer
idenom
,
inum
,
icount
,
jcount
! External
logical
equal
external
equal
integer
igcf
external
igcf
! Global Variables
double precision
den
(
maxgraphs
)
integer
neigen
,
nrows
,
jplace
(
maxgraphs
)
common
/
to_colormat
/
den
,
neigen
,
nrows
,
jplace
integer
flows
(
0
:
2
*
maxlines
,
0
:
maxfactors
,
0
:
24
,
0
:
maxflows
)
integer
graphcolor
(
0
:
2
,
0
:
maxflows
,
0
:
maxgraphs
)
common
/
to_color
/
graphcolor
,
flows
integer
alen
(
maxgraphs
)
integer
wnum
(
-
maxlines
:
maxlines
),
wcount
,
wlen
(
maxgraphs
*
maxlines
)
common
/
toopt
/
wcount
,
wnum
,
wlen
,
alen
integer
iline
(
-
maxlines
:
maxlines
),
idir
(
-
maxlines
:
maxlines
)
integer
this_coup
(
max_coup
)
,
goal_coup
(
max_coup
)
common
/
to_proc
/
iline
,
idir
,
this_coup
,
goal_coup
character
*
(
10
)
coup_name
(
max_coup
)
integer
ncoups
common
/
to_couplings
/
ncoups
,
coup_name
character
*
80
gname
integer
iname
common
/
to_name
/
iname
,
gname
integer
nincoming
common
/
to_proc2
/
nincoming
character
*
60
proc
integer
iproc
common
/
to_write
/
iproc
,
proc
integer
symfact
(
maxgraphs
)
common
/
tosymmetry
/
symfact
integer
info_p
(
5
,
max_particles
),
iposx
(
3
,
max_particles
)
character
*
(
max_string
)
iwave
(
max_particles
),
owave
(
max_particles
)
character
*
(
8
)
str
(
3
,
max_particles
)
common
/
to_external
/
iwave
,
owave
,
iposx
,
info_p
,
str
character
*
(
4
*
max_particles
)
particle
(
4
)
integer
charge_c
(
max_particles
)
integer
iparticle
(
0
:
max_particles
,
0
:
4
),
inverse
(
max_particles
)
common
/
to_model
/
iparticle
,
particle
,
inverse
,
charge_c
logical
lwp
,
lwm
,
lz
,
decay
,
cross
common
/
to_decay
/
lwp
,
lwm
,
lz
,
decay
,
cross
!-----------
! Begin Code
!-----------
ishift
=
iname
nlength
=
ishift
open
(
unit
=
92
,
file
=
gname
(
1
:
iname
)
//
'.f90'
,
status
=
'unknown'
)
name
=
gname
(
1
:
iname
)
call
lower_case
(
name
)
print
*
,
'Writing function '
,
name
(
1
:
iname
)
,
&
' in file '
,
gname
(
1
:
iname
)
//
'.f90'
ncolors
=
flows
(
0
,
0
,
0
,
0
)
! Header
write
(
92
,
5
)
'WHIZARD interface for Madgraph'
write
(
92
,
5
)
write
(
92
,
5
)
'Automatically generated file, do not edit'
write
(
92
,
5
)
! Driver module
write
(
92
,
5
)
'Driver module'
write
(
92
,
5
)
'Process: '
//
trim
(
proc
)
do
i
=
1
,
ncoups
write
(
92
,
fmt
=
5
,
advance
=
"no"
)
' PT order: '
//
&
trim
(
coup_name
(
i
))
//
' <= '
write
(
92
,
*
)
goal_coup
(
i
)
end do
write
(
92
,
5
)
write
(
92
,
10
)
'module '
//
trim
(
name
)
write
(
92
,
10
)
write
(
92
,
10
)
' use kinds, only: D=>default !NODEP!'
write
(
92
,
10
)
' use parameters'
write
(
92
,
10
)
' implicit none'
write
(
92
,
10
)
' private'
write
(
92
,
10
)
write
(
92
,
10
)
" public :: scatter_nonzero"
write
(
92
,
10
)
" public :: scatter_colored_nonzero"
write
(
92
,
10
)
" public :: scatter_diagonal_nonzero"
write
(
92
,
10
)
" public :: scatter_diagonal_colored_nz"
write
(
92
,
10
)
" public :: number_particles, "
//
&
"number_particles_in, number_particles_out"
write
(
92
,
10
)
" public :: number_spin_states, "
//
&
"number_spin_states_in, number_spin_states_out"
write
(
92
,
10
)
" public :: spin_states, "
//
&
"spin_states_in, spin_states_out"
write
(
92
,
10
)
" public :: number_flavor_states, "
//
&
"number_flavor_states_in, number_flavor_states_out"
write
(
92
,
10
)
" public :: flavor_states, "
//
&
"flavor_states_in, flavor_states_out"
write
(
92
,
10
)
" public :: number_flavor_zeros, "
//
&
"number_flavor_zeros_in, number_flavor_zeros_out"
write
(
92
,
10
)
" public :: flavor_zeros, "
//
&
"flavor_zeros_in, flavor_zeros_out"
write
(
92
,
10
)
" public :: number_color_flows"
write
(
92
,
10
)
" public :: color_flows, anticolor_flows"
write
(
92
,
10
)
" public :: create, reset, destroy"
write
(
92
,
10
)
! Code for decays, crossings (jets etc.) -- not activated yet
! c
! c These lines for artificial decay
! c
! if (decay) then
! call getparticle('w+',nwp,ngot)
! if (ngot .ne. 1) print*,'Error finding W+',nwp,ngot
! call getparticle('w-',nwm,ngot)
! if (ngot .ne. 1) print*,'Error finding W+',nwp,ngot
! call getparticle('z',nz,ngot)
! if (ngot .ne. 1) print*,'Error finding W+',nwp,ngot
! lwp=.false.
! lwm=.false.
! lz =.false.
! if (iline(iline(0)) .eq. nwp) lwp=.true.
! if (iline(iline(0)) .eq. nwm) lwm=.true.
! if (iline(iline(0)) .eq. nz ) lz =.true.
! if (lz) then
! call getparticle('e+e-',jline(1),jline(0))
! elseif(lwp) then
! call getparticle('e+ve',jline(1),jline(0))
! elseif(lwm) then
! call getparticle('ve~e-',jline(1),jline(0))
! else
! print*,'Warning last particle not Z or W+ or W-'
! print*,'Not decaying'
! decay=.false.
! endif
! endif
! if (decay) then
! iboson=iline(iline(0))
! iline(0)=iline(0)+1
! iline(iline(0)-1)=jline(1)
! iline(iline(0))=jline(2)
! c print*,'We have decay',decay,lz
! else
! c print*,'No decay'
! endif
!
if
(
cross
)
then
! call crossing(ncross,ipart,imatch)
stop
'Bug: crossing called in WHIZARD version'
else
ncross
=
1
imatch
(
0
)
=
1
imatch
(
1
)
=
1
do
i
=
1
,
iline
(
0
)
ipart
(
i
,
1
)
=
i
enddo
endif
c
c
set
up
all
the
different
crossings
which
are
not
identical
c
icross
=
0
do
i
=
1
,
ncross
if
(
imatch
(
i
)
.
gt
.
icross
)
then
call
getdenom
(
isym
,
initc
,
jhel
,
step
,
ipart
(
1
,
i
))
iden
(
imatch
(
i
))
=
isym
*
initc
*
jhel
c
write
(
*
,
*
)
'iden'
,
i
,
imatch
(
i
),
iden
(
imatch
(
i
))
if
(
imatch
(
i
)
.
ne
.
icross
+
1
)
then
print
*
,
'Error in imatch writesub.f'
,
i
,
imatch
(
i
),
icross
+
1
endif
icross
=
imatch
(
i
)
endif
enddo
if
(
icross
.
ne
.
imatch
(
0
))
then
print
*
,
'Warning icross .ne. imatch(0)'
,
icross
,
imatch
(
0
)
endif
c
write
(
*
,
*
)
'iden(1)'
,
iden
(
1
)
ishift
=
iname
c
nlength
=
ishift
! Number of particles
nexternal
=
iline
(
0
)
! Helicity states
ncomb_in
=
1
do
i
=
1
,
nexternal
if
(
info_p
(
2
,
iline
(
i
))
.
eq
.
1
)
then
step
(
i
)
=
3
!1 Helicity state 0 (scalar)
elseif
(
abs
(
info_p
(
2
,
iline
(
i
)))
.
eq
.
2
)
then
step
(
i
)
=
2
!2 Helicity states -1,1
if
(
i
.
le
.
nincoming
)
ncomb_in
=
ncomb_in
*
2
elseif
(
info_p
(
2
,
iline
(
i
))
.
eq
.
3
)
then
step
(
i
)
=
1
!3 Helicity state -1,0,1 W,Z
if
(
i
.
le
.
nincoming
)
ncomb_in
=
ncomb_in
*
3
else
print
*
,
'Unknown Helicity'
,
info_p
(
2
,
iline
(
i
)),
&
iline
(
i
),
i
stop
endif
enddo
call
sethel_f90
(
nhel
,
step
,
nexternal
,
ncomb
)
if
(
ncomb
/=
ncomb_in
*
(
ncomb
/
ncomb_in
))
then
stop
"Can't happen: mismatch in helicity counting"
end if
!!! Code for summing over partons (?) -- not activated yet
! c buff = 'REAL*8 FUNCTION S'
! c jchar = 18
! buff = 'SUBROUTINE S'
! jchar = 13
! buff(jchar:jchar+ishift) = name(1:iname)
! jchar = jchar+ishift
! ishift= 7
! buff(jchar:jchar+ishift) = '(P1,ANS)'
! jchar =jchar+ishift
! write(92,10) buff(1:jchar)
!
! call genparton(ncross,ipart,buff(12:jchar),nproc,imatch)
!
! jchar=jchar-4
! FNAME= buff(13:jchar)
! c FNAME(JCHAR-12:JCHAR-12+13)=',NHEL(1,IHEL))'
! FNAME(JCHAR-12:JCHAR-12+19)=',NHEL(1,IHEL),JC(1))'
! FLEN = jchar+19
!
!!! Code for crossings -- not activated yet
! do i=1,ncross
! buff = 'Crossing ??? is '
! write(buff(10:12),'(i3)') imatch(i)
! ichar = 17
! do j=1,2
! call part_string(inverse(iline(ipart(j,i))),str1,nlen)
! buff(ichar:ichar+nlen-1) = str1(1:nlen)
! ichar=ichar+nlen
! buff(ichar:ichar)=' '
! ichar=ichar+1
! enddo
! buff(ichar:ichar+3) = '-> '
! ichar=ichar+3
! do j=3,iline(0)
! call part_string(iline(ipart(j,i)),str1,nlen)
! buff(ichar:ichar+nlen-1) = str1(1:nlen)
! ichar=ichar+nlen
! buff(ichar:ichar)=' '
! ichar=ichar+1
! enddo
! write(92,5) buff(1:ichar)
! enddo
! Compute the color flow configurations and coefficients:
! Note that, in Madgraph2, the 'flow' amplitudes are actually color
! amplitudes (skeleton graphs of quarks and gluons),
! so we have to extract the color flow information by ourselves
! This routine is called here because we want to know the value of nflows
call
setflows
(
nflows
,
cflow
,
aflow
,
flow_c
)
! Now all constants are known
write
(
92
,
5
)
' '
write
(
92
,
5
)
'Constants '
write
(
92
,
5
)
' '
call
get_iformat
&
(
max
(
nexternal
,
ngraphs
,
wcount
,
&
nflows
,
ncolors
,
ncomb
,
icross
),
&
ifmt1
)
fmt1
=
"(2x,A,1x,"
//
ifmt1
//
")"
write
(
92
,
fmt1
)
'integer, parameter :: ncross ='
,
icross
write
(
92
,
fmt1
)
'integer, parameter :: nexternal ='
,
nexternal
write
(
92
,
fmt1
)
'integer, parameter :: nincoming ='
,
nincoming
write
(
92
,
fmt1
)
'integer, parameter :: noutgoing ='
,
&
nexternal
-
nincoming
write
(
92
,
fmt1
)
'integer, parameter :: ncomb ='
,
ncomb
write
(
92
,
fmt1
)
'integer, parameter :: ncomb_in ='
,
ncomb_in
write
(
92
,
fmt1
)
'integer, parameter :: ncomb_out ='
,
&
ncomb
/
ncomb_in
write
(
92
,
10
)
' integer, parameter :: thel = ncomb * ncross'
write
(
92
,
fmt1
)
'integer, parameter :: ncolors ='
,
ncolors
write
(
92
,
fmt1
)
'integer, parameter :: nflows ='
,
nflows
write
(
92
,
fmt1
)
'integer, parameter :: ngraphs ='
,
ngraphs
write
(
92
,
fmt1
)
'integer, parameter :: wcount ='
,
wcount
write
(
92
,
10
)
' integer, parameter :: SAMPLE = 10, REPEAT = 5'
write
(
92
,
5
)
' '
write
(
92
,
5
)
'Module Variables '
write
(
92
,
5
)
' '
write
(
92
,
10
)
' integer, dimension(nexternal,ncomb) :: nhel'
write
(
92
,
10
)
' integer, dimension(ncross) :: iden'
write
(
92
,
10
)
' real(D), dimension(ncross) :: cmult'
write
(
92
,
10
)
' integer, dimension(nexternal,ncross) :: ic'
write
(
92
,
10
)
' integer, dimension(nexternal,nflows) :: cflw'
write
(
92
,
10
)
' integer, dimension(nexternal,nflows) :: aflw'
write
(
92
,
10
)
' integer, dimension(ncolors,nflows) :: '
//
&
'flow_c_num'
write
(
92
,
10
)
' integer, dimension(ncolors,nflows) :: '
//
&
'flow_c_den'
write
(
92
,
10
)
' real(D), dimension(ncolors,nflows) :: flow_c'
write
(
92
,
10
)
' integer, dimension(ncolors,ncolors) :: cf_num'
write
(
92
,
10
)
' integer, dimension(ncolors) :: cf_den'
write
(
92
,
10
)
' real(D), dimension(ncolors,ncolors) :: cf'
write
(
92
,
10
)
' complex(D), dimension(ncolors,ncomb,ncross) ::'
&
//
' amp_color'
write
(
92
,
10
)
' complex(D), dimension(nflows,ncomb,ncross) ::'
&
//
' amp_flow'
c
Transfer
parameter
declarations
from
tmp
file
c
write
(
92
,
5
)
' '
write
(
92
,
5
)
'Physics parameters'
write
(
92
,
5
)
' '
rewind
(
unit_par_def
)
do
read
(
unit_par_def
,
'(A)'
,
end
=
70
)
buff
write
(
92
,
'(A)'
)
trim
(
buff
)
cycle
70
exit
end do
! write data tables
write
(
92
,
5
)
' '
write
(
92
,
5
)
'Helicity combination tables'
write
(
92
,
5
)
' '
call
get_iformat
(
ncomb
,
ifmt1
)
call
data_format
(
fmt1
,
2
,
",':,',"
//
ifmt1
,
nexternal
,
"I2"
)
do
i
=
1
,
ncomb
write
(
92
,
fmt1
)
"nhel"
,
i
,
nhel
(:
nexternal
,
i
)
end do
! Code for crossings
write
(
92
,
5
)
' '
write
(
92
,
5
)
'Crossing (permutation) data'
write
(
92
,
5
)
' '
call
get_iformat
(
nexternal
,
ifmt1
)
call
get_iformat
(
ncross
,
ifmt2
)
call
get_iformat
(
maxval
(
iden
(:
ncross
)),
ifmt3
)
call
data_format
(
fmt1
,
2
,
",':,',"
//
ifmt2
,
nexternal
,
ifmt1
)
call
data_format
(
fmt2
,
2
,
ifmt2
,
1
,
ifmt3
)
icross
=
0
do
i
=
1
,
ncross
if
(
imatch
(
i
)
.
gt
.
icross
)
then
write
(
92
,
fmt1
)
"ic"
,
i
,
ipart
(:
nexternal
,
i
)
write
(
92
,
fmt2
)
"iden"
,
i
,
iden
(
i
)
icross
=
icross
+
1
endif
enddo
! Write the color flow configurations, color ME, and coefficients
write
(
92
,
5
)
' '
write
(
92
,
5
)
'Color flow tables'
write
(
92
,
5
)
' '
call
get_iformat
(
nflows
,
ifmt1
)
call
data_format
(
fmt1
,
2
,
",':,',"
//
ifmt1
,
nexternal
,
"I2"
)
do
i
=
1
,
nflows
write
(
92
,
fmt1
)
"cflw"
,
i
,
cflow
(:
nexternal
,
i
)
end do
write
(
92
,
10
)
do
i
=
1
,
nflows
write
(
92
,
fmt1
)
"aflw"
,
i
,
aflow
(:
nexternal
,
i
)
end do
write
(
92
,
5
)
' '
write
(
92
,
5
)
'Coefficients for projecting color graphs onto '
//
&
'color flows'
write
(
92
,
5
)
' '
call
write_flow_coefficients_f90
(
ncolors
,
nflows
,
flow_c
)
write
(
92
,
5
)
' '
write
(
92
,
5
)
'Coefficients for squaring color graph amplitudes'
write
(
92
,
5
)
' '
call
write_color_matrix_f90
(
sflows
,
ncolors
)
write
(
92
,
5
)
' '
write
(
92
,
10
)
' save'
write
(
92
,
5
)
' '
write
(
92
,
10
)
'contains'
write
(
92
,
10
)
! User-level routines
call
write_interface
&
(
nexternal
,
nhel
,
unit_par_def
,
unit_par_set
)
!-----------------------------------------------------------------------------
! Write routines for computing eigenamplitudes (to be summed over eigenvalues)
! and optionally color flow amplitudes
!-----------------------------------------------------------------------------
! Subroutine: Loop over helicities, compute and store color amplitude values
write
(
92
,
5
)
write
(
92
,
5
)
'This subroutine calculates the amplitudes '
//
&
'for all crossings'
write
(
92
,
5
)
'and helicity combinations and stores them '
//
&
'in module variables '
write
(
92
,
5
)
' '
write
(
92
,
10
)
' subroutine calculate_amplitudes '
//
&
'(p, zero_ct, n, calc_flows)'
write
(
92
,
10
)
' '
//
&
'real(D), dimension(0:3,nexternal), intent(in) :: p'
write
(
92
,
10
)
' integer, dimension(:,:,:,:), intent(inout) '
//
&
':: zero_ct'
write
(
92
,
10
)
' integer, intent(in) :: n'
write
(
92
,
10
)
' logical, intent(in) :: calc_flows'
! write local declarations
write
(
92
,
10
)
' real(D), dimension(0:3,nexternal) :: pp'
write
(
92
,
10
)
' integer :: hi, ho, hel'
write
(
92
,
10
)
' integer :: icross, ipart'
write
(
92
,
10
)
' integer, dimension(nexternal) :: jc'
write
(
92
,
10
)
' do icross = 1, ncross'
! write(92,10) ' call switchmom (p,pp,ic(:,icross),jc,nexternal)
c
c
the
'switchmom'
routine
inlined
:
c
write
(
92
,
10
)
' do ipart = 1, nexternal'
write
(
92
,
10
)
' pp(:,ic(ipart,icross)) = p(:,ipart)'
write
(
92
,
10
)
' end do'
write
(
92
,
10
)
' jc = +1'
c
c
Careful
here
,
if
there
are
no
crossings
then
we
don
't want
c to flip the signs of the incoming particles because they already
c are flipped
c
if (nincoming .eq. 0) then
write(92,10) '
jc
(
ic
(
1
))
=
-
1
'
write(92,10) '
jc
(
ic
(
2
))
=
-
1
'
end if
write(92,10) '
hel
=
0
'
write(92,10) '
do
hi
=
1
,
number_spin_states_in
()
'
write(92,10) '
do
ho
=
1
,
number_spin_states_out
()
'
write(92,10) '
hel
=
hel
+
1
'
write(92,10) '
if
(
calc_flows
)
then
'
write(92,10) '
call
calculate_amplitude
&
'
write(92,10) '
'//
& '
&
(
pp
,
nhel
(:,
hel
),
jc
,
zero_ct
(
hi
,
1
,
ho
,
1
),
n
,
&
'
write(92,10) '
'//
& '
&
amp_color
(:,
hel
,
icross
),
amp_flow
(:,
hel
,
icross
))
'
write(92,10) '
else
'
write(92,10) '
call
calculate_amplitude
&
'
write(92,10) '
'//
& '
&
(
pp
,
nhel
(:,
hel
),
jc
,
zero_ct
(
hi
,
1
,
ho
,
1
),
n
,
&
'
write(92,10) '
'//
& '
&
amp_color
(:,
hel
,
icross
))
'
write(92,10) '
end if
'
write(92,10) '
end do
'
write(92,10) '
end do
'
write(92,10) '
end do
'
write(92,10) '
end subroutine
calculate_amplitudes
'
write(92,10)
write(92,10)
!------------------------------------------------------------------------
! Subroutine: Matrix element calculation, returning diagram amplitudes
! and transforming them into appropriate color bases
write(92,5)
write(92,5) '
This
subroutine
calculates
the
vector
'//
& '
of
color
amplitudes
'
write(92,5) '
for
a
given
helicity
combination
'
write(92,5) '
'
write(92,5) '
Process
'//proc(1:iproc)
write(92,5) '
'
write(92,10) '
subroutine
calculate_amplitude
'//
& '
(
p
,
nhel
,
ic
,
zero_ct
,
n
,
amp_color
,
amp_flow
)
'
write(92,5) '
'
write(92,10) '
'//
& '
real
(
D
),
dimension
(
0
:
3
,
nexternal
),
intent
(
in
)
::
p
'
write(92,10) '
'//
& '
integer
,
dimension
(
nexternal
),
intent
(
in
)
::
nhel
'
write(92,10) '
'//
& '
integer
,
dimension
(
nexternal
),
intent
(
in
)
::
ic
'
write(92,10) '
'//
& '
integer
,
intent
(
inout
)
::
zero_ct
'
write(92,10) '
integer
,
intent
(
in
)
::
n
'
write(92,10) '
'//
& '
complex
(
D
),
dimension
(:),
intent
(
out
)
'//
& '
::
amp_color
'
write(92,10) '
'//
& '
complex
(
D
),
dimension
(:),
intent
(
out
),
optional
'//
& '
::
amp_flow
'
! write local declarations
write(92,5) '
'
write(92,10) '
complex
(
D
),
dimension
(
ngraphs
)
::
amp
'
write(92,10) '
complex
(
D
),
dimension
(
6
,
wcount
)
::
w
'
write(92,5) '
'
! Begin executable code
write(92,10)'
if
(
zero_ct
>
REPEAT
)
then
'
write(92,10)'
amp_color
(:)
=
0
'
write(92,10)'
if
(
present
(
amp_flow
))
amp_flow
(:)
=
0
'
write(92,10)'
return
'
write(92,10)'
end if
'
! Write HELAS calls
do while (.true.)
read(91,'
(
a
)
',end=990) buff
call lower_case (buff)
if (buff(6:6)==" ") then
write (92,*)
write (92,'
(
a
)
',advance="no") '
' //
& trim (adjustl (buff))
else
buff(6:6) = " "
write (92,'
(
a
)
',advance="no") trim (adjustl (buff))
end if
end do
990 continue
write (92, *)
rewind(91)
! Write flow amplitudes
write (92,5)
call get_iformat (ncolors, ifmt1)
call get_iformat (ngraphs, ifmt3)
do i=1,ncolors
c
c First determine a common denominator for these amplitudes
c
idenom=1
inum = 999999
do j=1,ngraphs
do k=1,graphcolor(0,0,j) !Number of terms/flows
if (graphcolor(0,k,j) .eq. i) then
inum = min(inum,abs(graphcolor(1,k,j)))
idenom=abs(idenom*graphcolor(2,k,j)/
& igcf(idenom,graphcolor(2,k,j)))
end if
end do
end do
call get_iformat (inum, ifmt2)
if (idenom .ne. 1) inum=1
call write_amp_color_begin (i, inum, idenom, .false.)
icount = 5
jcount = 0
do j=1,ngraphs
do k=1,graphcolor(0,0,j) !Number of terms/flows
if (graphcolor(0,k,j) .eq. i) then
if (jcount == 40) then
call write_amp_color_end (inum, idenom)
icount = 5
jcount = 0
call write_amp_color_begin (i,inum,idenom,.true.)
end if
if (icount == 5) then
write (92, "(A)") " &"
write (92, "(12x)", advance="no")
icount = 0
jcount = jcount + 1
end if
den(j)=1d0
den(j)=dble(graphcolor(1,k,j))/graphcolor(2,k,j)*
& idenom*symfact(j)/inum
if (equal(den(j),1d0)) then
write (92, "(A,"//ifmt3//",A)", advance="no")
& '
+
amp
(
',j,'
)
'
icount = icount + 1
elseif (equal(den(j),-1d0)) then
write (92, "(A,"//ifmt3//",A)", advance="no")
& '
-
amp
(
',j,'
)
'
icount = icount + 1
elseif (den(j) .gt. 0d0) then
write (92, "(A,E10.3,A,"//ifmt3//",A)",
& advance="no")
& '
+
', den(j), '
*
amp
(
',j,'
)
'
icount = icount + 1
elseif (den(j) .lt. 0d0) then
write (92, "(A,E10.3,A,"//ifmt3//",A)",
& advance="no")
& '
-
', abs(den(j)), '
*
amp
(
',j,'
)
'
icount = icount + 1
end if
end if
end do
end do
call write_amp_color_end (inum, idenom)
enddo
write(92,10)
write(92,10)'
if
(
n
<=
SAMPLE
)
then
'
write(92,10)'
if
(
all
(
abs
(
amp_color
)
<=
tiny
(
1.0
_
D
)))
'//
& '
zero_ct
=
zero_ct
+
1
'
write(92,10)'
end if
'
write(92,10)
write(92,10)'
if
(
present
(
amp_flow
))
then
'
write(92,10)'
amp_flow
=
matmul
(
amp_color
,
flow_c
)
'
write(92,10)'
end if
'
write(92,10)
write(92,10)'
end subroutine
calculate_amplitude
'
write(92,10)
write(92, 10) '
end module
'//trim(name)
close(92)
5 format('
! ',a)
6
format
(
'! '
,
a
,
i5
)
10
format
(
a
)
11
format
(
a
,
i5
)
c
10
format
(
6
x
,
a
)
c
11
format
(
6
x
,
a
,
i5
)
15
format
(
'!'
,
6
x
,
a
)
20
format
(
',p'
,
i1
)
contains
subroutine
write_amp_color_begin
(
i
,
inum
,
idenom
,
add
)
integer
,
intent
(
in
)
::
i
,
inum
,
idenom
logical
,
intent
(
in
)
::
add
character
(
80
)
fmt1
,
str
str
=
""
if
(
add
)
then
fmt1
=
"(1x,A,"
//
ifmt1
//
",A),1x"
write
(
str
,
fmt1
)
"amp_color("
,
i
,
")"
end if
if
(
inum
.
eq
.
1
.
and
.
idenom
.
eq
.
1
)
then
fmt1
=
"(4x,A,"
//
ifmt1
//
",A)"
write
(
92
,
fmt1
,
advance
=
"no"
)
"amp_color("
,
i
,
") ="
&
//
trim
(
str
)
else if
(
inum
.
ne
.
1
)
then
fmt1
=
"(4x,A,"
//
ifmt1
//
",A,"
//
ifmt2
//
",A)"
write
(
92
,
fmt1
,
advance
=
"no"
)
"amp_color("
,
i
,
") = "
&
//
trim
(
str
),
inum
,
"*("
else if
(
idenom
.
ne
.
1
)
then
fmt1
=
"(4x,A,"
//
ifmt1
//
",A)"
write
(
92
,
fmt1
,
advance
=
"no"
)
"amp_color("
,
i
,
") = "
&
//
trim
(
str
)
//
"("
end if
end subroutine
subroutine
write_amp_color_end
(
inum
,
idenom
)
integer
,
intent
(
in
)
::
inum
,
idenom
if
(
idenom
.
ne
.
1
)
then
write
(
92
,
'(A)'
)
" &"
write
(
92
,
'(10x,A,I6)'
)
') / '
,
idenom
else if
(
inum
.
ne
.
1
)
then
write
(
92
,
'(A)'
)
" &"
write
(
92
,
'(10x,A)'
)
')'
else
write
(
92
,
*
)
end if
end subroutine
end
subroutine
sethel_f90
(
nhel
,
step
,
npart
,
ncomb
)
c
************************************************************************
c
Sets
up
all
possible
helicity
combinations
in
nhel
(
ipart
,
ncomb
)
c
given
the
number
of
particles
,
and
the
step
size
to
be
used
c
for
each
particle
c
The
first
state
is
always
-
1
,
the
next
is
0
or
1
depending
onstep
c
WK
:
c
************************************************************************
implicit none
! Constants
include
'params.inc'
c
integer
maxlines
c
parameter
(
maxlines
=
9
)
! Arguments
integer
nhel
(
maxlines
,
3
**
maxlines
),
step
(
maxlines
),
npart
,
ncomb
! Local
integer
ipart
,
iflip
,
i
logical
gotone
ncomb
=
1
do
ipart
=
1
,
npart
if
(
step
(
ipart
)
==
3
)
then
nhel
(
ipart
,
ncomb
)
=
0
else
nhel
(
ipart
,
ncomb
)
=
-
1
end if
enddo
gotone
=
.
true
.
do while
(
gotone
)
c
write
(
*
,
'(8I4)'
),
ncomb
,(
nhel
(
ipart
,
ncomb
),
ipart
=
1
,
npart
)
ncomb
=
ncomb
+
1
do
ipart
=
1
,
npart
nhel
(
ipart
,
ncomb
)
=
nhel
(
ipart
,
ncomb
-
1
)
enddo
gotone
=
.
false
.
iflip
=
npart
do while
(.
not
.
gotone
.
and
.
iflip
.
gt
.
0
)
nhel
(
iflip
,
ncomb
)
=
nhel
(
iflip
,
ncomb
)
+
step
(
iflip
)
if
(
nhel
(
iflip
,
ncomb
)
.
gt
.
1
)
then
if
(
step
(
iflip
)
==
3
)
then
nhel
(
iflip
,
ncomb
)
=
0
else
nhel
(
iflip
,
ncomb
)
=
-
1
end if
iflip
=
iflip
-
1
else
gotone
=
.
true
.
endif
enddo
enddo
ncomb
=
ncomb
-
1
end
Subroutine
write_flow_coefficients_f90
(
ncolors
,
nflows
,
flow_c
)
!***************************************************************************
! This routine writes out the coefficients to project color graphs onto
! color flows
!***************************************************************************
implicit none
! Constants
include
'params.inc'
integer
npl
parameter
(
npl
=
6
)
! Arguments
integer
ncolors
,
nflows
integer
flow_c
(
maxflows
,
maxflows
,
2
)
! Local Variables
integer
i
,
j
,
k
character
(
3
)
ifmt1
,
ifmt2
,
ifmt3
character
(
80
)
fmt1
,
fmt2
!-----------
! Begin Code
!-----------
call
get_iformat
(
ncolors
,
ifmt1
)
call
get_iformat
(
nflows
,
ifmt2
)
call
get_iformat
(
maxval
(
flow_c
)
*
10
,
ifmt3
)
do
i
=
1
,
nflows
do
j
=
0
,
(
ncolors
-
1
)
/
npl
k
=
min
(
npl
,
ncolors
-
j
*
npl
)
call
data_format
(
fmt2
,
2
,
&
ifmt1
//
",':',"
//
ifmt1
//
",',',"
//
ifmt2
,
k
,
ifmt3
)
write
(
92
,
fmt2
)
&
"flow_c_num"
,
j
*
npl
+
1
,
j
*
npl
+
k
,
i
,
&
flow_c
(
j
*
npl
+
1
:
j
*
npl
+
k
,
i
,
1
)
write
(
92
,
fmt2
)
&
"flow_c_den"
,
j
*
npl
+
1
,
j
*
npl
+
k
,
i
,
&
flow_c
(
j
*
npl
+
1
:
j
*
npl
+
k
,
i
,
2
)
end do
end do
end
Subroutine
write_color_matrix_f90
(
sflows
,
ncolors
)
!***************************************************************************
! This routine writes out the matrix of color coefficients
! for all of the graphs.
!***************************************************************************
implicit none
! Constants
include
'params.inc'
integer
npl
parameter
(
npl
=
6
)
! Arguments
double precision
sflows
(
maxflows
,
maxflows
)
integer
ncolors
! Local Variables
integer
i
,
j
,
k
character
(
3
)
ifmt1
,
ifmt2
,
ifmt3
character
(
80
)
fmt1
,
fmt2
integer
idenom
(
ncolors
)
integer
icf
(
ncolors
,
ncolors
)
! External
integer
igcf
external
igcf
! Global Variables
integer
isflows
(
2
,
maxflows
,
maxflows
)
common
/
to_isflow
/
isflows
integer
flows
(
0
:
2
*
maxlines
,
0
:
maxfactors
,
0
:
24
,
0
:
maxflows
)
integer
graphcolor
(
0
:
2
,
0
:
maxflows
,
0
:
maxgraphs
)
common
/
to_color
/
graphcolor
,
flows
!-----------
! Begin Code
!-----------
do
i
=
1
,
ncolors
idenom
(
i
)
=
1
do
j
=
1
,
ncolors
if
(
isflows
(
2
,
j
,
i
)
/=
0
)
then
idenom
(
i
)
=
abs
(
idenom
(
i
)
*
isflows
(
2
,
j
,
i
)
/
&
igcf
(
idenom
(
i
),
isflows
(
2
,
j
,
i
)))
endif
enddo
do
j
=
1
,
ncolors
if
(
isflows
(
2
,
j
,
i
)
/=
0
)
then
icf
(
j
,
i
)
=
isflows
(
1
,
j
,
i
)
*
(
idenom
(
i
)
/
isflows
(
2
,
j
,
i
))
else
icf
(
j
,
i
)
=
0
end if
end do
end do
call
get_iformat
(
ncolors
,
ifmt1
)
call
get_iformat
(
maxval
(
idenom
),
ifmt2
)
call
get_iformat
(
maxval
(
icf
)
*
10
,
ifmt3
)
call
data_format
(
fmt1
,
2
,
ifmt1
,
1
,
ifmt2
)
do
i
=
1
,
ncolors
do
j
=
0
,
(
ncolors
-
1
)
/
npl
k
=
min
(
npl
,
ncolors
-
j
*
npl
)
call
data_format
(
fmt2
,
2
,
&
ifmt1
//
",':',"
//
ifmt1
//
",',',"
//
ifmt1
,
k
,
ifmt3
)
write
(
92
,
fmt2
)
&
"cf_num"
,
j
*
npl
+
1
,
j
*
npl
+
k
,
i
,
icf
(
j
*
npl
+
1
:
j
*
npl
+
k
,
i
)
end do
write
(
92
,
fmt1
)
"cf_den"
,
i
,
idenom
(
i
)
end do
end
c
************************************************************************
c
*
WK
:
calculate
color
flows
in
PYTHIA
-
compatible
format
c
*
c
*
All
particles
are
taken
as
outgoing
(
so
incoming
codes
are
inverted
)
c
*
cflow
:
For
each
particle
,
the
particle
index
where
the
color
comes
from
c
*
aflow
:
For
each
particle
,
the
particle
index
where
the
anticolor
comes
from
c
*
flows
:
The
form
each
factor
in
each
flow
is
stored
by
MADGRAPH
:
c
*
first
the
quark
where
the
color
string
ends
,
c
*
then
the
antiquark
where
the
color
string
begins
,
c
*
then
the
gluons
from
end
to
beginning
of
the
color
string
.
c
*
For
a
pure
gluonic
string
,
the
string
is
implicitly
closed
.
c
************************************************************************
subroutine
setflows
(
nflows
,
cflow
,
aflow
,
flow_c
)
implicit none
! Constants
include
'params.inc'
c
integer
maxlines
,
maxfactors
,
maxflows
c
parameter
(
maxlines
=
9
,
maxfactors
=
10
,
maxflows
=
200
)
c
integer
maxgraphs
c
parameter
(
maxgraphs
=
9999
)
! Arguments
integer
nflows
integer
cflow
(
maxlines
,
maxflows
)
integer
aflow
(
maxlines
,
maxflows
)
integer
flow_c
(
maxflows
,
maxflows
,
2
)
! Local variables
integer
iflow
,
iterm
,
jterm
,
nterms
integer
icolor
,
ncolors
,
nfactors
,
ifact
,
nelements
,
ielem
integer
i
,
icu
,
ifr
,
ito
logical
gluons_only
integer
cflw
(
maxlines
),
aflw
(
maxlines
)
integer
flow_ptr
(
maxterms
)
! Global Variables
integer
flows
(
0
:
2
*
maxlines
,
0
:
maxfactors
,
0
:
24
,
0
:
maxflows
)
integer
graphcolor
(
0
:
2
,
0
:
maxflows
,
0
:
maxgraphs
)
common
/
to_color
/
graphcolor
,
flows
! Initialize flows
nflows
=
0
cflow
=
0
aflow
=
0
flow_c
(:,:,
1
)
=
0
! numerators
flow_c
(:,:,
2
)
=
1
! denominators
! Loop over MADGRAPH color diagrams
ncolors
=
flows
(
0
,
0
,
0
,
0
)
do
icolor
=
1
,
ncolors
cflw
=
0
aflw
=
0
flow_ptr
=
0
! Loop over all terms = flows associated to color diagram
nterms
=
flows
(
0
,
0
,
0
,
icolor
)
do
iterm
=
1
,
nterms
nfactors
=
flows
(
0
,
0
,
iterm
,
icolor
)
do
ifact
=
2
,
nfactors
gluons_only
=
(
flows
(
0
,
ifact
,
iterm
,
icolor
)
>
0
)
nelements
=
abs
(
flows
(
0
,
ifact
,
iterm
,
icolor
))
! Pure gluonic string
if
(
gluons_only
)
then
do
ielem
=
1
,
nelements
icu
=
flows
(
ielem
,
ifact
,
iterm
,
icolor
)
if
(
ielem
<
nelements
)
then
ifr
=
flows
(
ielem
+
1
,
ifact
,
iterm
,
icolor
)
else
ifr
=
flows
(
1
,
ifact
,
iterm
,
icolor
)
end if
if
(
ielem
>
1
)
then
ito
=
flows
(
ielem
-
1
,
ifact
,
iterm
,
icolor
)
else
ito
=
flows
(
nelements
,
ifact
,
iterm
,
icolor
)
end if
cflw
(
icu
)
=
ifr
aflw
(
icu
)
=
ito
end do
! Quark string with gluons attached
else if
(
nelements
>
2
)
then
icu
=
flows
(
1
,
ifact
,
iterm
,
icolor
)
ifr
=
flows
(
3
,
ifact
,
iterm
,
icolor
)
cflw
(
icu
)
=
ifr
icu
=
flows
(
2
,
ifact
,
iterm
,
icolor
)
ito
=
flows
(
nelements
,
ifact
,
iterm
,
icolor
)
aflw
(
icu
)
=
ito
do
ielem
=
3
,
nelements
icu
=
flows
(
ielem
,
ifact
,
iterm
,
icolor
)
if
(
ielem
<
nelements
)
then
ifr
=
flows
(
ielem
+
1
,
ifact
,
iterm
,
icolor
)
else
ifr
=
flows
(
2
,
ifact
,
iterm
,
icolor
)
end if
cflw
(
icu
)
=
ifr
if
(
ielem
>
3
)
then
ito
=
flows
(
ielem
-
1
,
ifact
,
iterm
,
icolor
)
else
ito
=
flows
(
1
,
ifact
,
iterm
,
icolor
)
end if
aflw
(
icu
)
=
ito
end do
! Pure quark string
else if
(
nelements
==
2
)
then
ito
=
flows
(
1
,
ifact
,
iterm
,
icolor
)
ifr
=
flows
(
2
,
ifact
,
iterm
,
icolor
)
cflw
(
ito
)
=
ifr
aflw
(
ifr
)
=
ito
else
stop
"Can't happen: Unterminated color string"
end if
end do
! ifact
! Check if flow already exists
iflow
=
0
do
i
=
1
,
nflows
if
(
all
(
cflw
==
cflow
(:,
i
))
.
and
.
&
all
(
aflw
==
aflow
(:,
i
)))
then
iflow
=
i
exit
end if
end do
if
(
iflow
==
0
)
then
nflows
=
nflows
+
1
iflow
=
nflows
cflow
(:,
iflow
)
=
cflw
aflow
(:,
iflow
)
=
aflw
end if
do
jterm
=
1
,
iterm
-
1
if
(
flow_ptr
(
jterm
)
==
iflow
)
then
stop
&
"Color flow assigned twice to color graph -- can't handle this"
end if
end do
flow_ptr
(
iterm
)
=
iflow
! Enter flow coefficient (num,den) for this term
flow_c
(
icolor
,
iflow
,
1
:
2
)
=
flows
(
1
:
2
,
1
,
iterm
,
icolor
)
end do
! iterm
end do
! icolor
return
end
!***************************************************************************
! This function returns the integer format appropriate for printing ival
! count is the field width as an integer, fmt is the format string
!***************************************************************************
subroutine
get_iformat
(
ival
,
fmt
)
integer
,
intent
(
in
)
::
ival
character
(
3
),
intent
(
out
)
::
fmt
integer
::
i
,
n
n
=
abs
(
ival
)
if
(
n
/=
0
)
then
i
=
0
do while
(
n
/=
0
)
i
=
i
+
1
n
=
n
/
10
end do
else
i
=
1
end if
if
(
ival
<
0
)
i
=
i
+
1
write
(
fmt
,
"('I',I2.2)"
)
i
end
!***************************************************************************
! This function returns the format appropriate for writing a DATA line
!***************************************************************************
subroutine
data_format
&
(
fmt
,
n_skip
,
arg_fmt
,
n_entries
,
entry_fmt
)
character
(
*
),
intent
(
out
)
::
fmt
character
(
*
),
intent
(
in
)
::
arg_fmt
,
entry_fmt
integer
,
intent
(
in
)
::
n_skip
,
n_entries
character
(
3
)
::
ifmt1
character
(
20
)
::
xfmt
,
rfmt
if
(
n_skip
>
0
)
then
call
get_iformat
(
n_skip
,
ifmt1
)
write
(
xfmt
,
"("
//
ifmt1
//
",'x,')"
)
n_skip
else
xfmt
=
""
end if
if
(
n_entries
>
1
)
then
call
get_iformat
(
n_entries
-
1
,
ifmt1
)
write
(
rfmt
,
"(',',"
//
ifmt1
//
",A)"
)
&
n_entries
-
1
,
"(',',"
//
trim
(
entry_fmt
)
//
")"
else
rfmt
=
""
end if
fmt
=
"("
//
trim
(
xfmt
)
//
"'data ',A,'(',"
//
trim
(
arg_fmt
)
//
&
",') /',"
//
trim
(
entry_fmt
)
//
trim
(
rfmt
)
//
",'/')"
end
!***************************************************************************
! *WK: writes out interfacing subroutines
!***************************************************************************
subroutine
write_interface
&
(
nexternal
,
nhel
,
unit_par_def
,
unit_par_set
)
! Constants
include
'params.inc'
c
integer
maxlines
c
parameter
(
maxlines
=
9
)
! Arguments
integer
,
intent
(
in
)
::
nexternal
integer
,
intent
(
in
)
::
nhel
(
maxlines
,
3
**
maxlines
)
integer
,
intent
(
in
)
::
unit_par_def
,
unit_par_set
! Global Variables
integer
iline
(
-
maxlines
:
maxlines
),
idir
(
-
maxlines
:
maxlines
)
integer
this_coup
(
max_coup
)
,
goal_coup
(
max_coup
)
common
/
to_proc
/
iline
,
idir
,
this_coup
,
goal_coup
! character*15 gname
character
*
80
gname
integer
iname
common
/
to_name
/
iname
,
gname
integer
nincoming
common
/
to_proc2
/
nincoming
character
*
(
4
*
max_particles
)
particle
(
4
)
integer
charge_c
(
max_particles
)
integer
iparticle
(
0
:
max_particles
,
0
:
4
),
inverse
(
max_particles
)
common
/
to_model
/
iparticle
,
particle
,
inverse
,
charge_c
integer
pdg_code
(
0
:
max_particles
)
common
/
to_pdg
/
pdg_code
! Local variables
integer
::
i
,
code
! character*20 name
character
*
80
name
character
(
120
)
buff
! Begin code
name
=
gname
(
1
:
iname
)
call
lower_case
(
name
)
! Write service functions
write
(
92
,
11
)
' function number_particles () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = nexternal'
write
(
92
,
11
)
' end function number_particles'
write
(
92
,
11
)
' function number_particles_in () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = nincoming'
write
(
92
,
11
)
' end function number_particles_in'
write
(
92
,
11
)
' function number_particles_out () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = noutgoing'
write
(
92
,
11
)
' end function number_particles_out'
write
(
92
,
10
)
write
(
92
,
11
)
' function number_spin_states () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = ncomb'
write
(
92
,
11
)
' end function number_spin_states'
write
(
92
,
11
)
' function number_spin_states_in () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = ncomb_in'
write
(
92
,
11
)
' end function number_spin_states_in'
write
(
92
,
11
)
' function number_spin_states_out () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = ncomb_out'
write
(
92
,
11
)
' end function number_spin_states_out'
write
(
92
,
10
)
! Assume that the out-spin index varies faster than the in-index
write
(
92
,
11
)
' subroutine spin_states (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
write
(
92
,
11
)
' a = nhel'
write
(
92
,
11
)
' end subroutine spin_states'
write
(
92
,
11
)
' subroutine spin_states_in (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
write
(
92
,
11
)
' a = nhel(:number_particles_in(), '
//
&
'::number_spin_states_out())'
write
(
92
,
11
)
' end subroutine spin_states_in'
write
(
92
,
11
)
' subroutine spin_states_out (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
write
(
92
,
11
)
' a = nhel(number_particles_in()+1:, '
//
&
':number_spin_states_out())'
write
(
92
,
11
)
' end subroutine spin_states_out'
write
(
92
,
10
)
write
(
92
,
11
)
' function number_flavor_states () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = 1'
write
(
92
,
11
)
' end function number_flavor_states'
write
(
92
,
11
)
' subroutine flavor_states (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
do
i
=
1
,
nincoming
code
=
inverse
(
iline
(
i
))
write
(
92
,
23
)
' a('
,
i
,
',1) ='
,
pdg_code
(
code
)
end do
do
i
=
nincoming
+
1
,
nexternal
code
=
iline
(
i
)
write
(
92
,
23
)
' a('
,
i
,
',1) ='
,
pdg_code
(
code
)
end do
write
(
92
,
11
)
' end subroutine flavor_states'
write
(
92
,
10
)
write
(
92
,
11
)
' function number_flavor_states_in () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = 1'
write
(
92
,
11
)
' end function number_flavor_states_in'
write
(
92
,
11
)
' subroutine flavor_states_in (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
do
i
=
1
,
nincoming
code
=
inverse
(
iline
(
i
))
write
(
92
,
23
)
' a('
,
i
,
',1) ='
,
pdg_code
(
code
)
end do
write
(
92
,
11
)
' end subroutine flavor_states_in'
write
(
92
,
10
)
write
(
92
,
11
)
' function number_flavor_states_out () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = 1'
write
(
92
,
11
)
' end function number_flavor_states_out'
write
(
92
,
11
)
' subroutine flavor_states_out (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
do
i
=
1
,
nexternal
-
nincoming
code
=
iline
(
i
+
nincoming
)
write
(
92
,
23
)
' a('
,
i
,
',1) ='
,
pdg_code
(
code
)
end do
write
(
92
,
11
)
' end subroutine flavor_states_out'
write
(
92
,
10
)
write
(
92
,
11
)
' function number_flavor_zeros () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = 0'
write
(
92
,
11
)
' end function number_flavor_zeros'
write
(
92
,
11
)
' subroutine flavor_zeros (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
write
(
92
,
11
)
' a = 0'
write
(
92
,
11
)
' end subroutine flavor_zeros'
write
(
92
,
10
)
write
(
92
,
11
)
' function number_flavor_zeros_in () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = 0'
write
(
92
,
11
)
' end function number_flavor_zeros_in'
write
(
92
,
11
)
' subroutine flavor_zeros_in (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
write
(
92
,
11
)
' a = 0'
write
(
92
,
11
)
' end subroutine flavor_zeros_in'
write
(
92
,
10
)
write
(
92
,
11
)
' function number_flavor_zeros_out () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = 0'
write
(
92
,
11
)
' end function number_flavor_zeros_out'
write
(
92
,
11
)
' subroutine flavor_zeros_out (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
write
(
92
,
11
)
' a = 0'
write
(
92
,
11
)
' end subroutine flavor_zeros_out'
write
(
92
,
10
)
write
(
92
,
11
)
' function number_color_flows () result (n)'
write
(
92
,
11
)
' integer :: n'
write
(
92
,
11
)
' n = nflows'
write
(
92
,
11
)
' end function number_color_flows'
write
(
92
,
11
)
' subroutine color_flows (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
write
(
92
,
11
)
' a = cflw'
write
(
92
,
11
)
' end subroutine color_flows'
write
(
92
,
11
)
' subroutine anticolor_flows (a)'
write
(
92
,
11
)
' integer, dimension(:,:), intent(inout) :: a'
write
(
92
,
11
)
' a = aflw'
write
(
92
,
11
)
' end subroutine anticolor_flows'
write
(
92
,
10
)
! Write the constructor for the constant block
! This initializes the REAL representation for the color data
! The multiplicity is renormalized by the number of in- spin states
! because this is taken into account by the in-state density matrix
write
(
92
,
10
)
" subroutine create"
write
(
92
,
10
)
" integer :: i"
write
(
92
,
10
)
" do i=1, ncolors"
write
(
92
,
10
)
" cf(:,i) = real(cf_num(:,i),D) / cf_den(i)"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" flow_c = real(flow_c_num,D) / flow_c_den"
write
(
92
,
10
)
" cmult = real(iden,D) / real(ncomb_in,D)"
write
(
92
,
10
)
" end subroutine create"
write
(
92
,
10
)
! Write the destructor for the constant block
! Nothing to do since no dynamical memory is used
write
(
92
,
10
)
" subroutine destroy"
write
(
92
,
10
)
" end subroutine destroy"
write
(
92
,
10
)
! Write the routine for setting constants:
! transfer from tmp file
! This replaces the corresponding helas routines
rewind
(
unit_par_set
)
do
read
(
unit_par_set
,
'(A)'
,
end
=
70
)
buff
write
(
92
,
'(A)'
)
trim
(
buff
)
cycle
70
exit
end do
! Write the sqme wrapper for general spin density matrices
! Assume that the out-spin index varies faster than the in-index
write
(
92
,
10
)
write
(
92
,
10
)
" subroutine scatter_nonzero "
//
&
"(p, rho_in, rho_out, zero_ct, n)"
write
(
92
,
10
)
" real(D), dimension(0:,:), intent(in) :: p"
write
(
92
,
10
)
" "
//
&
"complex(D), dimension(:,:,:,:), intent(in) :: rho_in"
write
(
92
,
10
)
" "
//
&
"complex(D), dimension(:,:,:,:), intent(inout) :: rho_out"
write
(
92
,
10
)
" "
//
&
"integer, dimension(:,:,:,:), intent(inout) :: zero_ct"
write
(
92
,
10
)
" integer, intent(in) :: n"
write
(
92
,
10
)
" integer :: hi1, ho1, hi2, ho2, h1, h2"
write
(
92
,
10
)
" "
//
&
"complex(D), dimension(ncolors) :: amp1, amp2"
write
(
92
,
10
)
" "
//
&
"call calculate_amplitudes "
//
&
"(p, zero_ct, n, calc_flows=.false.)"
write
(
92
,
10
)
" rho_out = 0"
write
(
92
,
10
)
" h1 = 0"
write
(
92
,
10
)
" do hi1 = 1, number_spin_states_in ()"
write
(
92
,
10
)
" do ho1 = 1, number_spin_states_out ()"
write
(
92
,
10
)
" h1 = h1 + 1"
write
(
92
,
10
)
" h2 = 0"
write
(
92
,
10
)
" do hi2 = 1, number_spin_states_in ()"
write
(
92
,
10
)
" do ho2 = 1, number_spin_states_out ()"
write
(
92
,
10
)
" h2 = h2 + 1"
write
(
92
,
10
)
" amp1 = "
//
&
"matmul (amp_color(:,h1,1), cf)"
write
(
92
,
10
)
" amp2 = amp_color(:,h2,1)"
write
(
92
,
10
)
" rho_out(ho1,1,ho2,1) = "
//
&
"rho_out(ho1,1,ho2,1) &"
write
(
92
,
10
)
" "
//
&
" & + dot_product (amp1, amp2) * rho_in(hi1,1,hi2,1) &"
write
(
92
,
10
)
" "
//
&
" & / cmult(1)"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end subroutine scatter_nonzero"
write
(
92
,
10
)
write
(
92
,
10
)
! Write the sqme wrapper for diagonal density matrices
! Again, assume that the out-spin index varies faster than the in-index
write
(
92
,
10
)
" subroutine scatter_diagonal_nonzero "
//
&
"(p, rho_in, rho_out, zero_ct, n)"
write
(
92
,
10
)
" real(D), dimension(0:,:), intent(in) :: p"
write
(
92
,
10
)
" real(D), dimension(:,:), intent(in) :: rho_in"
write
(
92
,
10
)
" "
//
&
"real(D), dimension(:,:), intent(inout) :: rho_out"
write
(
92
,
10
)
" "
//
&
"integer, dimension(:,:,:,:), intent(inout) :: zero_ct"
write
(
92
,
10
)
" integer, intent(in) :: n"
write
(
92
,
10
)
" integer :: hi, ho, h"
write
(
92
,
10
)
" "
//
&
"complex(D), dimension(ncolors) :: amp1, amp2"
write
(
92
,
10
)
" "
//
&
"call calculate_amplitudes "
//
&
"(p, zero_ct, n, calc_flows=.false.)"
write
(
92
,
10
)
" rho_out = 0"
write
(
92
,
10
)
" h = 0"
write
(
92
,
10
)
" do hi = 1, number_spin_states_in ()"
write
(
92
,
10
)
" do ho = 1, number_spin_states_out ()"
write
(
92
,
10
)
" h = h + 1"
write
(
92
,
10
)
" amp1 = matmul (amp_color(:,h,1), cf)"
write
(
92
,
10
)
" amp2 = amp_color(:,h,1)"
write
(
92
,
10
)
" rho_out(ho,1) = rho_out(ho,1) &"
write
(
92
,
10
)
" "
//
&
" & + dot_product (amp1, amp2) * rho_in(hi,1) &"
write
(
92
,
10
)
" "
//
&
" & / cmult(1)"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end subroutine scatter_diagonal_nonzero"
write
(
92
,
10
)
write
(
92
,
10
)
! Write the sqme wrapper for colored amplitudes (generic spin density)
! Again, assume that the out-spin index varies faster than the in-index
! The incoming color is assumed uncorrelated, therefore rho_in is
! four-dimensional (hel1, flv1, hel2, flv2),
! The outgoing color flow information is stored in rho_col_out
! which therefore is five-dimensional. The color-summed result
! is returned in the four-dimensional array rho_out.
write
(
92
,
10
)
" subroutine scatter_colored_nonzero &"
write
(
92
,
10
)
" & (p, rho_in, rho_out, rho_col_out, "
//
&
"zero_ct, n)"
write
(
92
,
10
)
" real(D), dimension(0:,:), intent(in) :: p"
write
(
92
,
10
)
" "
//
&
"complex(D), dimension(:,:,:,:), intent(in) :: rho_in"
write
(
92
,
10
)
" "
//
&
"complex(D), dimension(:,:,:,:), intent(inout) :: rho_out"
write
(
92
,
10
)
" "
//
&
"complex(D), dimension(:,:,:,:,:), intent(inout) :: "
//
&
"rho_col_out"
write
(
92
,
10
)
" "
//
&
"integer, dimension(:,:,:,:), intent(inout) :: zero_ct"
write
(
92
,
10
)
" integer, intent(in) :: n"
write
(
92
,
10
)
" integer :: hi1, ho1, hi2, ho2, h1, h2, c"
write
(
92
,
10
)
" "
//
&
"complex(D), dimension(ncolors) :: amp1, amp2"
write
(
92
,
10
)
" "
//
&
"call calculate_amplitudes "
//
&
"(p, zero_ct, n, calc_flows=.true.)"
write
(
92
,
10
)
" rho_out = 0"
write
(
92
,
10
)
" h1 = 0"
write
(
92
,
10
)
" do hi1 = 1, number_spin_states_in ()"
write
(
92
,
10
)
" do ho1 = 1, number_spin_states_out ()"
write
(
92
,
10
)
" h1 = h1 + 1"
write
(
92
,
10
)
" h2 = 0"
write
(
92
,
10
)
" do hi2 = 1, number_spin_states_in ()"
write
(
92
,
10
)
" do ho2 = 1, number_spin_states_out ()"
write
(
92
,
10
)
" h2 = h2 + 1"
write
(
92
,
10
)
" amp1 = "
//
&
"matmul (amp_color(:,h1,1), cf)"
write
(
92
,
10
)
" amp2 = amp_color(:,h2,1)"
write
(
92
,
10
)
" rho_out(ho1,1,ho2,1) = "
//
&
"rho_out(ho1,1,ho2,1) &"
write
(
92
,
10
)
" "
//
&
" & + dot_product (amp1, amp2) * rho_in(hi1,1,hi2,1) &"
write
(
92
,
10
)
" "
//
&
" & / cmult(1)"
write
(
92
,
10
)
" do c = 1, nflows"
write
(
92
,
10
)
" rho_col_out(c,ho1,1,ho2,1) = "
//
&
"rho_col_out(c,ho1,1,ho2,1) &"
write
(
92
,
10
)
" "
//
&
" & + amp_flow(c,h1,1) * conjg(amp_flow(c,h2,1)) &"
write
(
92
,
10
)
" "
//
&
" & * rho_in(hi1,1,hi2,1)"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end subroutine scatter_colored_nonzero"
write
(
92
,
10
)
write
(
92
,
10
)
! Write the sqme wrapper for colored amplitudes with diagonal spin density
! Again, assume that the out-spin index varies faster than the in-index
! The color flow index is the first index in rho_out which therefore
! is three-dimensional. For counting zeros, the color index is summed over.
write
(
92
,
10
)
" subroutine scatter_diagonal_colored_nz "
//
&
"(p, rho_in, rho_out, zero_ct, n)"
write
(
92
,
10
)
" real(D), dimension(0:,:), intent(in) :: p"
write
(
92
,
10
)
" real(D), dimension(:,:), intent(in) :: rho_in"
write
(
92
,
10
)
" "
//
&
"real(D), dimension(:,:,:), intent(inout) :: rho_out"
write
(
92
,
10
)
" "
//
&
"integer, dimension(:,:,:,:), intent(inout) :: zero_ct"
write
(
92
,
10
)
" integer, intent(in) :: n"
write
(
92
,
10
)
" integer :: hi, ho, h, c"
write
(
92
,
10
)
" real(D), dimension(nflows,ncross) :: pcol"
write
(
92
,
10
)
" real(D), dimension(ncross) :: sqme, sum_pcol"
write
(
92
,
10
)
" "
//
&
"complex(D), dimension(ncolors) :: amp1, amp2"
write
(
92
,
10
)
" "
//
&
"call calculate_amplitudes "
//
&
"(p, zero_ct, n, calc_flows=.true.)"
write
(
92
,
10
)
" rho_out = 0"
write
(
92
,
10
)
" h = 0"
write
(
92
,
10
)
" do hi = 1, number_spin_states_in ()"
write
(
92
,
10
)
" do ho = 1, number_spin_states_out ()"
write
(
92
,
10
)
" h = h + 1"
write
(
92
,
10
)
" amp1 = matmul (amp_color(:,h,1), cf)"
write
(
92
,
10
)
" amp2 = amp_color(:,h,1)"
write
(
92
,
10
)
" "
//
&
"sqme(1) = dot_product (amp1, amp2) * rho_in(hi,1) &"
write
(
92
,
10
)
" "
//
&
" & / cmult(1)"
write
(
92
,
10
)
" do c = 1, nflows"
write
(
92
,
10
)
" pcol(c,1) = abs (amp_flow(c,h,1))**2"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" sum_pcol(1) = sum (pcol(:,1))"
write
(
92
,
10
)
" if (sum_pcol(1) /= 0) then"
write
(
92
,
10
)
" rho_out(:,ho,1) = "
//
&
"rho_out(:,ho,1) &"
write
(
92
,
10
)
" "
//
&
" & + sqme(1) * pcol(:,1) / sum_pcol(1)"
write
(
92
,
10
)
" end if"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end do"
write
(
92
,
10
)
" end subroutine scatter_diagonal_colored_nz"
write
(
92
,
10
)
write
(
92
,
10
)
10
format
(
a
)
11
format
(
a
,
i1
)
12
format
(
a
,
i1
,
a
)
13
format
(
a
,
i1
,
a
,
i3
)
22
format
(
a
,
i2
,
a
,
i2
,
a
)
23
format
(
a
,
i2
,
a
,
i9
)
44
format
(
a
,
i4
)
c
10
format
(
6
x
,
a
)
c
11
format
(
6
x
,
a
,
i1
)
c
12
format
(
6
x
,
a
,
i1
,
a
)
c
13
format
(
6
x
,
a
,
i1
,
a
,
i3
)
c
22
format
(
6
x
,
a
,
i2
,
a
,
i2
,
a
)
c
23
format
(
6
x
,
a
,
i2
,
a
,
i9
)
c
44
format
(
6
x
,
a
,
i4
)
end subroutine
write_interface
File Metadata
Details
Attached
Mime Type
text/plain
Expires
Mon, Jan 20, 8:42 PM (22 h, 54 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242323
Default Alt Text
write_module.f (53 KB)
Attached To
rWHIZARDSVN whizardsvn
Event Timeline
Log In to Comment