Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F11221726
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
202 KB
Subscribers
None
View Options
Index: trunk/src/fks/fks.nw
===================================================================
--- trunk/src/fks/fks.nw (revision 6497)
+++ trunk/src/fks/fks.nw (revision 6498)
@@ -1,5292 +1,5367 @@
% -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*-
% WHIZARD code as NOWEB source: matrix elements and process libraries
@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{FKS Subtraction Scheme}
The code in this chapter implements the FKS subtraction scheme for use
with \whizard.
These are the modules:
\begin{description}
\item[fks\_regions]
Given a process definition, identify singular regions in the
associated phase space.
\item[virtual]
Handle the virtual correction matrix element.
\item[real\_subtraction]
Handle the real-subtraction matrix element.
\item[nlo\_data]
Manage the subtraction objects.
\item[phs\_fks]
Phase-space parameterization with modifications for the FKS scheme.
\end{description}
This chapter deals with next-to-leading order contributions to cross sections.
Basically, there are three major issues to be adressed: The creation
of the $N+1$-particle flavor structure, the construction of the
$N+1$-particle phase space and the actual calculation of the real- and
virtual-subtracted matrix elements. The first is dealt with using the
[[auto_components]] class, and it will be shown that the second
and third issue are connected in FKS subtraction.
\section{Brief outline of FKS subtraction}
{\em In the current state, this discussion is only concerned with
lepton collisions. For hadron collisions, renormalization of parton
distributions has to be taken into account. Further, for QCD
corrections, initial-state radiation is necessarily
present. However, most quantities have so far been only constructed
for final-state emissions}
The aim is to calculate the next-to-leading order cross section
according to
\begin{equation*}
d\sigma_{\rm{NLO}} = \mathcal{B} + \mathcal{V} +
\mathcal{R}d\Phi_{\rm{rad}}.
\end{equation*}
Analytically, the divergences, in terms of poles in the complex
quantity $\varepsilon = 2-d/2$, cancel. However, this is in general
only valid in an arbitrary, comlex number of dimensions. This is,
roughly, the content of the KLN-theorem. \whizard, as any
other numerical program, is confined to four dimensions. We will
assume that the KLN-theorem is valid and that there exist subtraction
terms $\mathcal{C}$ such that
\begin{equation*}
d\sigma_{\rm{NLO}} = \mathcal{B} + \underbrace{\mathcal{V} +
\mathcal{C}}_{\text{finite}} + \underbrace{\mathcal{R} -
\mathcal{C}}_{\text{finite}},
\end{equation*}
i.e. the subtraction terms correspond to the divergent limits of the
real and virtual matrix element.
Because $\mathcal{C}$ subtracts the divergences of $\mathcal{R}$ as
well as those of $\mathcal{V}$, it suffices to consider one of them,
so we focus on $\mathcal{R}$. For this purpose, $\mathcal{R}$ is
rewritten,
\begin{equation*}
\mathcal{R} = \frac{1}{\xi^2}\frac{1}{1-y} \left(\xi^2
(1-y)\mathcal{R}\right) =
\frac{1}{\xi^2}\frac{1}{1-y}\tilde{\mathcal{R}},
\end{equation*}
with $\xi = \left(2k_{\rm{rad}}^0\right)/\sqrt{s}$ and $y =
\cos\theta$, where $k_{\rm{rad}}^0$ denotes the energy of the radiated
parton and $\theta$ is the angle between emitter and radiated
parton. $\tilde{\mathcal{R}}$ is finite, therefore the whole
singularity structure is contained in the prefactor
$\xi^{-2}(1-y)^{-1}$. Combined with the d-dimensional phase space
element,
\begin{equation*}
\frac{d^{d-1}k}{2k^0(2\pi)^{d-1}} =
\frac{s^{1-\varepsilon}}{(4\pi)^{d-1}}\xi^{1-2\varepsilon}\left(1-y^2\right)^{-\varepsilon}
d\xi dy d\Omega^{d-2},
\end{equation*}
this yields
\begin{equation*}
d\Phi_{\rm{rad}} \mathcal{R} = dy (1-y)^{-1-\varepsilon} d\xi
\xi^{-1-2\varepsilon} \tilde{R}.
\end{equation*}
This can further be rewritten in terms of plus-distributions,
\begin{align*}
\xi^{-1-2\varepsilon} &= -\frac{1}{2\varepsilon}\delta(\xi) +
\left(\frac{1}{\xi}\right)_+ -
2\varepsilon\left(\frac{\log\xi}{\xi}\right)_+ +
\mathcal{O}(\varepsilon^2),\\
(1-y)^{-1-\varepsilon} &= -\frac{2^{-\varepsilon}}{\varepsilon}
\delta(1-y) + \left(\frac{1}{1-y}\right)_+ - \varepsilon
\left(\frac{1}{1-y}\right)_+\log(1-y) + \mathcal{O}(\varepsilon^2),
\end{align*}
(imagine that all this is written inside of integrals, which are
spared for ease of notation) such that
\begin{align*}
d\Phi_{\rm{rad}} \mathcal{R} &= -\frac{1}{2\varepsilon} dy
(1-y)^{-1-\varepsilon}\tilde{R} (0,y) -
d\xi\left[\frac{2^{-\varepsilon}}{\varepsilon}\left(\frac{1}{\xi}\right)_+
- 2\left(\frac{\log\xi}{\xi}\right)_+\right] \tilde{R}(\xi,1) \\
&+ dy d\xi \left(\frac{1}{\xi}\right)_+
\left(\frac{1}{1-y}\right)_+
\tilde{R}(\xi, y) +
\mathcal{O}(\varepsilon).\\
\end{align*}
The summand in the second line is of order $\mathcal{O}(1)$ and is the
only one to reproduce $\mathcal{R}{\xi,y}$. It thus constitutes the
sum of the real matrix element and the corresponding counterterms.
The first summand consequently consists of the subtraction terms to
the virtual matrix elements. Above formula thus allows to calculate
all quantities to render the matrix elements finite.
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Identifying singular regions}
In the FKS subtraction scheme, the phase space is decomposed into
disjoint singular regions, such that
\begin{equation}
\label{eq:S_complete}
\sum_i \mathcal{S}_i + \sum_{ij}\mathcal{S}_{ij} = 1.
\end{equation}
The quantities $\mathcal{S}_i$ and $\mathcal{S}_{ij}$ are functions of
phase space corresponding to a pair of particles indices which can
make up a divergent phase space region. We call such an index pair a
fundamental tuple. For example, the process $e^+ \, e^- \rightarrow u
\, \bar{u} \, g$ has two singular regions, $(3,5)$ and $(4,5)$,
indicating that the gluon can be soft or collinear with respect to
either the quark or the anti-quark. Therefore, the functions $S_{ij}$
have to be chosen in such a way that their contribution makes up most
of \eqref{eq:S_complete} in phase-space configurations where
(final-state) particle $j$ is collinear to particle $i$ or/and
particle $j$ is soft. The functions $S_i$ is the corresponding
quantity for initial-state divergences.
As a singular region we understand the collection of real flavor
structures associated with an emitter and a list of all possible
fundamental tuples. As an example, consider the process $e^+ \, e^-
\rightarrow u \, \bar{u} \, g$. At next-to-leading order, processes
with an additionally radiated particle have to be considered. In this
case, these are $e^+ \, e^- \rightarrow u \, \bar{u}, \, g \, g$,
and $e^+ \, e^- \rightarrow u \, \bar{u} \, u \, \bar{u}$ (or the same
process with any other quark). Table \ref{table:singular regions} sums
up all possible singular regions for this problem.
\begin{table}
\label{table:singular regions}
\begin{tabular}{|c|c|c|c|}
\hline
\texttt{alr} & \texttt{flst\_alr} & \texttt{emi} &
\texttt{ftuple\_list}\\ \hline
1 & [-11,11,2,-2,21,21] & 3 & {(3,5), (3,6), (4,5), (4,6), (5,6)} \\ \hline
2 & [-11,11,2,-2,21,21] & 4 & {(3,5), (3,6), (4,5), (4,6), (5,6)} \\ \hline
3 & [-11,11,2,-2,21,21] & 5 & {(3,5), (3,6), (4,5), (4,6), (5,6)} \\ \hline
4 & [-11,11,2,-2,2,-2] & 5 & {(5,6)} \\
\hline
\end{tabular}
\caption{List of singular regions. The particles are represented by
their PDG codes. The third column contains the emitter for the
specific singular region. For the process involving an additional
gluon, the gluon can either be emitted from one of the quarks or
from the first gluon. Each emitter yields the same list of
fundamental tuples, five in total. The last singular region
corresponds to the process where the gluon splits up into two
quarks. Here, there is only one fundamental tuple, corresponding to
a singular configuration of the momenta of the additional quarks.}
\end{table}
\\
\begin{table}
\label{table:ftuples and flavors}
\begin{tabular}{|c|c|c|c|}
\hline
\texttt{alr} & \texttt{ftuple} & \texttt{emitter} &
\texttt{flst\_alr} \\ \hline
1 & $(3,5)$ & 5 & [-11,11,-2,21,2,21] \\ \hline
2 & $(4,5)$ & 5 & [-11,11,2,21,-2,21] \\ \hline
3 & $(3,6)$ & 5 & [-11,11,-2,21,2,21] \\ \hline
4 & $(4,6)$ & 5 & [-11,11,2,21,-2,21] \\ \hline
5 & $(5,6)$ & 5 & [-11,11,2,-2,21,21] \\ \hline
6 & $(5,6)$ & 5 & [-11,11,2,-2,2,-2] \\ \hline
\end{tabular}
\caption{Initial list of singular regions}
\end{table}
Thus, during the preparation of a NLO-calculation, the possible
singular regions have to be identified. [[fks_regions.f90]] deals
with this issue.
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{FKS Regions}
<<[[fks_regions.f90]]>>=
<<File header>>
module fks_regions
use kinds
use io_units
<<Use strings>>
use constants
use diagnostics
use flavors
use process_constants
use lorentz
use pdg_arrays
use model_data
use physics_defs
<<Standard module head>>
<<fks regions: public>>
<<fks regions: parameters>>
<<fks regions: types>>
<<fks regions: interfaces>>
contains
<<fks regions: procedures>>
end module fks_regions
@ %def fks_regions
@ There are three fundamental splitting types: $q \rightarrow qg$, $g \rightarrow gg$ and
$g \rightarrow qq$.
<<fks regions: parameters>>=
integer, parameter :: UNDEFINED_SPLITTING = 0
integer, parameter :: Q_TO_QG = 1
integer, parameter :: G_TO_GG = 2
integer, parameter :: G_TO_QQ = 3
@
@ We group the indices of the emitting and the radiated particle in
the [[ftuple]]-object.
<<fks regions: public>>=
public :: ftuple_t
<<fks regions: types>>=
type :: ftuple_t
integer, dimension(2) :: ireg
integer :: splitting_type
contains
<<fks regions: ftuple: TBP>>
end type ftuple_t
@ %def ftuple_t
@
<<fks regions: ftuple: TBP>>=
procedure :: write => ftuple_write
<<fks regions: procedures>>=
subroutine ftuple_write (ftuple, unit)
class(ftuple_t), intent(in) :: ftuple
integer, intent(in), optional :: unit
integer :: u
u = given_output_unit (unit); if (u < 0) return
write (u, "(A1,I1,A1,I1,A1)") &
'(', ftuple%ireg(1), ',', ftuple%ireg(2), ')'
end subroutine ftuple_write
@ %def ftuple_write
@
<<fks regions: ftuple: TBP>>=
procedure :: get => ftuple_get
<<fks regions: procedures>>=
subroutine ftuple_get (ftuple, pos1, pos2)
class(ftuple_t), intent(in) :: ftuple
integer, intent(out) :: pos1, pos2
pos1 = ftuple%ireg(1)
pos2 = ftuple%ireg(2)
end subroutine ftuple_get
@ %def ftuple_get
@
<<fks regions: ftuple: TBP>>=
procedure :: set => ftuple_set
<<fks regions: procedures>>=
subroutine ftuple_set (ftuple, pos1, pos2)
class(ftuple_t), intent(inout) :: ftuple
integer, intent(in) :: pos1, pos2
ftuple%ireg(1) = pos1
ftuple%ireg(2) = pos2
end subroutine ftuple_set
@ %def ftuple_set
@
<<fks regions: ftuple: TBP>>=
procedure :: determine_splitting_type => ftuple_determine_splitting_type
<<fks regions: procedures>>=
subroutine ftuple_determine_splitting_type (ftuple, flv, i, j)
class(ftuple_t), intent(inout) :: ftuple
type(flv_structure_t), intent(in) :: flv
integer, intent(in) :: i, j
associate (flst => flv%flst)
if (flst(i) == GLUON .and. flst(j) == GLUON) then
ftuple%splitting_type = G_TO_GG
else if (flst(i)+flst(j) == 0 &
.and. is_quark (abs(flst(i)))) then
ftuple%splitting_type = G_TO_QQ
else if (is_quark(abs(flst(i))) .and. flst(j) == GLUON &
.or. is_quark(abs(flst(j))) .and. flst(i) == GLUON) then
ftuple%splitting_type = Q_TO_QG
else
ftuple%splitting_type = UNDEFINED_SPLITTING
end if
end associate
end subroutine ftuple_determine_splitting_type
@ %def ftuple_determine_splitting_type
@ Each singular region can have a different number of
emitter-radiation pairs. This is coped with using the linked list
[[ftuple_list]].
<<fks regions: types>>=
type :: ftuple_list_t
integer :: index
type(ftuple_t) :: ftuple
type(ftuple_list_t), pointer :: next
type(ftuple_list_t), pointer :: prev
type(ftuple_list_t), pointer :: equiv
contains
<<fks regions: ftuple list: TBP>>
end type ftuple_list_t
@ %def ftuple_list_t
@
<<fks regions: ftuple list: TBP>>=
procedure :: init => ftuple_list_init
<<fks regions: procedures>>=
subroutine ftuple_list_init (list)
class(ftuple_list_t), intent(inout) :: list
list%index = 0
nullify (list%next)
nullify (list%prev)
nullify (list%equiv)
end subroutine ftuple_list_init
@ %def ftuple_list_init
@
<<fks regions: ftuple list: TBP>>=
procedure :: write => ftuple_list_write
<<fks regions: procedures>>=
subroutine ftuple_list_write (list)
class(ftuple_list_t), intent(in), target :: list
type(ftuple_list_t), pointer :: current
select type (list)
type is (ftuple_list_t)
current => list
do
call current%ftuple%write
if (associated (current%next)) then
current => current%next
else
exit
end if
end do
end select
end subroutine ftuple_list_write
@ %def ftuple_list_write
@
<<fks regions: ftuple list: TBP>>=
procedure :: append => ftuple_list_append
<<fks regions: procedures>>=
subroutine ftuple_list_append (list, ftuple)
class(ftuple_list_t), intent(inout), target :: list
type(ftuple_t), intent(in) :: ftuple
type(ftuple_list_t), pointer :: current
select type (list)
type is (ftuple_list_t)
if (list%index == 0) then
nullify(list%next)
list%index = 1
list%ftuple = ftuple
else
current => list
do
if (associated (current%next)) then
current => current%next
else
allocate (current%next)
nullify (current%next%next)
nullify (current%next%equiv)
current%next%prev => current
current%next%index = current%index + 1
current%next%ftuple = ftuple
exit
end if
end do
end if
end select
end subroutine ftuple_list_append
@ %def ftuple_list_append
@
<<fks regions: ftuple list: TBP>>=
procedure :: get_n_tuples => ftuple_list_get_n_tuples
<<fks regions: procedures>>=
function ftuple_list_get_n_tuples (list) result(n_tuples)
class(ftuple_list_t), intent(inout), target :: list
integer :: n_tuples
type(ftuple_list_t), pointer :: current
select type (list)
type is (ftuple_list_t)
current => list
n_tuples = 1
do
if (associated (current%next)) then
current => current%next
n_tuples = n_tuples + 1
else
exit
end if
end do
end select
end function ftuple_list_get_n_tuples
@ %def ftuple_list_get_n_tuples
@
<<fks regions: ftuple list: TBP>>=
procedure :: get_entry => ftuple_list_get_entry
<<fks regions: procedures>>=
function ftuple_list_get_entry(list, index) result(entry)
class(ftuple_list_t), intent(inout), target :: list
integer, intent(in) :: index
type(ftuple_list_t), pointer :: entry
type(ftuple_list_t), pointer :: current
integer :: i
select type (list)
type is (ftuple_list_t)
current => list
if (index <= list%get_n_tuples ()) then
if (index == 1) then
entry => current
else
do i=1,index-1
current => current%next
end do
entry => current
end if
else
call msg_fatal &
("Index must be smaller or equal than the total number of regions!")
end if
end select
end function ftuple_list_get_entry
@ %def ftuple_list_get_entry
@
<<fks regions: ftuple list: TBP>>=
procedure :: get_ftuple => ftuple_list_get_ftuple
<<fks regions: procedures>>=
function ftuple_list_get_ftuple (list, index) result (ftuple)
class(ftuple_list_t), intent(inout) :: list
integer, intent(in) :: index
type(ftuple_t) :: ftuple
type(ftuple_list_t) :: entry
entry = list%get_entry (index)
ftuple = entry%ftuple
end function ftuple_list_get_ftuple
@ %def ftuple_list_get_ftuple
@
<<fks regions: ftuple list: TBP>>=
procedure :: set_equiv => ftuple_list_set_equiv
<<fks regions: procedures>>=
subroutine ftuple_list_set_equiv (list, i1, i2)
class(ftuple_list_t), intent(inout) :: list
integer, intent(in) :: i1, i2
type(ftuple_list_t), pointer :: list1, list2
select type (list)
type is (ftuple_list_t)
list1 => list%get_entry (i1)
list2 => list%get_entry (i2)
do
if (associated (list1%equiv)) then
list1 => list1%equiv
else
exit
end if
end do
list1%equiv => list2
end select
end subroutine ftuple_list_set_equiv
@ %def ftuple_list_set_equiv
@
<<fks regions: ftuple list: TBP>>=
procedure :: check_equiv => ftuple_list_check_equiv
<<fks regions: procedures>>=
function ftuple_list_check_equiv(list, i1, i2) result(eq)
class(ftuple_list_t), intent(inout) :: list
integer, intent(in) :: i1, i2
logical :: eq
type(ftuple_list_t), pointer :: current
select type (list)
type is (ftuple_list_t)
current => list%get_entry (i1)
do
if (associated (current%equiv)) then
current => current%equiv
if (current%index == i2) then
eq = .true.
exit
end if
else
eq = .false.
exit
end if
end do
end select
end function ftuple_list_check_equiv
@ %def ftuple_list_check_equiv
@ Class for working with the flavor specification arrays.
<<fks regions: public>>=
public :: flv_structure_t
<<fks regions: types>>=
type :: flv_structure_t
integer, dimension(:), allocatable :: flst
integer :: nlegs
logical, dimension(:), allocatable :: massive
logical, dimension(:), allocatable :: colored
contains
<<fks regions: flv structure: TBP>>
end type flv_structure_t
@ %def flv_structure_t
@
Returns \texttt{true} if the two particles at position \texttt{i}
and \texttt{j} in the flavor array can originate from the same
splitting. For this purpose, the function first checks whether the splitting is
allowed at all. If this is the case, the emitter is removed from the
flavor array. If the resulting array is equivalent to the Born flavor
structure \texttt{flv\_born}, the pair is accepted as a valid
splitting.
<<fks regions: flv structure: TBP>>=
procedure :: valid_pair => flv_structure_valid_pair
<<fks regions: procedures>>=
function flv_structure_valid_pair &
(flv_real,i,j, flv_born, model) result (valid)
class(flv_structure_t), intent(in) :: flv_real
integer, intent(in) :: i,j
type(flv_structure_t), intent(in) :: flv_born
type(model_data_t), intent(in) :: model
logical :: valid
integer :: k, n_orig
type(flv_structure_t) :: flv_test
integer, dimension(:), allocatable :: flv_orig, flv_orig2
valid = .false.
@
First check whether the splitting is possible. The array
[[flv_orig]] contains all particles which share a vertex with the
particles at position [[i]] and [[j]]. If its size is equal to zero,
no splitting is possible and the subroutine is exited.
<<fks regions: procedures>>=
call model%match_vertex &
(flv_real%flst(i), flv_real%flst(j), flv_orig)
n_orig = size (flv_orig)
if (n_orig == 0) then
return
else
@
For a quark emitting a gluon, [[flv_orig]] contains the PDG code of
the anti-quark. To be on the safe side, a second array is created,
which contains both the positively and negatively signed PDG
codes. Then, the origial tuple $(i,j)$ is removed from the real flavor
structure and the particles in [[flv_orig2]] are inserted.
If the resulting Born configuration is equal to the underlying Born
configuration, up to a permutation of final-state particles, the tuple
$(i,j)$ is accepted as valid.
<<fks regions: procedures>>=
allocate (flv_orig2 (2*n_orig))
flv_orig2 (1:n_orig) = flv_orig
flv_orig2 (n_orig+1:2*n_orig) = -flv_orig
do k = 1, 2*n_orig
flv_test = flv_real%insert_particle (i,j,flv_orig2(k))
valid = flv_born == flv_test
if (valid) return
end do
end if
end function flv_structure_valid_pair
@ %def flv_structure_valid_pair
@ This function checks whether two flavor arrays are the same up to a
permutation of the final-state particles
<<fks regions: procedures>>=
function flv_structure_equivalent (flv1, flv2) result(equiv)
type(flv_structure_t), intent(in) :: flv1, flv2
logical :: equiv
integer :: i, j, n
integer :: f1, f2
logical, dimension(:), allocatable :: present, checked
n = size (flv1%flst)
equiv = .true.
if (n /= size (flv2%flst)) then
call msg_fatal &
('flv_structure_equivalent: flavor arrays do not have equal lengths')
else
allocate (present(n))
allocate (checked(n))
do i=1,n
present(i) = .false.
checked(i) = .false.
end do
do i=1,n
do j=1,n
if (flv1%flst(i) == flv2%flst(j) .and. .not. checked(j)) then
present(i) = .true.
checked(j) = .true.
exit
end if
end do
end do
do i=1,n
if(.not.present(i)) equiv = .false.
end do
end if
end function flv_structure_equivalent
@ %def flv_structure_equivalent
@ Returs a new flavor array with the particle at position
\texttt{index} removed.
<<fks regions: flv structure: TBP>>=
procedure :: remove_particle => flv_structure_remove_particle
<<fks regions: procedures>>=
function flv_structure_remove_particle (flv1, index) result(flv2)
class(flv_structure_t), intent(in) :: flv1
integer, intent(in) :: index
type(flv_structure_t) :: flv2
integer :: n1, n2
n1 = size (flv1%flst)
n2 = n1-1
if (allocated (flv2%flst)) then
deallocate (flv2%flst)
end if
allocate (flv2%flst (n2))
if (index == 1) then
flv2%flst(1:n2) = flv1%flst(2:n1)
else if (index == n1) then
flv2%flst(1:n2) = flv1%flst(1:n2)
else
flv2%flst(1:index-1) = flv1%flst(1:index-1)
flv2%flst(index:n2) = flv1%flst(index+1:n1)
end if
end function flv_structure_remove_particle
@ %def flv_structure_remove_particle
@ Removes the paritcles at position i1 and i2 and inserts a new
particle at position i1.
<<fks regions: flv structure: TBP>>=
procedure :: insert_particle => flv_structure_insert_particle
<<fks regions: procedures>>=
function flv_structure_insert_particle (flv1, i1, i2, particle) result (flv2)
class(flv_structure_t), intent(in) :: flv1
integer, intent(in) :: i1, i2, particle
type(flv_structure_t) :: flv2
type(flv_structure_t) :: flv_tmp
integer :: n1, n2
n1 = size (flv1%flst)
n2 = n1-1
allocate (flv2%flst(n2))
if (i1 < i2) then
flv_tmp = flv1%remove_particle (i1)
flv_tmp = flv_tmp%remove_particle (i2-1)
else if(i2 < i1) then
flv_tmp = flv1%remove_particle(i2)
flv_tmp = flv_tmp%remove_particle(i1-1)
else
call msg_fatal ("Trying to set ftuple with i1 = i2!")
end if
if (i1 == 1) then
flv2%flst(1) = particle
flv2%flst(2:n2) = flv_tmp%flst(1:n2-1)
else if (i1 == n1 .or. i1 == n2) then
flv2%flst(1:n2-1) = flv_tmp%flst(1:n2-1)
flv2%flst(n2) = particle
else
flv2%flst(1:i1-1) = flv_tmp%flst(1:i1-1)
flv2%flst(i1) = particle
flv2%flst(i1+1:n2) = flv_tmp%flst(i1:n2-1)
end if
end function flv_structure_insert_particle
@ %def flv_structure_insert_particle
@ Returns the number of particles in a flavor array
<<fks regions: flv structure: TBP>>=
procedure :: get_nlegs => flv_structure_get_nlegs
<<fks regions: procedures>>=
function flv_structure_get_nlegs (flv) result(n)
class(flv_structure_t), intent(in) :: flv
integer :: n
n = size (flv%flst)
end function flv_structure_get_nlegs
@ %def flv_structure_get_nlegs
@ Counts the number of occurances of a particle in a
flavor array
<<fks regions: flv structure: TBP>>=
procedure :: count_particle => flv_structure_count_particle
<<fks regions: procedures>>=
function flv_structure_count_particle (flv, part) result (n)
class(flv_structure_t), intent(in) :: flv
integer, intent(in) :: part
integer :: n
n = count (flv%flst == part)
end function flv_structure_count_particle
@ %def flv_structure_count_particle
@ Initializer for flavor structures
<<fks regions: flv structure: TBP>>=
procedure :: init => flv_structure_init
<<fks regions: procedures>>=
subroutine flv_structure_init (flv, aval)
class(flv_structure_t), intent(inout) :: flv
integer, intent(in), dimension(:) :: aval
integer :: n
n = size (aval)
allocate (flv%flst (n))
flv%flst(1:n) = aval(1:n)
flv%nlegs = n
end subroutine flv_structure_init
@ %def flv_structure_init
@
<<fks regions: flv structure: TBP>>=
procedure :: write => flv_structure_write
<<fks regions: procedures>>=
subroutine flv_structure_write (flv, unit)
class(flv_structure_t), intent(in) :: flv
integer, intent(in), optional :: unit
integer :: i, u
u = given_output_unit (unit); if (u < 0) return
write (u, '(A1)',advance = 'no') '['
do i = 1, size(flv%flst)-1
write (u, '(I3,A1)', advance = 'no') flv%flst(i), ','
end do
write (u, '(I3,A1)') flv%flst(i), ']'
end subroutine flv_structure_write
@ %def flv_structure_write
@ Creates the underlying Born flavor structure for a given real flavor
structure if the particle at position \texttt{emitter} is removed
<<fks regions: flv structure: TBP>>=
procedure :: create_uborn => flv_structure_create_uborn
<<fks regions: procedures>>=
function flv_structure_create_uborn (flst_alr, emitter) result(flst_uborn)
class(flv_structure_t), intent(in) :: flst_alr
integer, intent(in) :: emitter
type(flv_structure_t) :: flst_uborn
integer n_alr, n_uborn
n_alr = size(flst_alr%flst)
n_uborn = n_alr-1
allocate (flst_uborn%flst (n_uborn))
if (emitter > 2) then
if (flst_alr%flst(n_alr) == 21) then
!!! Emitted particle is a gluon => just remove it
flst_uborn = flst_alr%remove_particle(n_alr)
!!! Emission type is a gluon splitting into two quars
else if (is_quark (abs(flst_alr%flst(n_alr))) .and. &
is_quark (abs(flst_alr%flst(n_alr-1))) .and. &
flst_alr%flst(n_alr) + flst_alr%flst(n_alr-1) == 0) then
flst_uborn = flst_alr%insert_particle(n_alr-1,n_alr,21)
end if
else
if (flst_alr%flst(n_alr) == 21) then
flst_uborn = flst_alr%remove_particle(n_alr)
else if (is_quark (abs(flst_alr%flst(n_alr))) .and. &
is_gluon (abs(flst_alr%flst(emitter)))) then
flst_uborn = &
flst_alr%insert_particle (emitter,n_alr,-flst_alr%flst(n_alr))
else if (is_quark (abs(flst_alr%flst(n_alr))) .and. &
is_quark (abs(flst_alr%flst(emitter))) .and. &
flst_alr%flst(n_alr) == flst_alr%flst(emitter)) then
flst_uborn = flst_alr%insert_particle(emitter,n_alr,21)
end if
end if
end function flv_structure_create_uborn
@ %def flv_structure_create_uborn
@
<<fks regions: flv structure: TBP>>=
procedure :: evaluate => flv_structure_evaluate
<<fks regions: procedures>>=
subroutine flv_structure_evaluate (flv, n, model)
class(flv_structure_t), intent(inout) :: flv
integer, intent(in) :: n
type(model_data_t), intent(in) :: model
integer :: i
type(flavor_t) :: flavor
allocate (flv%massive (n), flv%colored(n))
do i = 1, n
call flavor_init (flavor, flv%flst(i), model)
flv%massive(i) = flavor_get_mass (flavor) > 0
flv%colored(i) = is_quark (abs(flv%flst(i))) .or. &
is_gluon (flv%flst(i))
end do
end subroutine flv_structure_evaluate
@ %def flv_structure_evaluate
@
<<fks regions: public>>=
public :: singular_region_t
<<fks regions: types>>=
type :: singular_region_t
integer :: alr
type(flv_structure_t) :: flst_real
type(flv_structure_t) :: flst_uborn
integer :: mult
integer :: emitter
integer :: nregions
integer :: real_index
type(ftuple_t), dimension(:), allocatable :: flst_allreg
integer :: uborn_index
logical :: double_fsr
logical :: soft_divergence
logical :: coll_divergence
contains
<<fks regions: singular region: TBP>>
end type singular_region_t
@ %def singular_region_t
@ In case of a $g \rightarrow gg$ or $g \rightarrow qq$ splitting, the factor
\begin{equation*}
\frac{2E_{\rm{em}}}{E_{\rm{em}} + E_{\rm{rad}}}
\end{equation*}
is multiplied to the real matrix element. This way, the symmetry of the splitting is used
and only one singular region has to be taken into account. However, the factor ensures that
there is only a soft singularity if the radiated parton becomes soft.
<<fks regions: singular region: TBP>>=
procedure :: set_splitting_info => singular_region_set_splitting_info
<<fks regions: procedures>>=
subroutine singular_region_set_splitting_info (region)
class(singular_region_t), intent(inout) :: region
integer :: i1, i2
integer :: reg
region%double_fsr = .false.
associate (ftuple => region%flst_allreg)
do reg = 1, region%nregions
call ftuple(reg)%get (i1, i2)
if (i1 /= region%emitter) then
cycle
else
region%soft_divergence = &
ftuple(reg)%splitting_type /= G_TO_QQ
region%coll_divergence = &
.not. region%flst_real%massive(i1)
if (ftuple(reg)%splitting_type > 1) then
region%double_fsr = .true.
exit
else if (ftuple(reg)%splitting_type == UNDEFINED_SPLITTING) then
call msg_fatal ("All splittings should be defined!")
end if
end if
end do
end associate
end subroutine singular_region_set_splitting_info
@ %def singular_region_set_splitting_info
@
<<fks regions: singular region: TBP>>=
procedure :: double_fsr_factor => singular_region_double_fsr_factor
<<fks regions: procedures>>=
function singular_region_double_fsr_factor (region, p) result (val)
class(singular_region_t), intent(in) :: region
type(vector4_t), intent(in), dimension(:) :: p
real(default) :: val
real(default) :: E_rad, E_em
if (region%double_fsr) then
E_em = energy (p(region%emitter))
E_rad = energy (p(region%flst_real%get_nlegs()))
val = 2*E_em/(E_em + E_rad)
else
val = 1._default
end if
end function singular_region_double_fsr_factor
@ %def singular_region_double_fsr_factor
@
<<fks regions: singular region: TBP>>=
procedure :: has_soft_divergence => singular_region_has_soft_divergence
<<fks regions: procedures>>=
function singular_region_has_soft_divergence (region) result (div)
class(singular_region_t), intent(in) :: region
logical :: div
div = region%soft_divergence
end function singular_region_has_soft_divergence
@ %def singular_region_has_soft_divergence
@
<<fks regions: singular region: TBP>>=
procedure :: has_collinear_divergence => &
singular_region_has_collinear_divergence
<<fks regions: procedures>>=
function singular_region_has_collinear_divergence (region) result (div)
class(singular_region_t), intent(in) :: region
logical :: div
div = region%coll_divergence
end function singular_region_has_collinear_divergence
@ %def singular_region_has_collinear_divergence
@
<<fks regions: types>>=
type, abstract :: fks_mapping_t
real(default) :: sumdij
real(default) :: sumdij_soft
contains
<<fks regions: fks mapping: TBP>>
end type fks_mapping_t
@ %def fks_mapping_t
@
<<fks regions: public>>=
public :: fks_mapping_default_t
<<fks regions: types>>=
type, extends (fks_mapping_t) :: fks_mapping_default_t
real(default) :: exp_1, exp_2
contains
<<fks regions: fks mapping 1: TBP>>
end type fks_mapping_default_t
@ %def fks_mapping_default_t
@
<<fks regions: interfaces>>=
interface operator(==)
module procedure flv_structure_equivalent
end interface
@ %def operator_equiv
<<fks regions: public>>=
public :: region_data_t
<<fks regions: types>>=
type :: region_data_t
type(singular_region_t), dimension(:), allocatable :: regions
type(flv_structure_t), dimension(:), allocatable :: flv_born
type(flv_structure_t), dimension(:), allocatable :: flv_real
integer, dimension(:), allocatable :: emitters
integer :: n_regions
integer :: n_emitters
integer :: n_flv_born
integer :: n_flv_real
integer :: nlegs_born
integer :: nlegs_real
type(flavor_t) :: flv_extra
class(fks_mapping_t), allocatable :: fks_mapping
contains
<<fks regions: reg data: TBP>>
end type region_data_t
@ %def region_data_t
@
<<fks regions: reg data: TBP>>=
procedure :: init => region_data_init
<<fks regions: procedures>>=
subroutine region_data_init (reg_data, model, flavor_born, &
flavor_real, mapping_type)
class(region_data_t), intent(inout) :: reg_data
type(model_data_t), intent(in) :: model
integer, intent(inout), dimension(:,:), allocatable :: &
flavor_born, flavor_real
integer, intent(in) :: mapping_type
integer, dimension(:), allocatable :: current_flavor
type(ftuple_list_t), dimension(:), allocatable :: ftuples
integer, dimension(:), allocatable :: emitter
type(flv_structure_t), dimension(:), allocatable :: flst_alr
integer :: i
reg_data%n_flv_born = size(flavor_born(1,:))
reg_data%n_flv_real = size(flavor_real(1,:))
reg_data%nlegs_born = size(flavor_born(:,1))
reg_data%nlegs_real = reg_data%nlegs_born + 1
allocate (reg_data%flv_born (reg_data%n_flv_born))
allocate (reg_data%flv_real (reg_data%n_flv_real))
allocate (current_flavor (reg_data%n_flv_born))
do i = 1, reg_data%n_flv_born
current_flavor = flavor_born(:,i)
call reg_data%flv_born(i)%init (current_flavor)
end do
deallocate (current_flavor)
allocate (current_flavor (reg_data%n_flv_real))
do i = 1, reg_data%n_flv_real
current_flavor = flavor_real(:,i)
call reg_data%flv_real(i)%init (current_flavor)
end do
select case (mapping_type)
case (1)
allocate (fks_mapping_default_t :: reg_data%fks_mapping)
case default
call msg_fatal ("Init region_data: FKS mapping not implemented!")
end select
call flavor_init (reg_data%flv_extra, &
reg_data%flv_real(1)%flst(reg_data%nlegs_real), &
model)
call reg_data%find_regions (model, ftuples, emitter, flst_alr)
call reg_data%init_regions (ftuples, emitter, flst_alr)
call reg_data%evaluate_flavors (model)
call reg_data%set_splitting_info ()
call reg_data%find_emitters ()
call reg_data%write ()
end subroutine region_data_init
@ %def region_data_init
@
<<fks regions: reg data: TBP>>=
procedure :: evaluate_flavors => region_data_evaluate_flavors
<<fks regions: procedures>>=
subroutine region_data_evaluate_flavors (reg_data, model)
class(region_data_t), intent(inout) :: reg_data
type(model_data_t), intent(in) :: model
integer :: i
do i = 1, reg_data%n_regions
associate (region => reg_data%regions(i))
call region%flst_uborn%evaluate (reg_data%nlegs_born, model)
call region%flst_real%evaluate (reg_data%nlegs_real, model)
end associate
end do
do i = 1, reg_data%n_flv_born
call reg_data%flv_born(i)%evaluate (reg_data%nlegs_born, model)
end do
do i = 1, reg_data%n_flv_real
call reg_data%flv_real(i)%evaluate (reg_data%nlegs_real, model)
end do
end subroutine region_data_evaluate_flavors
@ %def region_data_evaluate_flavors
@ Creates a list containing the emitter of each singular region.
<<fks regions: reg data: TBP>>=
procedure :: get_emitters => region_data_get_emitters
<<fks regions: procedures>>=
function region_data_get_emitters (reg_data) result(emitters)
class(region_data_t), intent(inout) :: reg_data
integer, dimension(:), allocatable :: emitters
integer :: i
allocate (emitters (reg_data%n_regions))
do i = 1, reg_data%n_regions
emitters(i) = reg_data%regions(i)%emitter
end do
end function region_data_get_emitters
@ %def region_data_get_emitters
@ Returns $S_i = \frac{1}{\mathcal{D}d_i}$ or $S_{ij} =
\frac{1}{\mathcal{D}d_{ij}}$ for one particular singular region. At
this point, the flavor array should be rearranged in such a way that
the emitted particle is at the last position of
the flavor structure list.
<<fks regions: reg data: TBP>>=
procedure :: get_svalue => region_data_get_svalue
<<fks regions: procedures>>=
function region_data_get_svalue (reg_data, p, alr, emitter) result (sval)
class(region_data_t), intent(inout) :: reg_data
type(vector4_t), intent(in), dimension(:) :: p
integer, intent(in) :: alr, emitter
real(default) :: sval
associate (map => reg_data%fks_mapping)
map%sumdij = map%compute_sumdij (reg_data%regions(alr), p)
sval = map%svalue (p, emitter, reg_data%nlegs_real)
end associate
end function region_data_get_svalue
@ %def region_data_get_svalue
@ The same as above, but for the soft limit.
<<fks regions: reg data: TBP>>=
procedure :: get_svalue_soft => region_data_get_svalue_soft
<<fks regions: procedures>>=
function region_data_get_svalue_soft &
(reg_data, p, p_soft, alr, emitter) result (sval)
class(region_data_t), intent(inout) :: reg_data
type(vector4_t), intent(in), dimension(:) :: p
type(vector4_t), intent(in) :: p_soft
integer, intent(in) :: alr, emitter
real(default) :: sval
associate (map => reg_data%fks_mapping)
map%sumdij_soft = &
map%compute_sumdij_soft (reg_data%regions(alr), p, p_soft)
sval = map%svalue_soft (p, p_soft, emitter)
end associate
end function region_data_get_svalue_soft
@ %def region_data_get_svalue_soft
@ This subroutine starts with a specification of $N$- and
$N+1$-particle configurations, [[flst_born]] and [[flst_real]], saved
in [[reg_data]]. From these, it creates a list of fundamental tuples,
a list of emitters and a list containing the $N+1$-particle
configuration, rearranged in such a way that the emitter-radiation
pair is last ([[flst_alr]]). For the $e^+ \, e^- \, \rightarrow u \,
\bar{u} \, g$- example, the generated objects are shown in table
\ref{table:ftuples and flavors}. Note that at this point, [[flst_alr]]
is arranged in such a way that the emitter can only be equal to
$n_{legs}-1$ for final-state radiation or 0, 1, or 2 for initial-state
radiation. Further, it occurs that regions can be equivalent. For
example in table \ref{table:ftuples and flavors} the regions
correpsonding to \texttt{alr} = 1 and \texttt{alr} = 3 as well as
\texttt{alr} = 2 and \texttt{alr} = 4 describe the same physics and
are therefore equivalent.
@
<<fks regions: reg data: TBP>>=
procedure :: find_regions => region_data_find_regions
<<fks regions: procedures>>=
subroutine region_data_find_regions &
(reg_data, model, ftuples, emitter, flst_alr)
class(region_data_t), intent(in) :: reg_data
type(model_data_t), intent(in) :: model
type(ftuple_list_t), intent(out), dimension(:), allocatable :: ftuples
integer, intent(out), dimension(:), allocatable :: emitter
type(flv_structure_t), intent(out), dimension(:), allocatable :: flst_alr
type(ftuple_t) :: current_ftuple
integer, dimension(:), allocatable :: emitter_tmp
type(flv_structure_t), dimension(:), allocatable :: flst_alr_tmp
integer :: nreg, nborn, nreal
integer :: nlegreal
integer :: n_gluon_real, n_gluon_born
integer, parameter :: maxnregions = 100
integer :: i, j, k, l, n
associate (flv_born => reg_data%flv_born)
associate (flv_real => reg_data%flv_real)
nborn = size (flv_born)
nreal = size (flv_real)
nlegreal = size (flv_real(1)%flst)
allocate (ftuples (nreal))
allocate (emitter_tmp (maxnregions))
allocate (flst_alr_tmp (maxnregions))
n = 0
ITERATE_REAL_FLAVOR: do l = 1, nreal
call ftuples(l)%init
<<fks: check final state emissions>>
<<fks: check initial state emissions>>
end do ITERATE_REAL_FLAVOR
nreg = n
end associate
end associate
allocate (flst_alr (nreg))
allocate (emitter (nreg))
flst_alr(1:nreg) = flst_alr_tmp(1:nreg)
emitter(1:nreg) = emitter_tmp(1:nreg)
end subroutine region_data_find_regions
@ %def region_data_find_regions
@ First check final state emissions
<<fks: check final state emissions>>=
do i = 3, nlegreal
do j = i+1, nlegreal
do k = 1, nborn
if (flv_real(l)%valid_pair(i,j, flv_born(k), model) &
.or. flv_real(l)%valid_pair(j,i,flv_born(k), model)) then
n = n+1
n_gluon_real = flv_real(l)%count_particle (GLUON)
n_gluon_born = flv_born(k)%count_particle (GLUON)
if (n_gluon_born - n_gluon_real < 0) then
if(flv_real(l)%valid_pair(i,j, flv_born(k), model)) then
flst_alr_tmp(n) = create_alr (flv_real(l),i,j)
else
flst_alr_tmp(n) = create_alr (flv_real(l),j,i)
end if
else
flst_alr_tmp(n) = flv_real(l)
end if
call current_ftuple%set (i,j)
call current_ftuple%determine_splitting_type (flv_real(l), i, j)
call ftuples(l)%append (current_ftuple)
emitter_tmp(n) = nlegreal - 1
exit
end if
end do
end do
@ Check initial-state emissions. It suffices to just check the
final-state of the first and the initial-state of the second array.
<<fks: check initial state emissions>>=
do k = 1, nborn
if (flv_real(l)%valid_pair(1,i, flv_born(k), model) &
.and. flv_real(l)%valid_pair(2,i, flv_born(k), model)) then
n = n + 1
call current_ftuple%set (0,i)
call ftuples(l)%append (current_ftuple)
emitter_tmp(n) = 0
flst_alr_tmp(n) = create_alr (flv_real(l),0,i)
exit
else if (flv_real(l)%valid_pair(1,i, flv_born(k), model) &
.and. .not. &
flv_real(l)%valid_pair(2,i, flv_born(k), model)) then
n = n+1
call current_ftuple%set (1,i)
call ftuples(l)%append (current_ftuple)
emitter_tmp(n) = 1
flst_alr_tmp(n) = create_alr (flv_real(l),1,i)
exit
else if (flv_real(l)%valid_pair(2,i, flv_born(k), model) &
.and. .not. &
flv_real(l)%valid_pair(1,i, flv_born(k), model)) then
n = n+1
call current_ftuple%set(2,i)
call ftuples(l)%append (current_ftuple)
emitter_tmp(n) = 2
flst_alr_tmp(n) = create_alr (flv_real(l),2,i)
exit
end if
end do
end do
@ Creates singular regions according to table \ref{table:singular
regions}. It scans all regions in table \ref{table:ftuples and
flavors} and records the real flavor structures. If they are
equivalent, the flavor structure is not recorded, but the multiplicity
of the present one is increased.
<<fks regions: reg data: TBP>>=
procedure :: init_regions => region_data_init_singular_regions
<<fks regions: procedures>>=
subroutine region_data_init_singular_regions &
(reg_data, ftuples, emitter, flst_alr)
class(region_data_t), intent(inout) :: reg_data
type(ftuple_list_t), intent(inout), dimension(:), allocatable :: ftuples
type(ftuple_list_t) :: current_region
integer, intent(in), dimension(:), allocatable :: emitter
type(flv_structure_t), intent(in), dimension(:), allocatable :: flst_alr
type(flv_structure_t), dimension(:), allocatable :: flst_uborn, flst_alr2
integer, dimension(:), allocatable :: mult
integer, dimension(:), allocatable :: flst_emitter
integer :: nregions, maxregions
integer, dimension(:,:), allocatable :: perm_list
integer, dimension(:), allocatable :: index
integer :: i, j, k, l
integer :: nlegs
logical :: equiv
integer :: nreg, i1, i2
integer :: i_first, j_first
integer, dimension(:), allocatable :: &
region_to_ftuple, ftuple_limits, k_index
type(flv_structure_t) :: flst_save
maxregions = size (emitter)
nlegs = size(flst_alr(1)%flst)
allocate (flst_uborn (maxregions))
allocate (flst_alr2 (maxregions))
allocate (mult (maxregions))
allocate (flst_emitter (maxregions))
allocate (index (maxregions))
allocate (region_to_ftuple (maxregions))
allocate (ftuple_limits (size (ftuples)))
allocate (k_index (maxregions))
mult = 0
do i = 1, size(ftuples)
ftuple_limits(i) = ftuples(i)%get_n_tuples ()
end do
if (.not. (sum (ftuple_limits) == maxregions)) &
call msg_fatal ("Too many regions!")
k = 1
do j =1, size(ftuples)
do i = 1, ftuple_limits(j)
region_to_ftuple(k) = i
k = k + 1
end do
end do
i_first = 1
j_first = 1
j = 1
SCAN_REGIONS: do l = 1, size(ftuples)
SCAN_FTUPLES: do i = i_first, i_first + ftuple_limits (l) -1
equiv = .false.
if (i==i_first) then
flst_alr2(j)%flst = flst_alr(i)%flst
mult(j) = mult(j) + 1
flst_uborn(j) = flst_alr(i)%create_uborn (emitter(i))
flst_emitter(j) = emitter(i)
index (j) = region_to_index(ftuples, i)
k_index (j) = region_to_ftuple(i)
j = j+1
else
!!! Check for equivalent flavor structures
do k =j_first ,j-1
if (emitter(i) == emitter(k) .and. emitter(i) > 2) then
if (flst_alr(i) == flst_alr2(k) .and. &
flst_alr(i)%flst(nlegs-1) == flst_alr2(k)%flst(nlegs-1) &
.and. flst_alr(i)%flst(nlegs) == flst_alr2(k)%flst(nlegs)) then
mult(k) = mult(k) + 1
equiv = .true.
call ftuples (region_to_index(ftuples, i))%set_equiv &
(k_index(k), region_to_ftuple(i))
exit
end if
else if (emitter(i) == emitter(k) .and. emitter(i) <= 2) then
if (flst_alr(i) == flst_alr2(k)) then
mult(k) = mult(k) + 1
equiv = .true.
call ftuples (region_to_index(ftuples,i))%set_equiv &
(k_index(k), region_to_ftuple(i))
exit
end if
end if
end do
if (.not.equiv) then
flst_alr2(j)%flst = flst_alr(i)%flst
mult(j) = mult(j) + 1
flst_uborn(j) = flst_alr(i)%create_uborn (emitter(i))
flst_emitter(j) = emitter(i)
index (j) = region_to_index (ftuples, i)
k_index (j) = region_to_ftuple(i)
j = j+1
end if
end if
end do SCAN_FTUPLES
i_first = i_first + ftuple_limits(l)
j_first = j_first + j - 1
end do SCAN_REGIONS
nregions = j-1
allocate (reg_data%regions (nregions))
do j = 1, nregions
do i = 1, reg_data%n_flv_born
if (reg_data%flv_born (i) == flst_uborn (j)) then
if (allocated (perm_list)) deallocate (perm_list)
call fks_permute_born &
(reg_data%flv_born (i), flst_uborn (j), perm_list)
call fks_apply_perm (flst_alr2(j), flst_emitter(j), perm_list)
end if
end do
end do
!!! Check if new emitters require a rearrangement of ftuples
do i = 1, nregions
reg_data%regions(i)%alr = i
reg_data%regions(i)%flst_real = flst_alr2(i)
reg_data%regions(i)%mult = mult(i)
reg_data%regions(i)%flst_uborn = flst_uborn(i)
reg_data%regions(i)%emitter = flst_emitter(i)
nreg = ftuples (index(i))%get_n_tuples ()
reg_data%regions(i)%nregions = nreg
allocate (reg_data%regions(i)%flst_allreg (nreg))
do j = 1, nreg
current_region = ftuples (index(i))%get_entry (j)
if (.not. associated (current_region%equiv)) then
call current_region%ftuple%get (i1, i2)
if (i2 /= nlegs) &
call current_region%ftuple%set (i1, nlegs)
end if
reg_data%regions(i)%flst_allreg (j) = current_region%ftuple
end do
end do
!!! Find underlying Born index
do j = 1, nregions
do i = 1, reg_data%n_flv_born
if (reg_data%flv_born (i) == reg_data%regions(j)%flst_uborn) then
reg_data%regions(j)%uborn_index = i
exit
end if
end do
end do
k = 1
associate (regions => reg_data%regions)
do i = 1, nregions
if (i==1) then
regions(i)%real_index = 1
flst_save = flst_alr2(1)
cycle
end if
if (flst_alr2(i) == flst_save) then
regions(i)%real_index = k
else
k = k+1
regions(i)%real_index = k
flst_save = flst_alr2(i)
end if
end do
end associate
reg_data%n_regions = size (reg_data%regions)
end subroutine region_data_init_singular_regions
@ %def region_data_init_singular_regions
@ Create an array containing all emitters of a singular region.
<<fks regions: reg data: TBP>>=
procedure :: find_emitters => region_data_find_emitters
<<fks regions: procedures>>=
subroutine region_data_find_emitters (reg_data)
class(region_data_t), intent(inout) :: reg_data
integer :: i, j, n
integer :: em
integer, dimension(10) :: em_count
em_count = 0
n = 0
@ %def region_data_find_emitters
@ Count the number of different emitters.
<<fks regions: procedures>>=
do i = 1, reg_data%n_regions
em = reg_data%regions(i)%emitter
if (.not. any (em_count == em)) then
n = n+1
em_count(i) = em
end if
end do
if (n < 1) call msg_fatal ("region_data_find_emitters: No emitters found")
reg_data%n_emitters = n
allocate (reg_data%emitters (reg_data%n_emitters))
reg_data%emitters = 0
j = 1
do i = 1, size(reg_data%regions)
em = reg_data%regions(i)%emitter
if (.not. any (reg_data%emitters == em)) then
reg_data%emitters(j) = em
j = j+1
end if
end do
end subroutine region_data_find_emitters
@ %def region_data_find_emitters
@
<<fks regions: reg data: TBP>>=
procedure :: set_splitting_info => region_data_set_splitting_info
<<fks regions: procedures>>=
subroutine region_data_set_splitting_info (reg_data)
class(region_data_t), intent(inout) :: reg_data
integer :: reg
do reg = 1, reg_data%n_regions
call reg_data%regions(reg)%set_splitting_info ()
end do
end subroutine region_data_set_splitting_info
@ %def region_data_set_splitting_info
@
<<fks regions: reg data: TBP>>=
procedure :: write_regions => region_data_write_regions
<<fks regions: procedures>>=
subroutine region_data_write_regions (reg_data)
class(region_data_t), intent(inout) :: reg_data
integer :: n, i, j
n = size(reg_data%regions)
associate (regions => reg_data%regions)
do i = 1, n
print *, i, '//', regions(i)%flst_real%flst, '//', &
regions(i)%mult ,'//', regions(i)%flst_uborn%flst , &
'//', regions(i)%emitter
do j = 1, size (regions(i)%flst_allreg)
call regions(i)%flst_allreg(j)%write
end do
end do
end associate
end subroutine region_data_write_regions
@ %def region_data_write_regions
@ Creates a table with information about all singular regions and
writes it to a file.
@ Returns the index of the real flavor structure an ftuple belongs to.
<<fks regions: reg data: TBP>>=
procedure :: write => region_data_write
<<fks regions: procedures>>=
subroutine region_data_write (reg_data, proc)
class(region_data_t), intent(inout) :: reg_data
type(string_t), intent(inout), optional :: proc
integer :: u, i, j
integer :: nreal, nborn
integer :: i1, i2, nreg
integer :: maxnregions, nreg_diff
integer :: nleft, nright
type(singular_region_t) :: region
type(string_t) :: flst_title, ftuple_title
character(len=7) :: flst_format = "(I3,A1)"
character(len=16) :: ireg_format = "(A1,I3,A1,I3,A2)"
character(len=7) :: ireg_space_format = "(7X,A1)"
u = free_unit ()
open (u, file="region_data.log", action = "write", status="replace")
maxnregions = 1
do j = 1, reg_data%n_regions
if (size (reg_data%regions(j)%flst_allreg) > maxnregions) &
maxnregions = reg_data%regions(j)%nregions
end do
flst_title = '(A' // flst_title_format(reg_data%nlegs_real) // ')'
ftuple_title = '(A' // ftuple_title_format() // ')'
write (u,*) 'Total number of regions: ', size(reg_data%regions)
write (u, '(A3)', advance = 'no') 'alr'
call write_separator ()
write (u, char (flst_title), advance = 'no') 'flst_real'
call write_separator ()
write (u, '(A3)', advance = 'no') 'em'
call write_separator ()
write (u, '(A3)', advance = 'no') 'mult'
call write_separator ()
write (u, '(A4)', advance = 'no') 'nreg'
call write_separator ()
write (u, char (ftuple_title), advance = 'no') 'ftuples'
call write_separator ()
flst_title = '(A' // flst_title_format(reg_data%nlegs_born) // ')'
write (u, char (flst_title), advance = 'no') 'flst_born'
call write_separator ()
write (u, '(A7)') 'i_uborn'
do j = 1, reg_data%n_regions
region = reg_data%regions(j)
nreal = size (region%flst_real%flst)
nborn = size (region%flst_uborn%flst)
write (u, '(I3)', advance = 'no') j
call write_separator ()
write (u, '(A1)', advance = 'no') '['
do i = 1, nreal-1
write (u, flst_format, advance = 'no') region%flst_real%flst(i), ','
end do
write (u, flst_format, advance = 'no') region%flst_real%flst(nreal), ']'
call write_separator ()
write (u, '(I3)', advance = 'no') region%emitter
call write_separator ()
write (u, '(I3)', advance = 'no') region%mult
call write_separator ()
write (u, '(I4)', advance = 'no') region%nregions
call write_separator ()
!!! write ftuples
nreg = region%nregions
if (nreg == maxnregions) then
nleft = 0
nright = 0
else
nreg_diff = maxnregions - nreg
nleft = nreg_diff/2
if (mod(nreg_diff,2) == 0) then
nright = nleft
else
nright = nleft + 1
end if
end if
if (nleft > 0) then
do i=1,nleft
write(u,ireg_space_format, advance='no') ' '
end do
end if
write (u,'(A1)', advance = 'no') '{'
if (nreg > 1) then
do i=1,nreg-1
call region%flst_allreg(i)%get (i1, i2)
write(u,ireg_format,advance = 'no') '(', i1, ',', i2, '),'
end do
end if
call region%flst_allreg(nreg)%get (i1, i2)
write (u,ireg_format,advance = 'no') '(', i1, ',', i2, ')}'
if (nright > 0) then
do i=1,nright
write(u,ireg_space_format, advance='no') ' '
end do
end if
!!! end write ftuples
call write_separator ()
write (u,'(A1)',advance = 'no') '['
do i=1,nborn-1
write(u,flst_format,advance = 'no') region%flst_uborn%flst(i), ','
end do
write (u,flst_format, advance = 'no') region%flst_uborn%flst(nborn), ']'
call write_separator ()
write (u, '(I7)', advance = 'no') region%uborn_index
write(u,*) ''
end do
close (u)
contains
function flst_title_format (n) result (frmt)
integer, intent(in) :: n
type(string_t) :: frmt
character(len=2) :: frmt_char
write (frmt_char, '(I2)') 4*n+1
frmt = var_str (frmt_char)
end function flst_title_format
function ftuple_title_format () result (frmt)
type(string_t) :: frmt
character(len=2) :: frmt_char
write (frmt_char, '(I2)') 10*maxnregions+1
frmt = var_str (frmt_char)
end function ftuple_title_format
subroutine write_separator ()
character(len=10) :: sep_format = "(1X,A2,1X)"
write (u, sep_format, advance = 'no') '||'
end subroutine write_separator
end subroutine region_data_write
@ %def region_data_write
@ Returns the index of the real flavor structure an ftuple belogs to.
<<fks regions: procedures>>=
function region_to_index (list, i) result(index)
type(ftuple_list_t), intent(inout), dimension(:), allocatable :: list
integer, intent(in) :: i
integer :: index
integer :: nlist
integer :: j
integer, dimension(:), allocatable :: nreg
nlist = size(list)
allocate (nreg (nlist))
do j = 1, nlist
if (j == 1) then
nreg(j) = list(j)%get_n_tuples ()
else
nreg(j) = nreg(j-1) + list(j)%get_n_tuples ()
end if
end do
do j = 1, nlist
if (j == 1) then
if (i <= nreg(j)) then
index = j
exit
end if
else
if (i > nreg(j-1) .and. i <= nreg(j)) then
index = j
exit
end if
end if
end do
end function region_to_index
@ %def region_to_index
@ Rearrange the flavor array in such a way that the emitted particle
is last and the emitter is second last. [[i1]] is the index of the
emitter, [[i2]] is the index of the emitted particle. Only works for
final-state emitters.
<<fks regions: procedures>>=
function create_alr (flv1,i1,i2) result(flv2)
type(flv_structure_t), intent(in) :: flv1
integer, intent(in) :: i1, i2
type(flv_structure_t) :: flv2
integer :: n, i, j
n = size (flv1%flst)
allocate (flv2%flst (n))
if (i1 > 2) then
flv2%flst(1:2) = flv1%flst(1:2)
flv2%flst(n-1) = flv1%flst(i1)
flv2%flst(n) = flv1%flst(i2)
@ Order remaining particles according to their original position
<<fks regions: procedures>>=
j = 3
do i = 3,n
if (i /= i1 .and. i /= i2) then
flv2%flst(j) = flv1%flst(i)
j = j+1
end if
end do
else
call msg_fatal ("Create alr: Only works for final-state emissions!")
end if
end function create_alr
@ %def create_alr
@ Explain
<<fks regions: procedures>>=
subroutine fks_permute_born (flv_in, flv_out, perm_list)
type(flv_structure_t), intent(in) :: flv_in
type(flv_structure_t), intent(inout) :: flv_out
integer, intent(out), dimension(:,:), allocatable :: perm_list
integer, dimension(:,:), allocatable :: perm_list_tmp
integer :: n_perms, n_perms_max
integer :: nlegs
integer :: flv1, flv2, tmp
integer :: i, j, j_min
n_perms_max = 100
!!! actually (n-1)!, but there seems to be no intrinsic function
!!! of this type in fortran
if (allocated (perm_list_tmp)) deallocate (perm_list_tmp)
allocate (perm_list_tmp (n_perms_max,2))
n_perms = 0
j_min = 3
nlegs = size (flv_in%flst)
do i = 3, nlegs
flv1 = flv_in%flst(i)
do j = j_min, nlegs
flv2 = flv_out%flst(j)
if (flv1 == flv2 .and. i /= j) then
n_perms = n_perms + 1
tmp = flv_out%flst(i)
flv_out%flst(i) = flv2
flv_out%flst(j) = tmp
perm_list_tmp (n_perms, 1) = j
perm_list_tmp (n_perms, 2) = i
j_min = j_min + 1
exit
end if
end do
end do
allocate (perm_list (n_perms, 2))
perm_list (1:n_perms, :) = perm_list_tmp (1:n_perms, :)
end subroutine fks_permute_born
@ %def fks_permute_born
@ Explain
<<fks regions: procedures>>=
subroutine fks_apply_perm (flv, emitter, perm_list)
type(flv_structure_t), intent(inout) :: flv
integer, intent(inout) :: emitter
integer, intent(in), dimension(:,:), allocatable :: perm_list
integer :: i
integer :: i1, i2
integer :: tmp
do i = 1, size (perm_list (:,1))
i1 = perm_list (i,1)
i2 = perm_list (i,2)
tmp = flv%flst (i1)
flv%flst (i1) = flv%flst (i2)
flv%flst (i2) = tmp
if (i1 == emitter) emitter = i2
end do
end subroutine fks_apply_perm
@ %def fks_apply_perm
@ Translates the tree code of the emitter branch into the position of
the emitter in the flavor structure array.
<<fks regions: public>>=
public :: fks_tree_to_position
<<fks regions: procedures>>=
function fks_tree_to_position (k, n_tot) result(pos)
integer, intent(in) :: k, n_tot
integer :: pos
integer :: k_tot
k_tot = 2**(n_tot - 1)
!!! Inital-state particles
if (k == k_tot) then
pos = 1
else if (k == k_tot/2) then
pos = 2
!!! Final-state particles
else
! pos = 3 + nint(log(k)/log(2))
pos = 3 + dual_log (k)
end if
contains
recursive function dual_log (x) result (ld)
integer, intent(in) :: x
integer :: ld
if (x == 1) then
ld = 0
else
ld = 1 + dual_log (x/2)
end if
end function dual_log
end function fks_tree_to_position
@ %def fks_tree_to_position
@
<<fks regions: fks mapping: TBP>>=
procedure (fks_mapping_dij), deferred :: dij
<<fks regions: interfaces>>=
abstract interface
function fks_mapping_dij (map, p, i, j) result (d)
import
class(fks_mapping_t), intent(in) :: map
type(vector4_t), intent(in), dimension(:) :: p
integer, intent(in) :: i, j
real(default) :: d
end function fks_mapping_dij
end interface
@ %def fks_mapping_dij
@
<<fks regions: fks mapping: TBP>>=
procedure (fks_mapping_compute_sumdij), deferred :: compute_sumdij
<<fks regions: interfaces>>=
abstract interface
function fks_mapping_compute_sumdij (map, sregion, p) result (d)
import
class(fks_mapping_t), intent(in) :: map
type(singular_region_t), intent(inout) :: sregion
type(vector4_t), intent(in), dimension(:) :: p
real(default) :: d
end function fks_mapping_compute_sumdij
end interface
@ %def fks_mapping_compute_sumdij
@
<<fks regions: fks mapping: TBP>>=
procedure (fks_mapping_svalue), deferred :: svalue
<<fks regions: interfaces>>=
abstract interface
function fks_mapping_svalue (map, p, i, j) result (value)
import
class(fks_mapping_t), intent(in) :: map
type(vector4_t), intent(in), dimension(:) :: p
integer, intent(in) :: i, j
real(default) :: value
end function fks_mapping_svalue
end interface
@ %def fks_mapping_svalue
<<fks regions: fks mapping: TBP>>=
procedure (fks_mapping_dij_soft), deferred :: dij_soft
<<fks regions: interfaces>>=
abstract interface
function fks_mapping_dij_soft (map, p_born, p_soft, em) result (d)
import
class(fks_mapping_t), intent(in) :: map
type(vector4_t), intent(in), dimension(:) :: p_born
type(vector4_t), intent(in) :: p_soft
integer, intent(in) :: em
real(default) :: d
end function fks_mapping_dij_soft
end interface
@ %def fks_mapping_dij_soft
@
<<fks regions: fks mapping: TBP>>=
procedure (fks_mapping_compute_sumdij_soft), deferred :: compute_sumdij_soft
<<fks regions: interfaces>>=
abstract interface
function fks_mapping_compute_sumdij_soft (map, sregion, p_born, p_soft) result (d)
import
class(fks_mapping_t), intent(in) :: map
type(singular_region_t), intent(inout) :: sregion
type(vector4_t), intent(in), dimension(:) :: p_born
type(vector4_t), intent(in) :: p_soft
real(default) :: d
end function
end interface
@ %def fks_mapping_compute_sumdij_soft
@
<<fks regions: fks mapping: TBP>>=
procedure (fks_mapping_svalue_soft), deferred :: svalue_soft
<<fks regions: interfaces>>=
abstract interface
function fks_mapping_svalue_soft (map, p_born, p_soft, em) result (value)
import
class(fks_mapping_t), intent(in) :: map
type(vector4_t), intent(in), dimension(:) :: p_born
type(vector4_t), intent(in) :: p_soft
integer, intent(in) :: em
real(default) :: value
end function fks_mapping_svalue_soft
end interface
@ %def fks_mapping_svalue_soft
@
<<fks regions: fks mapping 1: TBP>>=
procedure :: set_parameter => fks_mapping_default_set_parameter
<<fks regions: procedures>>=
subroutine fks_mapping_default_set_parameter (map, dij_exp1, dij_exp2)
class(fks_mapping_default_t), intent(inout) :: map
real(default), intent(in) :: dij_exp1, dij_exp2
map%exp_1 = dij_exp1
map%exp_2 = dij_exp2
end subroutine fks_mapping_default_set_parameter
@ %def fks_mapping_default_set_parameter
@
<<fks regions: fks mapping 1: TBP>>=
procedure :: dij => fks_mapping_default_dij
<<fks regions: procedures>>=
function fks_mapping_default_dij (map, p, i, j) result (d)
class(fks_mapping_default_t), intent(in) :: map
type(vector4_t), intent(in), dimension(:) :: p
integer, intent(in) :: i, j
real(default) :: d
real(default) :: sqrts
real :: y
real :: E1, E2
if (i /= j .and. (i > 2 .or. j > 2)) then
if (i == 0 .or. j == 0) then
if (j == 0) then
E1 = energy (p(i))
y = polar_angle_ct (p(i))
else
E1 = energy (p(j))
y = polar_angle_ct(p(j))
end if
d = (E1**2 * (1-y**2))**map%exp_2
else
E1 = energy(p(i))
E2 = energy(p(j))
y = enclosed_angle_ct (p(i), p(j))
sqrts = (p(1)+p(2))**1
d = (2*p(i)*p(j) * E1*E2 / (E1 + E2)**2)**map%exp_1
end if
else if (i == j) then
call msg_fatal ("Invalid FKS region: Emitter equals FKS parton!")
else
!!! case i,j <= 2 not yet implemented
d = 0
end if
end function fks_mapping_default_dij
@ %def fks_mapping_default_dij
@
<<fks regions: fks mapping 1: TBP>>=
procedure :: compute_sumdij => fks_mapping_default_compute_sumdij
<<fks regions: procedures>>=
function fks_mapping_default_compute_sumdij (map, sregion, p) result (d)
class(fks_mapping_default_t), intent(in) :: map
type(singular_region_t), intent(inout) :: sregion
type(vector4_t), intent(in), dimension(:) :: p
real(default) :: d
integer :: i, k, l
associate (ftuples => sregion%flst_allreg)
d = 0
do i = 1, sregion%nregions
call ftuples(i)%get (k, l)
d = d + 1.0/map%dij (p, k, l)
end do
end associate
end function fks_mapping_default_compute_sumdij
@ %def fks_mapping_default_compute_sumdij
@
<<fks regions: fks mapping 1: TBP>>=
procedure :: svalue => fks_mapping_default_svalue
<<fks regions: procedures>>=
function fks_mapping_default_svalue (map, p, i, j) result (value)
class(fks_mapping_default_t), intent(in) :: map
type(vector4_t), intent(in), dimension(:) :: p
integer, intent(in) :: i, j
real(default) :: value
value = 1._default / (map%dij (p, i, j) * map%sumdij)
end function fks_mapping_default_svalue
@ %def fks_mapping_default_svalue
@
<<fks regions: fks mapping 1: TBP>>=
procedure :: dij_soft => fks_mapping_default_dij_soft
<<fks regions: procedures>>=
function fks_mapping_default_dij_soft (map, p_born, p_soft, em) result (d)
class(fks_mapping_default_t), intent(in) :: map
type(vector4_t), intent(in), dimension(:) :: p_born
type(vector4_t), intent(in) :: p_soft
integer, intent(in) :: em
real(default) :: d
d = (2*p_born(em)*p_soft / energy(p_born(em)))**map%exp_1
end function fks_mapping_default_dij_soft
@ %def fks_mapping_default_dij_soft
@
<<fks regions: fks mapping 1: TBP>>=
procedure :: compute_sumdij_soft => fks_mapping_default_compute_sumdij_soft
<<fks regions: procedures>>=
function fks_mapping_default_compute_sumdij_soft (map, sregion, p_born, p_soft) result (d)
class(fks_mapping_default_t), intent(in) :: map
type(singular_region_t), intent(inout) :: sregion
type(vector4_t), intent(in), dimension(:) :: p_born
type(vector4_t), intent(in) :: p_soft
real(default) :: d
integer :: i, k, l
integer :: nlegs
d = 0
nlegs = size (sregion%flst_real%flst)
associate (ftuples => sregion%flst_allreg)
do i = 1, sregion%nregions
call ftuples(i)%get (k,l)
if (l == nlegs) then
d = d + 1._default/map%dij_soft (p_born, p_soft, k)
end if
end do
end associate
end function fks_mapping_default_compute_sumdij_soft
@ %def fks_mapping_default_compute_sumdij_soft
@
<<fks regions: fks mapping 1: TBP>>=
procedure :: svalue_soft => fks_mapping_default_svalue_soft
<<fks regions: procedures>>=
function fks_mapping_default_svalue_soft (map, p_born, p_soft, em) result (value)
class(fks_mapping_default_t), intent(in) :: map
type(vector4_t), intent(in), dimension(:) :: p_born
type(vector4_t), intent(in) :: p_soft
integer, intent(in) :: em
real(default) :: value
value = 1._default/(map%sumdij_soft*map%dij_soft (p_born, p_soft, em))
end function fks_mapping_default_svalue_soft
@ %def fks_mapping_default_svalue_soft
@
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Virtual contribution to the cross section}
<<[[virtual.f90]]>>=
<<File header>>
module virtual
<<Use kinds>>
<<Use strings>>
use constants, only: pi, twopi
use diagnostics
use pdg_arrays
use model_data
use physics_defs
use sm_physics
use lorentz
use flavors
use fks_regions
<<Standard module head>>
<<virtual: public>>
<<virtual: types>>
contains
<<virtual: procedures>>
end module virtual
@ %def virtual
<<virtual: public>>=
public :: virtual_t
<<virtual: types>>=
type :: virtual_t
real(default) :: Q
real(default), dimension(:,:), allocatable :: I
real(default) :: vfin
real(default) :: sqme_cc
real(default) :: sqme_virt
real(default), dimension(:,:), allocatable :: gamma_0, gamma_p, c_flv
real(default) :: ren_scale2
integer, dimension(:), allocatable :: n_is_neutrinos
integer :: nlegs, nflv
logical :: bad_point
logical :: use_internal_color_correlations
contains
<<virtual: virtual: TBP>>
end type virtual_t
@ %def virtual_t
<<virtual: virtual: TBP>>=
procedure :: init => virtual_init
<<virtual: procedures>>=
subroutine virtual_init (object, flv_born)
class(virtual_t), intent(inout) :: object
integer, intent(in), dimension(:,:) :: flv_born
integer :: nlegs, nflv
integer :: i_flv
type(flavor_t) :: flv
object%nlegs = size (flv_born, 1); object%nflv = size (flv_born, 2)
allocate (object%I (object%nlegs, object%nlegs))
allocate (object%gamma_0 (object%nlegs, object%nflv), &
object%gamma_p (object%nlegs, object%nflv), &
object%c_flv (object%nlegs, object%nflv))
call object%init_constants (flv_born)
allocate (object%n_is_neutrinos (object%nflv))
object%n_is_neutrinos = 0
do i_flv = 1, object%nflv
if (is_neutrino (flv_born(1, i_flv))) &
object%n_is_neutrinos(i_flv) = object%n_is_neutrinos(i_flv) + 1
if (is_neutrino (flv_born(2, i_flv))) &
object%n_is_neutrinos(i_flv) = object%n_is_neutrinos(i_flv) + 1
end do
contains
function is_neutrino (flv) result (neutrino)
integer, intent(in) :: flv
logical :: neutrino
neutrino = (abs(flv)==12 .or. abs(flv)==14 .or. abs(flv)==16)
end function is_neutrino
end subroutine virtual_init
@ %def virtual_init
@ Write down constant definition somewhere
<<virtual: virtual: TBP>>=
procedure :: init_constants => virtual_init_constants
<<virtual: procedures>>=
subroutine virtual_init_constants (object, flv_born)
class(virtual_t), intent(inout) :: object
integer, intent(in), dimension(:,:) :: flv_born
integer :: i_part, i_flv
integer, parameter :: nf = 1
do i_flv = 1, size (flv_born, 2)
do i_part = 1, size (flv_born, 1)
if (is_gluon (flv_born(i_part, i_flv))) then
object%gamma_0(i_part, i_flv) = (11*ca - 2*nf)/6
object%gamma_p(i_part, i_flv) = (67.0/9 - 2*pi**2/3)*ca - 23.0/18*nf
object%c_flv(i_part, i_flv) = ca
else if (is_quark (abs(flv_born(i_part, i_flv)))) then
object%gamma_0(i_part, i_flv) = 1.5*cf
object%gamma_p(i_part, i_flv) = (6.5 - 2*pi**2/3)*cf
object%c_flv(i_part, i_flv) = cf
else
object%gamma_0(i_part, i_flv) = 0
object%gamma_p(i_part, i_flv) = 0
object%c_flv(i_part, i_flv) = 0
end if
end do
end do
end subroutine virtual_init_constants
@ %def virtual_init_constants
@ Set the renormalization scale. If the input is zero, use the
center-of-mass energy.
<<virtual: virtual: TBP>>=
procedure :: set_ren_scale => virtual_set_ren_scale
<<virtual: procedures>>=
subroutine virtual_set_ren_scale (object, p, ren_scale)
class(virtual_t), intent(inout) :: object
type(vector4_t), dimension(:), intent(in) :: p
real(default), intent(in) :: ren_scale
if (ren_scale > 0) then
object%ren_scale2 = ren_scale**2
else
object%ren_scale2 = (p(1)+p(2))**2
end if
end subroutine virtual_set_ren_scale
@ %def virtual_set_ren_scale
@ The virtual-subtracted matrix element is given by the equation
\begin{equation}
\label{virt_sub}
\mathcal{V} = \frac{\alpha_s}{2\pi}\left(\mathcal{Q}\mathcal{B} +
\sum \mathcal{I}_{ij}\mathcal{B}_{ij} + \mathcal{V}_{fin}\right),
\end{equation}
where the quantity $\mathcal{Q}$ is given by
\begin{equation}
\begin{split}
\mathcal{Q} = \sum_{i=1}^n&\left[\gamma'_{f_i} - \log\frac{s}{Q^2}
\left(\gamma_{f_i} - 2C_{f_i}\log\frac{2E_i}{\sqrt{s}}\right)\right.\\
&\left. +2C_{f_i}\log^2\frac{2E_i}{\sqrt{s}} - 2\gamma_{f_i}\frac{2E_i}{\sqrt{s}}\right]\\
&- \log\frac{\mu_F^2}{Q^2}\left(\gamma_{f+} + \gamma_{f-}\right).
\end{split}
\label{virt_Q}
\end{equation}
The expressions for $\mathcal{I}_{ij}$ can be found in equations (\ref{I_00}), (\ref{I_mm}), (\ref{I_0m}),
depending on whether the particles involved in the radiation process are massive or massless.
<<virtual: virtual: TBP>>=
procedure :: evaluate => virtual_evaluate
<<virtual: procedures>>=
subroutine virtual_evaluate &
(object, reg_data, i_flv, alpha_s, p_born, born, b_ij)
class(virtual_t), intent(inout) :: object
type(region_data_t), intent(in) :: reg_data
integer, intent(in) :: i_flv
real(default), intent(in) :: alpha_s
type(vector4_t), intent(inout), dimension(:), allocatable :: p_born
real(default), intent(in) :: born
real(default), intent(in), dimension(:,:,:), allocatable :: b_ij
integer :: i, j, alr
integer :: nlegs
real(default) :: BI
if (object%bad_point) then
object%sqme_virt = 0
else
BI = 0
alr = find_first_matching_uborn (reg_data, i_flv)
associate (flst_born => reg_data%regions(alr)%flst_uborn)
call object%compute_Q (p_born, i_flv, flst_born%massive)
do i = 1, object%nlegs
do j = 1, object%nlegs
if (i /= j) then
if (flst_born%colored(i) .and. flst_born%colored(j)) then
call object%compute_I (p_born, flst_born%massive, i, j)
BI = BI + b_ij (i,j,reg_data%regions(alr)%uborn_index) * &
object%I(i,j)
end if
end if
end do
end do
end associate
if (object%use_internal_color_correlations) BI = BI*born
!!! A factor of alpha_s/twopi is assumed to be included in vfin
object%sqme_virt = alpha_s/twopi * (object%Q*born + BI) + object%vfin
if (object%n_is_neutrinos(i_flv) > 0) &
object%sqme_virt = object%sqme_virt * object%n_is_neutrinos(i_flv)*2
end if
contains
function find_first_matching_uborn (reg_data, i_proc) result (alr_out)
type(region_data_t), intent(in) :: reg_data
integer, intent(in) :: i_proc
integer :: alr_out
integer :: k
alr_out = 0
do k = 1, reg_data%n_regions
alr_out = alr_out+1
if (reg_data%regions(k)%uborn_index == i_proc) exit
end do
end function find_first_matching_uborn
end subroutine virtual_evaluate
@ %def virtual_evaluate
@
<<virtual: virtual: TBP>>=
procedure :: compute_vfin_test => virtual_compute_vfin_test
<<virtual: procedures>>=
subroutine virtual_compute_vfin_test (object, p_born, sqme_born)
class(virtual_t), intent(inout) :: object
type(vector4_t), intent(in), dimension(:) :: p_born
real(default), intent(in) :: sqme_born
real(default) :: s
s = (p_born(1)+p_born(2))**2
!!! ----NOTE: Test implementation for e+ e- -> u ubar
object%vfin = sqme_born * cf * &
(pi**2 - 8 + 3*log(s/object%ren_scale2) - log(s/object%ren_scale2)**2)
object%bad_point = .false.
end subroutine virtual_compute_vfin_test
@ %def virtual_compute_vfin
@
<<virtual: virtual: TBP>>=
procedure :: set_vfin => virtual_set_vfin
<<virtual: procedures>>=
subroutine virtual_set_vfin (object, vfin)
class(virtual_t), intent(inout) :: object
real(default) :: vfin
object%vfin = vfin
end subroutine virtual_set_vfin
@ %def virtual_set_vfin
@
<<virtual: virtual: TBP>>=
procedure :: set_bad_point => virtual_set_bad_point
<<virtual: procedures>>=
subroutine virtual_set_bad_point (object, value)
class(virtual_t), intent(inout) :: object
logical, intent(in) :: value
object%bad_point = value
end subroutine virtual_set_bad_point
@ %def virtual_set_bad_point
@
<<virtual: virtual: TBP>>=
procedure :: compute_Q => virtual_compute_Q
<<virtual: procedures>>=
subroutine virtual_compute_Q (object, p_born, i_flv, massive)
class(virtual_t), intent(inout) :: object
type(vector4_t), intent(in), dimension(:), allocatable :: p_born
integer, intent(in) :: i_flv
logical, dimension(:), intent(in) :: massive
real(default) :: sqrts, E
real(default) :: s1, s2, s3, s4
integer :: i
sqrts = sqrt ((p_born(1)+p_born(2))**2)
!!! ---- NOTE: Implementation only works for lepton collisions.
!!! This implies that both the summand containing log(s/q**2)
!!! and (γ_fp + γ…vfm) vanish.
!!! Also, s = (p_born(1)+p_born(2))**2
object%Q = 0
do i = 1, object%nlegs
if (.not. massive (i)) then
s1 = object%gamma_p(i, i_flv)
E = vector4_get_component (p_born(i), 0)
s2 = log(sqrts**2/object%ren_scale2)*&
(object%gamma_0(i, i_flv)-2 * object%c_flv(i, i_flv) * log(2*E/sqrts))
s3 = 2*log(2*E/sqrts)**2*object%c_flv(i, i_flv)
s4 = 2*log(2*E/sqrts)*object%gamma_0(i, i_flv)
object%Q = object%Q + s1 - s2 + s3 - s4
else
s1 = log(sqrts**2/object%ren_scale2)
s2 = 0.5*I_m_eps (p_born(i))
object%Q = object%Q - object%c_flv(i, i_flv) * (s1-s2)
end if
end do
end subroutine virtual_compute_Q
@ %def virtual_compute_Q
@ The following code implements the $\mathcal{I}_{ij}$-function appearing in eq. blub.
They are defined as follows:
\begin{itemize}
\item[Massles-Massles Case]
\begin{equation}
\begin{split}
\mathcal{I}_{ij} &= \frac{1}{2}\log^2\frac{s}{Q^2} + \log\frac{s}{Q^2}\log\frac{k_ik_j}{2E_iE_j}
- \rm{Li}_2\left(\frac{k_ik_j}{2E_iE_j}\right) \\
&+ \frac{1}{2}\log^2\frac{k_ik_j}{2E_iE_j} - \log\left(1-\frac{k_ik_j}{2E_iE_j}\right)
\log\frac{k_ik_j}{2E_iE_j}.
\end{split}
\label{I_00}
\end{equation}
\item[Massive-Massive Case]
\begin{equation}
\mathcal{I}_{ij} = -\frac{1}{2}I_0(k_i, k_j)\log\frac{Q^2}{s} - \frac{1}{2}I_\epsilon(k_i,k_j)
\label{I_mm}
\end{equation}
with
\begin{equation*}
I_0(k_i, k_j) = \frac{1}{\beta}\log\frac{1+\beta}{1-\beta}, \qquad
\beta = \sqrt{1-\frac{k_i^2k_j^2}{(k_i \cdot k_j)^2}}
\end{equation*}
and a rather involved expression for $I_\epsilon$:
\begin{align*}
I_\epsilon(k_i, k_j) &= \left(K(z_j)-K(z_i)\right) \frac{1-\vec{\beta_i}\cdot\vec{\beta_j}}{\sqrt{a(1-b)}}, \\
\vec{\beta_i} &= \frac{\vec{k}_i}{k_i^0}, \\
a &= \beta_i^2 + \beta_j^2 - 2\vec{\beta}_i \cdot \vec{\beta}_j, \\
x_i &= \frac{\beta_i^2 -\vec{\beta}_i \cdot \vec{\beta}_j}{a}, \\
x_j &= \frac{\beta_j^2 -\vec{\beta}_i \cdot \vec{\beta}_j}{a} = 1-x_j, \\
b &= \frac{\beta_i^2\beta_j^2 - (\vec{\beta}_i\cdot\vec{\beta}_j)^2}{a}, \\
c &= \sqrt{\frac{b}{4a}}, \\
z_+ &= \frac{1+\sqrt{1-b}}{\sqrt{b}}, \\
z_- &= \frac{1-\sqrt{1-b}}{\sqrt{b}}, \\
z_i &= \frac{\sqrt{x_i^2 + 4c^2} - x_i}{2c}, \\
z_j &= \frac{\sqrt{x_j^2 + 4c^2} + x_j}{2c}, \\
K(z) = &-\frac{1}{2}\log^2\frac{(z-z_-)(z_+-z)}{(z_++z)(z_-+z)} - 2Li_2\left(\frac{2z_-(z_+-z)}{(z_+-z_-)(z_-+z)}\right) \\
&-2Li_2\left(-\frac{2z_+(z_-+z)}{(z_+-z_-)(z_+-z)}\right)
\end{align*}
\item[Massive-Massless Case]
\begin{equation}
\mathcal{I}_{ij} = \frac{1}{2}\left[\log^2\frac{Q}{s} - \frac{\pi^2}{6}\right] -\frac{1}{2}I_0(k_i,k_j)\log\frac{Q^2}{s}
-\frac{1}{2}I_\epsilon(k_i,k_j)
\label{I_0m}
\end{equation}
with
\begin{align*}
I_0(p,k) &= \log\frac{(\hat{p}\cdot\hat{k})^2}{\hat{k}^2}, \\
I_\varepsilon(p,k) &= -2\left[\frac{1}{4}\log^2\frac{1-\beta}{1+\beta} + \log\frac{\hat{p}\cdot\hat{k}}{1+\beta}\log\frac{\hat{p}\cdot\hat{k}}{1-\beta}
+Li_2\left(1-\frac{\hat{p}\cdot\hat{k}}{1+\beta}\right) + Li_2\left(1-\frac{\hat{p}\cdot\hat{k}}{1-\beta}\right)\right]
\end{align*}
using
\begin{equation*}
\hat{p} = \frac{p}{p^0}, \quad \hat{k} = \frac{k}{k^0}, \quad \beta = \frac{|\vec{k}|}{k_0}.
\end{equation*}
\end{itemize}
<<virtual: virtual: TBP>>=
procedure :: compute_I => virtual_compute_I
<<virtual: procedures>>=
subroutine virtual_compute_I (object, p_born, massive, i, j)
class(virtual_t), intent(inout) :: object
type(vector4_t), intent(in), dimension(:) :: p_born
logical, dimension(:), intent(in) :: massive
integer, intent(in) :: i, j
real(default) :: somu2
somu2 = (p_born(1)+p_born(2))**2/object%ren_scale2
if (massive(i) .and. massive(j)) then
object%I(i,j) = compute_Imm (p_born(i), p_born(j), somu2)
else if (.not.massive(i) .and. massive(j)) then
object%I(i,j) = compute_I0m (p_born(i), p_born(j), somu2)
else if (massive(i) .and. .not.massive(j)) then
object%I(i,j) = compute_I0m (p_born(j), p_born(i), somu2)
else
object%I(i,j) = compute_I00 (p_born(i), p_born(j), somu2)
end if
end subroutine virtual_compute_I
function compute_I00 (pi, pj, somu2) result (I)
type(vector4_t), intent(in) :: pi, pj
real(default), intent(in) :: somu2
real(default) :: I
real(default) :: Ei, Ej
real(default) :: pij, Eij
real(default) :: s
real(default) :: s1, s2, s3, s4, s5
real(default) :: arglog
real(default), parameter :: tiny_value = epsilon(1.0)
!!! ----NOTE: As above, only lepton collisions. Therefore, the
!!! first and second summand are not taken into account.
s1 = 0; s2 = 0; s3 = 0; s4 = 0; s5 = 0
Ei = vector4_get_component (pi, 0)
Ej = vector4_get_component (pj, 0)
pij = pi*pj; Eij = Ei*Ej
s1 = 0.5*log(somu2)**2
s2 = log(somu2)*log(pij/(2*Eij))
s3 = Li2 (pij / (2*Eij))
s4 = 0.5*log (pij / (2*Eij))**2
arglog = 1 - pij/(2*Eij)
if (arglog > tiny_value) then
s5 = log(arglog) * log(pij / (2*Eij))
else
s5 = 0
end if
I = s1 + s2 - s3 + s4 - s5
end function compute_I00
function compute_I0m (ki, kj, somu2) result (I)
type(vector4_t), intent(in) :: ki, kj
real(default), intent(in) :: somu2
real(default) :: I
real(default) :: logsomu
real(default) :: s1, s2, s3
s1 = 0; s2 = 0; s3 = 0
logsomu = log(somu2)
s1 = 0.5*(0.5*logsomu - pi**2/6)
s2 = 0.5*I_0m_0 (ki, kj)*logsomu
s3 = 0.5*I_0m_eps (ki, kj)
I = s1 + s2 - s3
end function compute_I0m
function compute_Imm (pi, pj, somu2) result (I)
type(vector4_t), intent(in) :: pi, pj
real(default), intent(in) :: somu2
real(default) :: I
real(default) :: s1, s2
s1 = 0.5*log(somu2)*I_mm_0(pi, pj)
s2 = 0.5*I_mm_eps(pi, pj)
I = s1 - s2
end function compute_Imm
@ %def virtual_compute_I
@
<<virtual: procedures>>=
function I_m_eps (p) result (I)
type(vector4_t), intent(in) :: p
real(default) :: I
real(default) :: beta
beta = space_part_norm (p)/p%p(0)
I = 2*log((1+beta)/(1-beta))/beta
end function I_m_eps
@ %def I_m_eps
@For $p^2=0$ and $k^2 \neq 0$, this computes the expression
\begin{equation*}
I_\epsilon(p,k) = -2\left[\frac{1}{4}\log^2\frac{1-\beta}{1+\beta}
+\log \frac{\hat{p} \cdot k}{1+\beta}\log{\hat{p} \cdot k}{1-\beta}
+Li_2\left(1-\frac{\hat{p} \cdot k}{1+\beta}\right)
+Li_2\left(1-\frac{\hat{p} \cdot k}{1-\beta}\right)\right],
\end{equation*}
with $\hat{p} = \frac{p}{p^0}$ and $\beta = \frac{|\vec{k}|}{k^0}$.
<<virtual: procedures>>=
function I_0m_eps (p, k) result (I)
type(vector4_t), intent(in) :: p, k
real(default) :: I
type(vector4_t) :: pp, kp
real(default) :: beta
pp = p/energy(p); kp = k/energy(k)
beta = sqrt (1-pp*kp)
I = -2*(log((1-beta)/(1+beta))**2/4 + log((pp*kp)/(1+beta))*log((pp*kp)/(1-beta)) &
+ Li2(1-(pp*kp)/(1+beta)) + Li2(1-(pp*kp)/(1-beta)))
end function I_0m_eps
@ %def I_0m_eps
@For $p^2=0$ and $k^2 \neq 0$, computes the expression
\begin{equation*}
I_0(p,k) = \log\frac{(\hat{p}\cdot\hat{k})^2}{\hat{k}^2}
\end{equation*}
<<virtual: procedures>>=
function I_0m_0 (p, k) result (I)
type(vector4_t), intent(in) :: p, k
real(default) :: I
type(vector4_t) :: pp, kp
pp = p/energy(p); kp = k/energy(k)
I = log((pp*kp)**2/kp**2)
end function I_0m_0
@ %def I_0m_0
@ For $k_1^2 \neq 0$ and $k_2^2 \neq 0$, computes the expression
\begin{equation*}
I_\epsilon(k_1, k_2) = \left[K(z_2) - K(z_1)\right]
\frac{1-\vec{\beta}_1\cdot\vec{\beta}_2}{\sqrt{a(1-b)}},
\end{equation*}
where $\vec{\beta}_i = \frac{\vec{k}_i}{k_i^0}$. Further
\begin{align*}
a &= \beta_1^2 + \beta_2^2 - 2\vec{\beta}_1\cdot\vec{\beta}_2, \\
b &= \frac{\beta_1^2\beta_2^2 - \left(\vec{\beta}_1\cdot\vec{\beta}_2\right)^2}{a}, \\
c &= \sqrt{\frac{b}{4a}}, \\
K(z) = &-\frac{1}{2}\log^2\frac{(z-z_-)(z_+-z)}{(z_++z)(z_-+z)}
- 2Li_2\left(\frac{2z_-(z_+-z)}{(z_+-z_-)(z_-+z)}\right) \\
&-2Li_2\left(-\frac{2z_+(z_-+z)}{(z_+-z_-)(z_+-z)}\right),\\
z_+ &= \frac{1+\sqrt{1-b}}{\sqrt{b}}, \\
z_- &= \frac{1-\sqrt{1-b}}{\sqrt{b}}, \\
x_1 &= \frac{\beta_1^2-\vec{\beta}_1\cdot \vec{\beta}_2}{a}, \\
x_2 &= 1-x_1
z_1 &= \frac{\sqrt{x_1^2+4c^2}-x_1}{2c}, \\
z_2 &= \frac{\sqrt{x_2^2+4c^2}+x_2}{2c}.
\end{align*}
<<virtual: procedures>>=
function I_mm_eps (p1, p2) result (I)
type(vector4_t), intent(in) :: p1, p2
real(default) :: I
type(vector3_t) :: beta1, beta2
real(default) :: a, b
real(default) :: zp, zm, z1, z2, x1, x2
real(default) :: zmb, z1b
real(default) :: K1, K2
beta1 = space_part (p1)/energy(p1)
beta2 = space_part (p2)/energy(p2)
a = beta1**2 + beta2**2 - 2*beta1*beta2
b = beta1**2*beta2**2 - (beta1*beta2)**2
x1 = beta1**2 - beta1*beta2
x2 = beta2**2 - beta1*beta2
zp = sqrt(a) + sqrt(a-b)
zm = sqrt(a) - sqrt(a-b)
zmb = 1/zp
z1 = sqrt(x1**2+b) - x1
z2 = sqrt(x2**2+b) + x2
z1b = 1/(sqrt(x1**2+b)+x1)
K1 = -0.5*log(((z1b-zmb)*(zp-z1))/((zp+z1)*(z1b+zmb)))**2 &
-2*Li2((2*zmb*(zp-z1))/((zp-zm)*(zmb+z1b))) &
-2*Li2((-2*zp*(zm+z1))/((zp-zm)*(zp-z1)))
K2 = -0.5*log(((z2-zm)*(zp-z2))/((zp+z2)*(z2+zm)))**2 &
-2*Li2((2*zm*(zp-z2))/((zp-zm)*(zm+z2))) &
-2*Li2((-2*zp*(zm+z2))/((zp-zm)*(zp-z2)))
I = (K2 - K1) * (1-beta1*beta2)/sqrt(a-b)
end function I_mm_eps
@ %def I_mm_eps
@For $k_1^2 \neq 0$ and $k_2^2 \neq 0$ ,computes the expression
\begin{equation*}
I_0(k_1, k_2) = \frac{1}{\beta} \log\frac{1+\beta}{1-\beta},
\quad
\beta = \sqrt{1-\frac{k_1^2k_2^2}{(k_1 \cdot k_2)^2}}
\end{equation*}
<<virtual: procedures>>=
function I_mm_0 (k1, k2) result (I)
type(vector4_t), intent(in) :: k1, k2
real(default) :: I
real(default) :: beta
beta = sqrt (1-k1**2*k2**2/(k1*k2)**2)
I = log((1+beta)/(1-beta))/beta
end function I_mm_0
@ %def I_mm_0
@
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Real Subtraction}
<<[[real_subtraction.f90]]>>=
<<File header>>
module real_subtraction
<<Use kinds>>
<<Use strings>>
- use constants, only: pi, twopi
+ use constants
use diagnostics
use pdg_arrays
use model_data
use physics_defs
use sm_physics
use lorentz
use flavors
use fks_regions
<<Standard module head>>
<<real subtraction: public>>
<<real subtraction: types>>
contains
<<real subtraction: procedures>>
end module real_subtraction
@ %def real_subtraction
@
\subsubsection{Soft subtraction terms}
In the soft limit, the real matrix element behaves as
\begin{equation*}
\mathcal{R}_{\rm{soft}} = 4\pi\alpha_s \left[\sum_{i \neq j} \mathcal{B}_{ij} \frac{k_i \cdot k_j}{(k_i \cdot k)(k_j \cdot k)} - \mathcal{B} \sum_{i} \frac{k_i^2}{(k_i \cdot k)^2}C_i\right],
\end{equation*}
where $k$ denotes the momentum of the emitted parton. The quantity $\mathcal{B}_{ij}$ is called the color-correlated Born matrix element defined as
\begin{equation*}
\mathcal{B}_{ij} = \frac{1}{2s} \sum_{\stackrel{colors}{spins}} \mathcal{M}_{\{c_k\}}\left(\mathcal{M}^\dagger_{\{c_k\}}\right)_{\stackrel{c_i \rightarrow c_i'}{c_j \rightarrow c_j'}} T^a_{c_i,c_i'} T^a_{c_j,c_j'}.
\end{equation*}
<<real subtraction: types>>=
type :: soft_subtraction_t
real(default), dimension(:), allocatable :: value
type(region_data_t) :: reg_data
integer :: nlegs_born, nlegs_real
real(default), dimension(:,:), allocatable :: momentum_matrix
logical :: use_internal_color_correlations = .true.
logical :: use_internal_spin_correlations = .false.
contains
<<real subtraction: soft sub: TBP>>
end type soft_subtraction_t
@ %def soft_subtraction_t
<<real subtraction: soft sub: TBP>>=
procedure :: init => soft_subtraction_init
<<real subtraction: procedures>>=
subroutine soft_subtraction_init (sub_soft, reg_data, nlegs_born, &
nlegs_real)
class(soft_subtraction_t), intent(inout) :: sub_soft
type(region_data_t), intent(in) :: reg_data
integer, intent(in) :: nlegs_born, nlegs_real
sub_soft%reg_data = reg_data
sub_soft%nlegs_born = nlegs_born
sub_soft%nlegs_real = nlegs_real
allocate (sub_soft%value (reg_data%n_regions))
allocate (sub_soft%momentum_matrix &
(nlegs_born, nlegs_born))
end subroutine soft_subtraction_init
@ %def soft_subtraction_init
@ The treatment of the momentum $k$ follows the discussion about the soft limit of the partition functions (ref????).
The parton momentum is pulled out, $k = E \hat{k}$. In fact, we will substitute $\hat{k}$ for $k$ throughout the code,
because the energy will factor out of the equation when the soft $\mathcal{S}$-function is multiplied.
The soft momentum is a unit vector, because $k^2 = \left(k^0\right)^2 - \left(k^0\right)^2\hat{\vec{k}}^2 = 0$.
<<real subtraction: soft sub: TBP>>=
procedure :: create_softvec_fsr => soft_subtraction_create_softvec_fsr
<<real subtraction: procedures>>=
function soft_subtraction_create_softvec_fsr &
(sub_soft, p_born, y, phi, emitter) result (p_soft)
class(soft_subtraction_t), intent(inout) :: sub_soft
type(vector4_t), intent(in), dimension(:) :: p_born
real(default), intent(in) :: y, phi
integer, intent(in) :: emitter
real(default) :: cpsi
type(vector4_t) :: p_soft
type(vector3_t) :: dir
type(lorentz_transformation_t) :: rot
@ The soft momentum is constructed by first creating a unit vector parallel to the emitter's Born momentum. This unit vector
is then rotated about the corresponding angles $y$ and $\phi$.
<<real subtraction: procedures>>=
p_soft%p(0) = 1._default
p_soft%p(1:3) = p_born(emitter)%p(1:3) / space_part_norm (p_born(emitter))
dir = create_orthogonal (space_part (p_born(emitter)))
rot = rotation (y, sqrt(1-y**2), dir)
p_soft = rot*p_soft
if (phi /= 0) then
dir = space_part (p_born(emitter)) / &
space_part_norm (p_born(emitter))
rot = rotation (cos(phi), sin(phi), dir)
p_soft = rot*p_soft
end if
end function soft_subtraction_create_softvec_fsr
@ %def soft_subtraction_create_softvec_fsr
@ Computation of $\mathcal{R}_{\rm{soft}}$:
<<real subtraction: soft sub: TBP>>=
procedure :: compute => soft_subtraction_compute
<<real subtraction: procedures>>=
subroutine soft_subtraction_compute (sub_soft, p_born, p_real, &
sqme_born, born_ij, y, y_soft, phi, alpha_s_born, alr, emitter)
class(soft_subtraction_t), intent(inout) :: sub_soft
type(vector4_t), intent(in), dimension(:) :: p_born
type(vector4_t), intent(in), dimension(:) :: p_real
real(default), intent(in) :: sqme_born
real(default), intent(in), dimension(:,:) :: born_ij
real(default), intent(in) :: y, y_soft, phi
real(default), intent(in) :: alpha_s_born
integer, intent(in) :: alr, emitter
type(vector4_t) :: p_soft
real(default) :: s_alpha_soft
real(default) :: q02
real(default) :: kb
integer :: i, j
p_soft = sub_soft%create_softvec_fsr (p_born, y_soft, phi, emitter)
s_alpha_soft = sub_soft%reg_data%get_svalue_soft &
(p_born, p_soft, alr, emitter)
call sub_soft%compute_momentum_matrix (p_born, p_soft)
sub_soft%value(alr) = 4*pi*alpha_s_born * s_alpha_soft
kb = 0._default
do i = 1, size (p_born)
do j = 1, size (p_born)
kb = kb + sub_soft%momentum_matrix (i,j) * &
born_ij (i,j)
end do
end do
sub_soft%value(alr) = sub_soft%value(alr)*kb
q02 = 4* vector4_get_component (p_born(1), 0) * &
vector4_get_component (p_born(2), 0)
sub_soft%value(alr) = 4/q02 * (1-y) * sub_soft%value(alr)
end subroutine soft_subtraction_compute
@ %def soft_subtraction_compute
@ We have to multiply this with $\xi^2(1-y)$. Further, when applying
the soft $\mathcal{S}$-function, the energy of the radiated particle
is factored out. Thus we have $\xi^2/E_{em}^2(1-y) = 4/q_0^2(1-y)$.
<<nlo data: soft subtraction computation>>=
@ Computes the quantity $\mathcal{K}_{ij} = \frac{k_i \cdot
k_j}{(k_i\cdot k)(k_j\cdot k)}$.
<<real subtraction: soft sub: TBP>>=
procedure :: compute_momentum_matrix => &
soft_subtraction_compute_momentum_matrix
<<real subtraction: procedures>>=
subroutine soft_subtraction_compute_momentum_matrix &
(sub_soft, p_born, p_soft)
class(soft_subtraction_t), intent(inout) :: sub_soft
type(vector4_t), intent(in), dimension(:) :: p_born
type(vector4_t), intent(in) :: p_soft
real(default) :: num, deno1, deno2
integer :: i, j
do i = 1, sub_soft%nlegs_born
do j = 1, sub_soft%nlegs_born
if (i <= j) then
num = p_born(i) * p_born(j)
deno1 = p_born(i)*p_soft
deno2 = p_born(j)*p_soft
sub_soft%momentum_matrix(i,j) = num/(deno1*deno2)
else
!!! momentum matrix is symmetric.
sub_soft%momentum_matrix(i,j) = sub_soft%momentum_matrix(j,i)
end if
end do
end do
end subroutine soft_subtraction_compute_momentum_matrix
@ %def soft_subtraction_compute_momentum_matrx
@
\subsection{Collinear and soft-collinear subtraction terms}
This data type deals with the calculation of the collinear and
soft-collinear contribution to the cross section.
<<real subtraction: types>>=
type :: coll_subtraction_t
real(default), dimension(:), allocatable :: value
real(default), dimension(:), allocatable :: value_soft
integer :: n_alr
real(default), dimension(0:3,0:3) :: b_munu
real(default) , dimension(0:3,0:3) :: k_perp_matrix
contains
<<real subtraction: coll sub: TBP>>
end type coll_subtraction_t
@ %def coll_subtraction_t
@
<<real subtraction: coll sub: TBP>>=
procedure :: init => coll_subtraction_init
<<real subtraction: procedures>>=
subroutine coll_subtraction_init (coll_sub, n_alr)
class(coll_subtraction_t), intent(inout) :: coll_sub
integer, intent(in) :: n_alr
coll_sub%n_alr = n_alr
allocate (coll_sub%value (n_alr))
allocate (coll_sub%value_soft (n_alr))
end subroutine coll_subtraction_init
@ %def coll_subtraction_init
@ To compute the collinear limit of $\mathcal{R}$, we follow the original FKS-paper. Here, the real amplitude is supposed to factorize
in the collinear limit. considering a splitting $g \rightarrow g(i)g(j)$,
\begin{equation*}
\mathcal{A}^{(n)}(h_i,h_j,\{h_l\} \stackrel{i \parallel j}{\rightarrow} g_s \sum_{d_e} \sum_{h_e} C(d_e, b, c)
S_{gg}^{h_eh_ih_j}(z)\mathcal{A}_{d_e}^{(n-1)} (h_e, \{h_l\})
\end{equation*}
\textbf{Explain quantities.} Evaluating this expression leads to
\begin{equation*}
|\mathcal{A}^{(n)}(h_i,h_j,\{h_l\})|^2 \stackrel{i \parallel j}{\rightarrow} |\mathcal{N}^{(n)}(h_i,h_j,\{h_l\}|^2 + \mathcal{R}(h_i,h_j,\{h_l\})|^2,
\end{equation*}
where $\mathcal{R}$ contains all contributions with spin-correlated amplitudes, i.e. terms like $\mathcal{A}_{d_e}^{(n-1)}(+,\{h_l\})\left(\mathcal{A}_{d_e}^{(n-1)}(-,\{h_l\})\right)^*$. Explicitly,
\begin{equation*}
\sum_{h_i,h_j,\{h_l\}}\sum_{\{d_l\}} |\mathcal{N}^{(n)}(h_i,h_j,\{h_l\})|^2 = \frac{4\pi\alpha_s}{k_i \cdot k_j} P_{gg}(z)|\mathcal{A}^{(n-1)}|^2,
\end{equation*}
and
\begin{align*}
\sum_{h_i,h_j,\{h_l\}}\sum_{\{d_l\}} &\mathcal{R}(h_i,h_j,\{h_l\}) = 32\pi\alpha_s C_A z(1-z) \\
&\times \underbrace{Re\left\{\frac{\langle ij \rangle}{[ij]}\sum_{\{h_l\}}\sum_{d_e, \{d_l\}} \mathcal{A}^{(n-1)}_{d_e}(+,\{h_l\})\left(\mathcal{A}^{(n-1)}_{d_e}(-,\{h_l\})\right)^*\right\}}
_{-\tilde{\mathcal{M}}^{(n-1)}/(2k_i\cdot k_j)}.
\end{align*}
This yields
\begin{equation}
\mathcal{M}^{(n)} \stackrel{i \parallel j}{\rightarrow} \frac{4\pi\alpha_s}{k_i \cdot k_j} P_{gg}(z)\mathcal{M}^{(n-1)} - \frac{16\pi\alpha_s}{k_i \cdot k_j} C_A z(1-z)\tilde{\mathcal{M}}^{(n-1)}.
\label{coll1}
\end{equation}
The equivalent expression for a $g \rightarrow qq$-splitting is given by
\begin{equation}
\mathcal{M}^{(n)} \stackrel{i \parallel j}{\rightarrow} \frac{4\pi\alpha_s}{k_i \cdot k_j} P_{qg}(z)\mathcal{M}^{(n-1)} + \frac{16\pi\alpha_s}{k_i \cdot k_j} T_F z(1-z)\tilde{\mathcal{M}}^{(n-1)}.
\label{coll2}
\end{equation}
<<real subtraction: coll sub: TBP>>=
procedure :: compute => coll_subtraction_compute
<<real subtraction: procedures>>=
subroutine coll_subtraction_compute &
(coll_sub, sregion, p_born, sqme_born, sqme_born_sc, &
xi, alpha_s, alr, soft_in)
class(coll_subtraction_t), intent(inout) :: coll_sub
type(singular_region_t), intent(in) :: sregion
type(vector4_t), intent(in), dimension(:) :: p_born
real(default), intent(in) :: sqme_born
real(default), intent(in) :: sqme_born_sc
real(default), intent(in) :: xi, alpha_s
integer, intent(in) :: alr
@ The function obtains a flag to indicate whether the limit to be computed is also soft. Moreover, this consideration
explains the structure of the code. In the soft limit, we find $z \rightarrow 0$ as well as $\xi \rightarrow 0$. However,
the quantity $z/\xi$ is finite, because
\begin{equation*}
\frac{z}{\xi} = \frac{p_{rad}^0}{\bar{p}_{em}^0} \frac{q^0}{2p_{em}^0} = \frac{q^0}{2\bar{p}^0_{em}}.
\end{equation*}
Thus, all expressions are written in terms of this well-behaved quantity. Recalling that there is a prefactor of $(\xi^2z)^{-1}$,
it is necessary to expand with $z$ once.
<<real subtraction: procedures>>=
logical, intent(in), optional :: soft_in
real(default) :: res
real(default) :: q0, z, p0
real(default) :: zoxi, onemz
real(default) :: pggz, pqgz
integer :: nlegs, emitter
integer :: flv_em, flv_rad
logical :: soft
if (present (soft_in)) then
soft = soft_in
else
soft = .false.
end if
nlegs = size (sregion%flst_real%flst)
emitter = sregion%emitter
flv_rad = sregion%flst_real%flst(nlegs)
flv_em = sregion%flst_real%flst(emitter)
p0 = vector4_get_component (p_born(emitter),0)
if (sregion%emitter <= 2) then
coll_sub%value(alr) = 0
return
else
q0 = vector4_get_component (p_born(1), 0) + &
vector4_get_component (p_born(2), 0)
!!! Here, z corresponds to 1-z in the formulas of arXiv:1002.2581;
!!! the integrand is symmetric under this variable change
zoxi = q0/(2*p0)
z = xi*zoxi; onemz = 1-z
if (is_gluon(flv_em) .and. is_gluon(flv_rad)) then
@ Implementation of equation (\ref{coll1}). Note that an additional factor $z$, so that in the last step, the whole expression is divided by $z/\xi$.
<<real subtraction: procedures>>=
pggz = 2*CA*(z**2*onemz + z**2/onemz + onemz)
res = pggz*sqme_born - 4*CA*z**2*onemz*sqme_born_sc
res = res/zoxi
else if (is_quark(abs(flv_em)) .and. is_quark (abs(flv_rad))) then
@ Equation \ref{coll2}
<<real subtraction: procedures>>=
pqgz = TR*z*(1-2*z*onemz)
res = pqgz*sqme_born + 4*TR*z**2*onemz*sqme_born_sc
res = res/zoxi
else if (is_quark (abs(flv_em)) .and. is_gluon (flv_rad)) then
res = sqme_born*CF*(1+onemz**2)/zoxi
else
call msg_fatal ('Impossible flavor structure in collinear counterterm!')
end if
end if
res = res /(p0**2*onemz*zoxi)
res = res * 4*pi*alpha_s
if (soft) then
coll_sub%value_soft (alr) = res
else
coll_sub%value (alr) = res
end if
end subroutine coll_subtraction_compute
@ %def coll_subtraction_compute
@ Here, $\xi = 0$ is already required.
<<real subtraction: coll sub: TBP>>=
procedure :: compute_soft_limit => coll_subtraction_compute_soft_limit
<<real subtraction: procedures>>=
subroutine coll_subtraction_compute_soft_limit &
(coll_sub, sregion, p_born, sqme_born, &
sqme_born_sc, xi, alpha_s, alr)
class(coll_subtraction_t), intent(inout) :: coll_sub
type(singular_region_t), intent(in) :: sregion
type(vector4_t), intent(in), dimension(:) :: p_born
real(default), intent(in) :: sqme_born
real(default), intent(in) :: sqme_born_sc
real(default) :: xi, alpha_s
integer, intent(in) :: alr
call coll_sub%compute (sregion, p_born, sqme_born, &
sqme_born_sc, xi, alpha_s, alr, .true.)
end subroutine coll_subtraction_compute_soft_limit
@ %def coll_subtraction_compute_soft_limit
@
\subsection{Real Subtraction}
Just a container for the real kinematics variables.
<<real subtraction: public>>=
public :: real_kinematics_t
<<real subtraction: types>>=
type :: real_kinematics_t
real(default) :: xi_tilde
real(default) :: phi
real(default), dimension(:), allocatable :: xi_max, y
real(default), dimension(3) :: jac
type(vector4_t), dimension(:), allocatable :: p_real
real(default), dimension(:), allocatable :: jac_rand
real(default), dimension(:), allocatable :: y_soft
end type real_kinematics_t
@ %def real_kinematics_t
<<real subtraction: public>>=
public :: real_subtraction_t
<<real subtraction: types>>=
type :: real_subtraction_t
type(region_data_t) :: reg_data
+ type(vector4_t), dimension(:), allocatable :: p_born, p_real
+ type(real_kinematics_t) :: fks_kinematics
+ integer :: current_alr = 0
real(default), dimension(:), allocatable :: sqme_real_bare
real(default), dimension(:), allocatable :: sqme_born
real(default), dimension(:,:,:), allocatable :: sqme_born_cc
complex(default) :: me_sc
type(soft_subtraction_t) :: sub_soft
type(coll_subtraction_t) :: sub_coll
logical, dimension(:), allocatable :: sc_required
contains
<<real subtraction: real subtraction: TBP>>
end type real_subtraction_t
@ %def real_subtraction_t
@ Initializer
<<real subtraction: real subtraction: TBP>>=
procedure :: init => real_subtraction_init
<<real subtraction: procedures>>=
subroutine real_subtraction_init (rsub, reg_data, nlegs_born, &
nlegs_real)
class(real_subtraction_t), intent(inout) :: rsub
type(region_data_t), intent(in) :: reg_data
integer, intent(in) :: nlegs_born, nlegs_real
integer :: alr, i_uborn
rsub%reg_data = reg_data
+ allocate (rsub%p_born (nlegs_born))
+ allocate (rsub%p_real (nlegs_real))
allocate (rsub%sqme_real_bare (reg_data%n_flv_real))
allocate (rsub%sqme_born (reg_data%n_flv_born))
allocate (rsub%sqme_born_cc (nlegs_born, nlegs_born, &
reg_data%n_flv_born))
allocate (rsub%sc_required (reg_data%n_regions))
do alr = 1, reg_data%n_regions
i_uborn = reg_data%regions(alr)%uborn_index
rsub%sc_required(alr) = &
reg_data%flv_born(i_uborn)%count_particle (GLUON) > 0
end do
call rsub%sub_soft%init (reg_data, nlegs_born, nlegs_real)
call rsub%sub_coll%init (reg_data%n_regions)
end subroutine real_subtraction_init
@ %def real_subtraction_init
+@
+<<real subtraction: real subtraction: TBP>>=
+ procedure :: set_momenta => real_subtraction_set_momenta
+<<real subtraction: procedures>>=
+ subroutine real_subtraction_set_momenta (rsub, p_born, p_real)
+ class(real_subtraction_t), intent(inout) :: rsub
+ type(vector4_t), dimension(:) :: p_born, p_real
+ rsub%p_born = p_born; rsub%p_real = p_real
+ end subroutine real_subtraction_set_momenta
+
+@ %def real_subtraction_set_momenta
+@
+<<real subtraction: real subtraction: TBP>>=
+ procedure :: set_fks_kinematics => real_subtraction_set_fks_kinematics
+<<real subtraction: procedures>>=
+ subroutine real_subtraction_set_fks_kinematics (rsub, fks_kinematics)
+ class(real_subtraction_t), intent(inout) :: rsub
+ type(real_kinematics_t), intent(in) :: fks_kinematics
+ rsub%fks_kinematics = fks_kinematics
+ end subroutine real_subtraction_set_fks_kinematics
+
+@ %def real_subtraction_set_fks_kinematics
+@
+<<real subtraction: real subtraction: TBP>>=
+ procedure :: set_alr => real_subtraction_set_alr
+<<real subtraction: procedures>>=
+ subroutine real_subtraction_set_alr (rsub, alr)
+ class(real_subtraction_t), intent(inout) :: rsub
+ integer, intent(in) :: alr
+ rsub%current_alr = alr
+ end subroutine real_subtraction_set_alr
+
+@ %def real_subtraction_set_alr
@\subsection{The real contribution to the cross section}
In each singular region $\alpha$, the real contribution to $\sigma$ is
given by the second summand of eqn. \ref{fks: sub: complete},
\begin{equation}
\sigma^\alpha_{\text{real}} = \int d\Phi_n \int_0^{2\pi} d\phi
\int_{-1}^1 dy \int_0^{\xi_{\text{max}}} d\xi
\left(\frac{1}{\xi}\right)_+ \left(\frac{1}{1-y}\right)_+
\underbrace{\frac{J(\Phi_n, \xi, y, \phi)}{\xi}
\left[(1-y)\xi^2\mathcal{R}^\alpha(\Phi_{n+1})\right]}_{g^\alpha(\xi,y)}.
\end{equation}
Writing out the plus-distribution and introducing $\tilde{\xi} =
\xi/\xi_{\text{max}}$ to set the upper integration limit to 1, this
turns out to be equal to
\begin{equation}
\begin{split}
\sigma^\alpha_{\rm{real}} &= \int d\Phi_n \int_0^{2/pi}d\phi
\int_{-1}^1 \frac{dy}{1-y} \Bigg\{\int_0^1
d\tilde{\xi}\Bigg[\frac{g^\alpha(\tilde{\xi}\xi_{\rm{max}},y)}{\tilde{\xi}}
- \underbrace{\frac{g^\alpha(0,y)}{\tilde{\xi}}}_{\text{soft}} -
\underbrace{\frac{g^\alpha(\tilde{\xi}\xi_{\rm{max}},1)}{\tilde{\xi}}}_{\text{coll.}}
+
\underbrace{\frac{g^\alpha(0,1)}{\tilde{\xi}}}_{\text{coll.+soft}}\Bigg]
\\
&+ \left[\log\xi_{\rm{max}}(y)g^\alpha(0,y) - \log\xi_{\rm{max}}(1)g^\alpha(0,1)\right]\Bigg\}.
\end{split}
\end{equation}
This formula is implemented in \texttt{compute\_sqme\_real\_fin}
<<real subtraction: real subtraction: TBP>>=
procedure :: compute => real_subtraction_compute
<<real subtraction: procedures>>=
- function real_subtraction_compute (rsub, emitter, p_real, p_born, &
- fks_kinematics, alpha_s) result (sqme)
+ function real_subtraction_compute (rsub, emitter, alpha_s) result (sqme)
class(real_subtraction_t), intent(inout) :: rsub
integer, intent(in) :: emitter
- type(vector4_t), dimension(:), intent(in) :: p_real, p_born
- type(real_kinematics_t), intent(in) :: fks_kinematics
real(default) :: alpha_s
real(default) :: sqme
- real(default) :: xi, y, xi_max, phi
- real(default) :: xi_tilde
- real(default), dimension(3) :: jac
- real(default) :: jac_rand
- real(default) :: s_alpha
- real(default) :: sqme0
- real(default) :: sqme_soft, sqme_coll, sqme_cs, sqme_remn
- integer :: alr, i_real, i_uborn
+ integer :: alr
+
sqme = 0._default
- xi_tilde = fks_kinematics%xi_tilde
- xi_max = fks_kinematics%xi_max(emitter)
- y = fks_kinematics%y(emitter)
- phi = fks_kinematics%phi
- jac = fks_kinematics%jac
- jac_rand = fks_kinematics%jac_rand(emitter)
- LOOP_OVER_ALPHA_REGIONS: do alr = 1, size (rsub%reg_data%regions)
- associate (region => rsub%reg_data%regions(alr))
- i_uborn = region%uborn_index
+ !!! Loop over all singular regions
+ do alr = 1, size (rsub%reg_data%regions)
if (emitter == rsub%reg_data%regions(alr)%emitter) then
- i_real = region%real_index
- sqme0 = rsub%sqme_real_bare (i_real)
- xi = xi_tilde * xi_max
- sqme0 = sqme0 * xi**2/xi_tilde
- s_alpha = rsub%reg_data%get_svalue (p_real, alr, emitter)
- sqme0 = sqme0 * s_alpha * jac(1)
- sqme0 = sqme0 * region%mult
- sqme0 = sqme0 * region%double_fsr_factor (p_real)
- call rsub%compute_sub_soft (p_born, p_real, fks_kinematics, &
- alr, emitter, alpha_s)
- call rsub%compute_sub_coll (p_born, p_real, fks_kinematics, &
- alr, emitter, alpha_s)
- call rsub%compute_sub_coll_soft (p_born, p_real, &
- alr, emitter, alpha_s)
- sqme_soft = rsub%sub_soft%value(alr)
- sqme_coll = rsub%sub_coll%value(alr)
- sqme_cs = rsub%sub_coll%value_soft(alr)
+ call rsub%set_alr (alr)
+ sqme = sqme + rsub%evaluate_region (emitter, alpha_s)
+ end if
+ end do
+ sqme = sqme*rsub%get_phs_factor ()
+ end function real_subtraction_compute
+
+@ %def real_subtraction_compute
+@
+<<real subtraction: real subtraction: TBP>>=
+ procedure :: evaluate_region => real_subtraction_evaluate_region
+<<real subtraction: procedures>>=
+ function real_subtraction_evaluate_region (rsub, emitter, &
+ alpha_s) result (sqme)
+ class(real_subtraction_t), intent(inout) :: rsub
+ integer, intent(in) :: emitter
+ real(default), intent(in) :: alpha_s
+ real(default) :: sqme
+ integer :: i_real
+ real(default) :: sqme0, sqme_soft, sqme_coll, sqme_cs, sqme_remn
+ real(default) :: s_alpha
+ real(default) :: xi, xi_max, xi_tilde, y, phi
+ xi_tilde = rsub%fks_kinematics%xi_tilde
+ xi_max = rsub%fks_kinematics%xi_max(emitter)
+ xi = xi_tilde * xi_max
+ y = rsub%fks_kinematics%y(emitter)
+ phi = rsub%fks_kinematics%phi
+ associate (region => rsub%reg_data%regions(rsub%current_alr))
+ i_real = region%real_index
+ sqme0 = rsub%sqme_real_bare (i_real)
+ s_alpha = rsub%reg_data%get_svalue (rsub%p_real, rsub%current_alr, emitter)
+ sqme0 = sqme0 * s_alpha
+ sqme0 = sqme0 * region%mult
+ sqme0 = sqme0 * region%double_fsr_factor (rsub%p_real)
+ call rsub%compute_sub_soft (emitter, alpha_s)
+ call rsub%compute_sub_coll (emitter, alpha_s)
+ call rsub%compute_sub_coll_soft (emitter, alpha_s)
+ call rsub%get_limit_sqme (sqme_soft, sqme_coll, sqme_cs)
@ Supply the $\tilde{\xi}^{-1}(1-y)^{-1}$-factor; note that this is not done for the pure real matrix element because
it directly cancels out with the numerator. Further, the jacobians are supplied in their appropriate limits.
<<real subtraction: procedures>>=
- sqme_soft = sqme_soft/(1-y)/xi_tilde*jac(2)
- sqme_coll = sqme_coll/(1-y)/xi_tilde*jac(3)
- sqme_cs = sqme_cs/(1-y)/xi_tilde*jac(2)
- sqme_remn = (sqme_soft - sqme_cs)*log(xi_max)*xi_tilde
- sqme = sqme + sqme0 - sqme_soft - sqme_coll + sqme_cs + sqme_remn
- !!! This occurs if the result is NaN.
- if (sqme /= sqme) &
- call write_computation_status ()
- sqme = sqme*jac_rand
- end if
+ associate (jac => rsub%fks_kinematics%jac)
+ sqme0 = sqme0 * xi**2/xi_tilde * jac(1)
+ sqme_soft = sqme_soft/(1-y)/xi_tilde * jac(2)
+ sqme_coll = sqme_coll/(1-y)/xi_tilde * jac(3)
+ sqme_cs = sqme_cs/(1-y)/xi_tilde * jac(2)
end associate
- end do LOOP_OVER_ALPHA_REGIONS
+ sqme_remn = (sqme_soft - sqme_cs)*log(xi_max)*xi_tilde
+ sqme = sqme0 - sqme_soft - sqme_coll + sqme_cs + sqme_remn
+ sqme = sqme * rsub%fks_kinematics%jac_rand (emitter)
+ end associate
+ !!! This occurs if the result is NaN.
+ if (sqme /= sqme) &
+ call write_computation_status ()
contains
subroutine write_computation_status ()
- print *, 'alr: ', alr
+ integer :: i_uborn
+ i_uborn = rsub%reg_data%regions(rsub%current_alr)%uborn_index
+ print *, 'alr: ', rsub%current_alr
print *, 'emitter: ', emitter
print *, 'xi_max: ', xi_max
print *, 'xi: ', xi, 'y: ', y
print *, 'sqme_born: ', rsub%sqme_born(i_uborn), &
'sqme_real: ', sqme0, &
'sqme_soft: ', sqme_soft, &
'sqme_coll: ', sqme_coll, &
'sqme_coll-soft: ', sqme_cs, &
'sqme_remn: ', sqme_remn
end subroutine write_computation_status
- end function real_subtraction_compute
-@ %def real_subtraction_compute
-@
+ end function real_subtraction_evaluate_region
+
+@ %def real_subtraction_evalute_region
+@
+<<real subtraction: real subtraction: TBP>>=
+ procedure :: get_limit_sqme => real_subtraction_get_limit_sqme
+<<real subtraction: procedures>>=
+ subroutine real_subtraction_get_limit_sqme (rsub, &
+ sqme_soft, sqme_coll, sqme_cs)
+ class(real_subtraction_t), intent(in) :: rsub
+ real(default), intent(out) :: sqme_soft, sqme_coll, sqme_cs
+ integer :: alr
+ alr = rsub%current_alr
+ sqme_soft = rsub%sub_soft%value(alr)
+ sqme_coll = rsub%sub_coll%value(alr)
+ sqme_cs = rsub%sub_coll%value_soft(alr)
+ end subroutine real_subtraction_get_limit_sqme
+
+@ %def real_subtraction_get_limit_sqme
+@
+<<real subtraction: real subtraction: TBP>>=
+ procedure :: get_phs_factor => real_subtraction_get_phs_factor
+<<real subtraction: procedures>>=
+ function real_subtraction_get_phs_factor (rsub) result (factor)
+ class(real_subtraction_t), intent(in) :: rsub
+ real(default) :: factor
+ real(default) :: s
+ s = (rsub%p_born(1)+rsub%p_born(2))**2
+ factor = s / (8*twopi2)
+ end function real_subtraction_get_phs_factor
+
+@ %def real_subtraction_get_phs_factor
+@
\subsection{xxx}
<<real subtraction: real subtraction: TBP>>=
procedure :: compute_sub_soft => real_subtraction_compute_sub_soft
<<real subtraction: procedures>>=
subroutine real_subtraction_compute_sub_soft &
- (rsub, p_born, p_real, fks_kinematics, alr, emitter, alpha_s)
+ (rsub, emitter, alpha_s)
class(real_subtraction_t), intent(inout) :: rsub
- type(vector4_t), intent(in), dimension(:) :: p_born, p_real
- type(real_kinematics_t), intent(in) :: fks_kinematics
- integer, intent(in) :: alr, emitter
+ integer, intent(in) :: emitter
real(default), intent(in) :: alpha_s
+ integer :: alr
+ alr = rsub%current_alr
associate (sregion => rsub%reg_data%regions(alr))
if (sregion%has_soft_divergence ()) then
- call rsub%sub_soft%compute (p_born, p_real, &
+ call rsub%sub_soft%compute (rsub%p_born, rsub%p_real, &
rsub%sqme_born(sregion%uborn_index), &
rsub%sqme_born_cc(:,:,sregion%uborn_index), &
- fks_kinematics%y(emitter), &
- fks_kinematics%y_soft(emitter), &
- fks_kinematics%phi, &
+ rsub%fks_kinematics%y(emitter), &
+ rsub%fks_kinematics%y_soft(emitter), &
+ rsub%fks_kinematics%phi, &
alpha_s, alr, emitter)
else
rsub%sub_soft%value(alr) = 0._default
end if
end associate
end subroutine real_subtraction_compute_sub_soft
@ %def real_subtraction_compute_sub_soft
@
<<real subtraction: real subtraction: TBP>>=
procedure :: compute_sub_coll => real_subtraction_compute_sub_coll
<<real subtraction: procedures>>=
- subroutine real_subtraction_compute_sub_coll (rsub, p_born, p_real, fks_kinematics, &
- alr, em, alpha_s)
+ subroutine real_subtraction_compute_sub_coll (rsub, em, alpha_s)
class(real_subtraction_t), intent(inout) :: rsub
- type(vector4_t), intent(in), dimension(:) :: p_born, p_real
- type(real_kinematics_t), intent(in) :: fks_kinematics
- integer, intent(in) :: alr, em
+ integer, intent(in) :: em
real(default), intent(in) :: alpha_s
real(default) :: xi
real(default) :: sqme_sc
complex(default) :: prod1, prod2
- xi = fks_kinematics%xi_tilde * fks_kinematics%xi_max (em)
+ integer :: alr
+ alr = rsub%current_alr
+ xi = rsub%fks_kinematics%xi_tilde * rsub%fks_kinematics%xi_max (em)
associate (sregion => rsub%reg_data%regions(alr))
if (sregion%has_collinear_divergence ()) then
if (rsub%sc_required(alr)) then
- call spinor_product (p_real(em), p_real(rsub%reg_data%nlegs_real), &
+ call spinor_product (rsub%p_real(em), rsub%p_real(rsub%reg_data%nlegs_real), &
prod1, prod2)
sqme_sc = real (prod1/prod2*rsub%me_sc)
else
sqme_sc = 0._default
end if
- call rsub%sub_coll%compute (sregion, p_born, &
+ call rsub%sub_coll%compute (sregion, rsub%p_born, &
rsub%sqme_born(sregion%uborn_index), &
sqme_sc, xi, alpha_s, alr)
else
rsub%sub_coll%value(alr) = 0._default
end if
end associate
end subroutine real_subtraction_compute_sub_coll
@ %def real_subtraction_compute_sub_coll
@
<<real subtraction: real subtraction: TBP>>=
procedure :: compute_sub_coll_soft => real_subtraction_compute_sub_coll_soft
<<real subtraction: procedures>>=
- subroutine real_subtraction_compute_sub_coll_soft (rsub, p_born, p_real, &
- alr, em, alpha_s)
+ subroutine real_subtraction_compute_sub_coll_soft (rsub, em, alpha_s)
class(real_subtraction_t), intent(inout) :: rsub
- type(vector4_t), intent(in), dimension(:) :: p_born, p_real
- integer, intent(in) :: alr, em
+ integer, intent(in) :: em
real(default), intent(in) :: alpha_s
complex(default) :: prod1, prod2
real(default) :: sqme_sc
real(default), parameter :: xi = 0
+ integer :: alr
+ alr = rsub%current_alr
associate (sregion => rsub%reg_data%regions(alr))
if (sregion%has_collinear_divergence ()) then
if (rsub%sc_required(alr)) then
- call spinor_product (p_real(em), p_real(rsub%reg_data%nlegs_real), &
+ call spinor_product (rsub%p_real(em), rsub%p_real(rsub%reg_data%nlegs_real), &
prod1, prod2)
sqme_sc = real (prod1/prod2*rsub%me_sc)
else
sqme_sc = 0._default
end if
- call rsub%sub_coll%compute_soft_limit (sregion, p_born, &
+ call rsub%sub_coll%compute_soft_limit (sregion, rsub%p_born, &
rsub%sqme_born(sregion%uborn_index), &
sqme_sc, xi, alpha_s, alr)
else
rsub%sub_coll%value_soft(alr) = 0._default
end if
end associate
end subroutine real_subtraction_compute_sub_coll_soft
@ %def real_subtraction_compute_sub_coll_soft
@
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Combining the FKS Pieces}
<<[[nlo_data.f90]]>>=
<<File header>>
module nlo_data
<<Use kinds>>
<<Use strings>>
use io_units
use constants, only: pi, twopi
use unit_tests
use diagnostics
use physics_defs
use process_constants !NODEP!
use sm_physics
use os_interface
use model_data
use parser
use pdg_arrays
use particle_specifiers
use phs_single
use state_matrices
use interactions
use lorentz
use prc_core
use pdg_arrays
use sf_base
use colors
use flavors
use fks_regions
use phs_fks
use virtual
use real_subtraction
<<Standard module head>>
<<nlo data: public>>
<<nlo data: types>>
contains
<<nlo data: procedures>>
end module nlo_data
@ %def nlo_data
@ We need to associate singular regions with compatible color flows.
<<nlo data: public>>=
public :: fks_template_t
<<nlo data: types>>=
type :: fks_template_t
integer :: mapping_type
logical :: count_kinematics = .false.
real(default) :: fks_dij_exp1
real(default) :: fks_dij_exp2
contains
<<nlo data: fks template: TBP>>
end type fks_template_t
@ %def fks_template_t
@
<<nlo data: fks template: TBP>>=
procedure :: write => fks_template_write
<<nlo data: procedures>>=
subroutine fks_template_write (object, unit)
class(fks_template_t), intent(in) :: object
integer, intent(in), optional :: unit
integer :: u
u = given_output_unit (unit)
write (u,'(1x,A)') 'FKS Template: '
write (u,'(1x,A,I0)') 'Mapping Type: ', object%mapping_type
write (u,'(1x,A,ES4.3,ES4.3)') 'd_ij exponentials: ', object%fks_dij_exp1, object%fks_dij_exp2
end subroutine fks_template_write
@ %def fks_template_write
@
<<nlo data: fks template: TBP>>=
procedure :: set_dij_exp => fks_template_set_dij_exp
<<nlo data: procedures>>=
subroutine fks_template_set_dij_exp (object, exp1, exp2)
class(fks_template_t), intent(inout) :: object
real(default), intent(in) :: exp1, exp2
object%fks_dij_exp1 = exp1
object%fks_dij_exp2 = exp2
end subroutine fks_template_set_dij_exp
@ %def fks_template_set_dij_exp
@
<<nlo data: fks template: TBP>>=
procedure :: set_mapping_type => fks_template_set_mapping_type
<<nlo data: procedures>>=
subroutine fks_template_set_mapping_type (object, val)
class(fks_template_t), intent(inout) :: object
integer, intent(in) :: val
object%mapping_type = val
end subroutine fks_template_set_mapping_type
@ %def fks_template_set_mapping_type
@
<<nlo data: fks template: TBP>>=
procedure :: set_counter => fks_template_set_counter
<<nlo data: procedures>>=
subroutine fks_template_set_counter (object)
class(fks_template_t), intent(inout) :: object
object%count_kinematics = .true.
end subroutine fks_template_set_counter
@ %def fks_template_set_counter
@
<<nlo data: types>>=
type :: ftuple_color_map_t
integer :: index
integer :: color_index
type(ftuple_color_map_t), pointer :: next
type(ftuple_color_map_t), pointer :: prev
contains
<<nlo data: color map: TBP>>
end type ftuple_color_map_t
@ %def ftuple_color_map_t
@
<<nlo data: color map: TBP>>=
procedure :: init => ftuple_color_map_init
<<nlo data: procedures>>=
subroutine ftuple_color_map_init (icm)
class(ftuple_color_map_t), intent(inout), target :: icm
icm%index = 0
icm%color_index = 0
nullify (icm%next)
nullify (icm%prev)
end subroutine ftuple_color_map_init
@ %def ftuple_color_map_init
@ Explain
<<nlo data: color map: TBP>>=
procedure :: present => ftuple_color_map_present
<<nlo data: procedures>>=
function ftuple_color_map_present (icm, color_index) result(pres)
class(ftuple_color_map_t), intent(in), target :: icm
integer, intent(in) :: color_index
logical :: pres
type(ftuple_color_map_t), pointer :: current
select type (icm)
type is (ftuple_color_map_t)
current => icm
pres = .false.
do
if (current%color_index == color_index) then
pres = .true.
exit
else
if (associated (current%next)) then
current => current%next
else
exit
end if
end if
end do
end select
end function ftuple_color_map_present
@ %def ftuple_color_map_present
@ Appends a color flow to the list
<<nlo data: color map: TBP>>=
procedure :: append => ftuple_color_map_append
<<nlo data: procedures>>=
subroutine ftuple_color_map_append (icm, val)
class(ftuple_color_map_t), intent(inout), target :: icm
integer, intent(in) :: val
type(ftuple_color_map_t), pointer :: current
select type (icm)
type is (ftuple_color_map_t)
if (.not. icm%present (val)) then
if (icm%index == 0) then
nullify(icm%next)
icm%index = 1
icm%color_index = val
else
current => icm
do
if (associated (current%next)) then
current => current%next
else
allocate (current%next)
nullify (current%next%next)
current%next%prev => current
current%next%index = current%index + 1
current%next%color_index = val
exit
end if
end do
end if
end if
end select
end subroutine ftuple_color_map_append
@ %def ftuple_color_map_append
@
<<nlo data: color map: TBP>>=
procedure :: get_n_entries => ftuple_color_map_get_n_entries
<<nlo data: procedures>>=
function ftuple_color_map_get_n_entries (icm) result(n_entries)
class(ftuple_color_map_t), intent(in), target :: icm
integer :: n_entries
type(ftuple_color_map_t), pointer :: current
select type (icm)
type is (ftuple_color_map_t)
current => icm
n_entries = 0
do
if (associated (current%next)) then
current => current%next
else
n_entries = current%index
exit
end if
end do
end select
end function ftuple_color_map_get_n_entries
@ %def ftuple_color_map_get_n_entries
@ Explain
<<nlo data: color map: TBP>>=
procedure :: get_index_array => ftuple_color_map_get_index_array
<<nlo data: procedures>>=
function ftuple_color_map_get_index_array (icm) result(iarr)
class(ftuple_color_map_t), intent(in), target :: icm
integer, dimension(:), allocatable :: iarr
type(ftuple_color_map_t), pointer :: current
integer :: n_entries
integer :: i
select type (icm)
type is (ftuple_color_map_t)
n_entries = icm%get_n_entries ()
allocate (iarr(n_entries))
do i = 1, n_entries
if (i == 1) then
current => icm
else
current => current%next
end if
iarr(i) = current%color_index
end do
end select
end function ftuple_color_map_get_index_array
@ %def ftuple_color_map_get_index_array
@
<<nlo data: color map: TBP>>=
procedure :: get_entry => ftuple_color_map_get_entry
<<nlo data: procedures>>=
function ftuple_color_map_get_entry (icm, index) result(entry)
class(ftuple_color_map_t), intent(in), target :: icm
integer, intent(in) :: index
integer :: entry
type(ftuple_color_map_t), pointer :: current
integer :: i
select type (icm)
type is (ftuple_color_map_t)
if (index <= icm%get_n_entries ()) then
do i = 1, icm%get_n_entries ()
if (i == 1) then
current => icm
else
current => current%next
end if
if (i == index) entry = current%color_index
end do
else
entry = 0
end if
end select
end function ftuple_color_map_get_entry
@ %def ftuple_color_map_get_entry
@ Explain
<<nlo data: color map: TBP>>=
procedure :: create_map => ftuple_color_map_create_map
<<nlo data: procedures>>=
recursive subroutine ftuple_color_map_create_map (icm, flst, &
emitter, allreg, color_states_born, color_states_real, p_rad_in)
class(ftuple_color_map_t), intent(inout) :: icm
type(flv_structure_t), intent(in) :: flst
integer, intent(in) :: emitter
type(ftuple_t), intent(in), dimension(:) :: allreg
integer, intent(in), dimension(:,:,:) :: color_states_born
integer, intent(in), dimension(:,:,:) :: color_states_real
integer, intent(in), optional :: p_rad_in
integer :: nreg, region
integer :: p1, p2, p_rad
integer :: flv_em, flv_rad
integer :: n_col_real, n_col_born
integer, dimension(2) :: col_em, col_rad
integer :: i
!!! splitting type: 1 - q -> qg
!!! 2 - g -> qq
!!! 3 - g -> gg
integer :: splitting_type_flv, splitting_type_col
nreg = size (allreg)
n_col_real = size (color_states_real (1,1,:))
n_col_born = size (color_states_born (1,1,:))
do region = 1, nreg
call allreg(region)%get (p1, p2)
if (p1 == emitter .or. p2 == emitter .or. present (p_rad_in)) then
if (.not. present (p_rad_in)) then
if (p1 == emitter) then
p_rad = p2
else
p_rad = p1
end if
else
p_rad = p_rad_in
end if
if (emitter /= 0) then
flv_em = flst%flst (emitter)
else
call icm%create_map &
(flst, 1, allreg, color_states_born, color_states_real, p_rad)
call icm%create_map &
(flst, 2, allreg, color_states_born, color_states_real, p_rad)
end if
flv_rad = flst%flst (p_rad)
if (is_quark (abs(flv_em)) .and. is_gluon (flv_rad)) then
splitting_type_flv = 1
else if (is_quark (abs(flv_em)) .and. flv_em + flv_rad == 0) then
splitting_type_flv = 2
else if (is_gluon (flv_em) .and. is_gluon (flv_rad)) then
splitting_type_flv = 3
else
splitting_type_flv = 0
end if
do i = 1, n_col_real
col_em = color_states_real(:,emitter,i)
col_rad = color_states_real(:,p_rad,i)
if (is_color_singlet (col_em(1), col_em(2)) &
.and. (is_color_doublet (col_rad(1), col_rad(2)) &
.or. is_color_ghost (col_rad(1), col_rad(2)))) then
splitting_type_col = 1
else if (is_color_singlet (col_em(1), col_em(2)) .and. &
is_color_singlet (col_rad(1), col_rad(2))) then
splitting_type_col = 2
else if (is_color_doublet (col_em(1), col_em(2)) .and. &
is_color_doublet (col_rad(1), col_rad(2))) then
splitting_type_col = 3
else
splitting_type_col = 0
end if
if (splitting_type_flv == splitting_type_col .and. &
splitting_type_flv /= 0) then
call icm%append (i)
end if
end do
end if
end do
contains
function is_color_singlet (c1, c2) result (singlet)
integer, intent(in) :: c1, c2
logical :: singlet
singlet = (c1 == 0 .and. c2 /= 0) .or. (c1 /= 0 .and. c2 == 0)
end function is_color_singlet
function is_color_doublet (c1, c2) result (doublet)
integer, intent(in) :: c1, c2
logical :: doublet
doublet = c1 /= 0 .and. c2 /= 0
end function is_color_doublet
function is_color_ghost (c1, c2) result (ghost)
integer, intent(in) :: c1, c2
logical :: ghost
ghost = c1 == 0 .and. c2 == 0
end function is_color_ghost
end subroutine ftuple_color_map_create_map
@ %def ftuple_color_map_create_map
@ This data type contains color information, necessary for both soft
and virtual counterterms.
<<nlo data: types>>=
type color_data_t
type(ftuple_color_map_t), dimension(:), allocatable :: icm
integer, dimension(:,:,:), allocatable :: col_state_born, col_state_real
logical, dimension(:,:), allocatable :: ghost_flag_born, ghost_flag_real
integer :: n_col_born, n_col_real
type(color_t), dimension(:,:), allocatable :: color_real, color_born
integer, dimension(:), allocatable :: col_born
complex(default), dimension(:), allocatable :: color_factors_born
integer, dimension(:,:), allocatable :: cf_index_real
real(default), dimension(:,:,:), allocatable :: beta_ij
logical :: color_is_conserved
contains
<<nlo data: color data: TBP>>
end type color_data_t
@ %def color_data_t
@
<<nlo data: color data: TBP>>=
procedure :: init => color_data_init
<<nlo data: procedures>>=
subroutine color_data_init (color_data, reg_data, prc_constants)
class(color_data_t), intent(inout) :: color_data
type(region_data_t), intent(inout) :: reg_data
type(process_constants_t), intent(inout), dimension(2) :: prc_constants
integer :: nlegs_born, nlegs_real
integer :: i
nlegs_born = reg_data%nlegs_born
nlegs_real = reg_data%nlegs_real
call prc_constants(1)%get_col_state (color_data%col_state_born)
call prc_constants(2)%get_col_state (color_data%col_state_real)
call prc_constants(2)%get_cf_index (color_data%cf_index_real)
call prc_constants(1)%get_color_factors (color_data%color_factors_born)
color_data%n_col_born = size (color_data%col_state_born(1,1,:))
color_data%n_col_real = size (color_data%col_state_real(1,1,:))
color_data%ghost_flag_born = prc_constants(1)%get_ghost_flag ()
color_data%ghost_flag_real = prc_constants(2)%get_ghost_flag ()
allocate (color_data%color_real (nlegs_real, color_data%n_col_real))
allocate (color_data%icm (reg_data%n_regions))
do i = 1, color_data%n_col_real
call color_init_from_array (color_data%color_real (:,i), &
color_data%col_state_real (:,:,i), &
color_data%ghost_flag_real (:,i))
end do
do i = 1, size(reg_data%regions)
call color_data%icm(i)%init
associate (region => reg_data%regions(i))
call color_data%icm(i)%create_map (region%flst_real, region%emitter, &
region%flst_allreg, color_data%col_state_born, &
color_data%col_state_real)
end associate
end do
call color_data%init_betaij (reg_data)
end subroutine color_data_init
@ %def color_data_init
@ Allocate and compute $\beta_{ij}$:
<<nlo data: color data: TBP>>=
procedure :: init_betaij => color_data_init_betaij
<<nlo data: procedures>>=
subroutine color_data_init_betaij (color_data, reg_data)
class(color_data_t), intent(inout) :: color_data
type(region_data_t), intent(inout) :: reg_data
integer :: i
allocate (color_data%beta_ij (reg_data%nlegs_born, &
reg_data%nlegs_born, reg_data%n_flv_born))
do i = 1, reg_data%n_flv_born
call color_data%fill_betaij_matrix (reg_data%nlegs_born, i, &
reg_data%regions(1)%flst_real, reg_data)
end do
end subroutine color_data_init_betaij
@ %def color_data_init_betaij
@ Actual computation of $\beta_{ij}$.
<<nlo data: color data: TBP>>=
procedure :: fill_betaij_matrix => color_data_fill_betaij_matrix
<<nlo data: procedures>>=
subroutine color_data_fill_betaij_matrix &
(color_data, n_legs, uborn_index, flst_real, reg_data)
class(color_data_t), intent(inout) :: color_data
integer, intent(in) :: n_legs, uborn_index
type(flv_structure_t), intent(in) :: flst_real
type(region_data_t), intent(inout) :: reg_data
integer :: em1, em2
associate (flv_born => reg_data%flv_born (uborn_index))
do em1 = 1, n_legs
do em2 = 1, n_legs
if (flv_born%colored(em1) .and. flv_born%colored(em2)) then
if (em1 < em2) then
color_data%beta_ij (em1, em2, uborn_index) &
= color_data%compute_bij &
(reg_data, uborn_index, flst_real, em1, em2)
else if (em1 > em2) then
!!! B_ij is symmetric
color_data%beta_ij (em1, em2, uborn_index) = &
color_data%beta_ij (em2, em1, uborn_index)
else
if (is_quark (abs (flv_born%flst (em1)))) then
color_data%beta_ij (em1, em2, uborn_index) = -cf
else
color_data%beta_ij (em1, em2, uborn_index) = -ca
end if
end if
else
color_data%beta_ij (em1, em2, uborn_index) = 0.0
end if
end do
end do
end associate
call check_color_conservation (color_data%beta_ij (:,:,uborn_index), &
n_legs, color_data%color_is_conserved)
contains
subroutine check_color_conservation (bij_matrix, n_legs, success)
real(default), intent(in), dimension(:,:) :: bij_matrix
integer, intent(in) :: n_legs
logical, intent(out) :: success
logical, dimension(:), allocatable :: check
integer :: i, j
real(default) :: bcheck
real(default), parameter :: tol = 0.0001_default
allocate (check (n_legs))
do i = 1, n_legs
bcheck = 0.0
do j = 1, n_legs
if (i /= j) bcheck = bcheck + bij_matrix (i, j)
end do
if (is_quark (abs(flst_real%flst (i))) .or. &
is_gluon (flst_real%flst (i))) then
if (is_quark (abs(flst_real%flst (i))) .and. &
(bcheck - cf) < tol) then
check (i) = .true.
else if (is_gluon (flst_real%flst (i)) .and. &
(bcheck - ca) < tol) then
check (i) = .true.
else
check (i) = .false.
end if
else
if (bcheck < tol) then
check (i) = .true.
else
check (i) = .false.
end if
end if
end do
if (.not. all (check)) then
success = .false.
! call msg_fatal ("Color conservation violated!")
else
success = .true.
end if
end subroutine check_color_conservation
end subroutine color_data_fill_betaij_matrix
@ %def color_data_fill_betaij_matrix
@ Explain
<<nlo data: color data: TBP>>=
procedure :: compute_bij => color_data_compute_bij
<<nlo data: procedures>>=
function color_data_compute_bij &
(color_data, reg_data, uborn_index, flst_real, em1, em2) result (bij)
class(color_data_t), intent(inout) :: color_data
type(region_data_t), intent(inout) :: reg_data
integer, intent(in) :: uborn_index
type(flv_structure_t), intent(in) :: flst_real
integer, intent(in) :: em1, em2
real(default) :: bij
logical, dimension(:,:), allocatable :: cf_present
type(singular_region_t), dimension(2,100) :: reg
integer :: i, j, k, l
type(ftuple_color_map_t) :: icm1, icm2
integer :: i1, i2
real(default) :: color_factor, color_factor_born
integer, dimension(2) :: i_reg
logical , dimension(2) :: found
integer, dimension(2,100) :: map_em_col_tmp
integer, dimension(:), allocatable :: map_em_col1, map_em_col2
integer, dimension(2) :: col1, col2
integer, dimension(:), allocatable :: iarray1, iarray2
integer, dimension(:), allocatable :: iisec1, iisec2
integer :: sign
color_factor = 0.0; color_factor_born = 0.0
found = .false.
!!! Include distinction between Born flavors
do i = 1, size (color_data%color_factors_born)
color_factor_born = color_factor_born + color_data%color_factors_born (i)
end do
i1 = 1
i2 = 1
!!! Catch case em = 0
if (em1 == 0 .or. em2 == 0) then
!!! What to do?
bij = 0.0
else
do i = 1, color_data%n_col_real
col1 = color_data%col_state_real (:, em1, i)
col2 = color_data%col_state_real (:, reg_data%nlegs_real, i)
if (share_line (col1, col2)) then
map_em_col_tmp(1,i1) = i
i1 = i1+1
end if
col1 = color_data%col_state_real (:, em2, i)
if (share_line (col1, col2)) then
map_em_col_tmp(2,i2) = i
i2 = i2 + 1
end if
end do
allocate (map_em_col1 (i1), map_em_col2 (i2))
map_em_col1 = map_em_col_tmp (1,1:i1-1)
map_em_col2 = map_em_col_tmp (2,1:i2-1)
i_reg = 1
do i = 1, reg_data%n_regions
if (uborn_index == reg_data%regions(i)%uborn_index) then
if (em1 == reg_data%regions(i)%emitter .or. &
(em1 <= 2 .and. reg_data%regions(i)%emitter == 0)) then
reg(1,i_reg(1)) = reg_data%regions(i)
i_reg(1) = i_reg(1)+1
found(1) = .true.
end if
if (em2 == reg_data%regions(i)%emitter .or. &
(em2 <= 2 .and. reg_data%regions(i)%emitter == 0)) then
reg(2,i_reg(2)) = reg_data%regions(i)
i_reg(2) = i_reg(2)+1
found(2) = .true.
end if
end if
end do
if (.not. (found(1).and.found(2))) then
bij = 0
return
end if
do i = 1, i_reg(1)-1
do j = 1, i_reg(2)-1
icm1 = color_data%icm (reg(1,i)%alr)
icm2 = color_data%icm (reg(2,j)%alr)
iarray1 = icm1%get_index_array ()
iarray2 = icm2%get_index_array ()
iisec1 = pack (iarray1, [ (any(iarray1(i) == map_em_col1), &
i = 1, size(iarray1)) ])
iisec2 = pack (iarray2, [ (any(iarray2(i) == map_em_col2), &
i = 1, size(iarray2)) ])
cf_present = color_index_present (color_data%cf_index_real)
do k = 1, size (iisec1)
do l = 1, size (iisec2)
i1 = iisec1(k)
i2 = iisec2(l)
if (cf_present (i1, i2)) then
if (is_gluon (flst_real%flst (em1)) .or. &
is_gluon (flst_real%flst (em2))) then
sign = get_sign (color_data%col_state_real (:,:,i1)) * &
get_sign (color_data%col_state_real (:,:,i2))
else
sign = 1
end if
color_factor = color_factor + sign*compute_color_factor &
(color_data%color_real(:,i1), &
color_data%color_real(:,i2))
end if
end do
end do
end do
end do
!!! The real color factor always differs from the Born one
!!! by one vertex factor. Thus, apply the factor 1/2
bij = color_factor / (2 * color_factor_born)
end if
contains
function share_line (col1, col2) result (share)
integer, intent(in), dimension(2) :: col1, col2
logical :: share
logical :: id1, id2, id3
id1 = (abs(col1(1)) == abs(col2(1)) .and. col1(1) /= 0) .or. &
(abs(col1(2)) == abs(col2(2)) .and. col1(2) /= 0)
id2 = (abs(col1(1)) == abs(col2(2)) .and. col1(1) /= 0) .or. &
(abs(col1(2)) == abs(col2(1)) .and. col1(2) /= 0)
id3 = col2(1) == 0 .and. col2(2) == 0
if (id1 .or. id2 .or. id3) then
share = .true.
else
share = .false.
end if
end function share_line
function get_sign (col) result (sign)
integer, intent(in), dimension(:,:) :: col
integer :: sign
integer, dimension(:), allocatable :: iref, iperm
integer :: iref1, iperm1
integer :: n, i, i_first, j
integer :: i1, i2
integer :: p1, p2
p1 = 2; p2 = 2
iref1 = 0; iperm1 = 0; i_first = 0
do i = 1, size(col(1,:))
if (.not. all (col(:,i) == 0)) then
if (col(1,i) == 0) then
i1 = col(2,i)
iref1 = i; iperm1 = i
i_first = i
else
i1 = col(1,i)
iref1 = i; iperm1 = i
i_first = i
end if
exit
end if
end do
if (iref1 == 0 .or. iperm1 == 0 .or. i_first == 0) &
call msg_fatal ("Invalid color structure")
n = size(col(1,:)) - i_first + 1
allocate (iref(n), iperm(n))
iref(1) = iref1; iperm(1) = iperm1
do i = i_first+1, size(col(1,:))
if (all (col(:,i) == 0)) cycle
if (i == size(col(1,:))) then
iref(p1) = i_first + 1
else
iref(p1) = i + 1
p1 = p1 + 1
end if
do j = i_first+1, size(col(1,:))
if (col(1,j) == -i1) then
i1 = col(2,j)
iperm(p2) = j
p2 = p2 + 1
exit
else if (col(2,j) == -i1) then
i1 = col(1,j)
iperm(p2) = j
p2 = p2 + 1
exit
end if
end do
end do
sign = 1
do i = 1, n
if (iperm(i) == iref(i)) then
cycle
else
do j = i+1, n
if (iperm(j) == iref(i)) then
i1 = j
exit
end if
end do
i2 = iperm(i)
iperm(i) = iperm(i1)
iperm(i1) = i2
sign = -sign
end if
end do
end function get_sign
function color_index_present (cf_index) result (cf_present)
integer, intent(in), dimension(:,:), allocatable :: cf_index
logical, dimension(:,:), allocatable :: cf_present
integer :: n_col
integer :: c, i1, i2
n_col = size (cf_index(1,:))
allocate (cf_present (n_col, n_col))
cf_present = .false.
do c = 1, n_col
i1 = cf_index (1, c)
i2 = cf_index (2, c)
cf_present (i1, i2) = .true.
if (i1 /= i2) cf_present(i2, i1) = .true.
end do
end function color_index_present
end function color_data_compute_bij
@ %def color_data_compute_bij
@
<<nlo data: color data: TBP>>=
procedure :: write => color_data_write
<<nlo data: procedures>>=
subroutine color_data_write (color_data, unit)
class(color_data_t), intent(in) :: color_data
integer, intent(in), optional :: unit
integer :: u, i, i1, i2
integer :: n_legs
u = given_output_unit (unit); if (u < 0) return
n_legs = size (color_data%beta_ij, dim=2)
write (u, "(1x,A)") "Color information: "
write (u, "(1x,A,1x,I1)") "Number of Born color states: ", &
color_data%n_col_born
write (u, "(1x,A,1x,I1)") "Number of real color states: ", &
color_data%n_col_real
write (u, "(1x,A)") "Color correlation: "
do i = 1, size (color_data%beta_ij, dim=3)
write (u, "(1x,A,1x,I1)") "State nr. ", i
write (u, "(1x,A)") "-------------"
write (u, "(1x,A,1x,A,1x,A)") "i1", "i2", "color factor"
do i1 = 1, n_legs
do i2 = 1, i1
write (u, "(1x,I1,1x,I1,1x,F5.2)") &
i1, i2, color_data%beta_ij (i1,i2,i)
end do
end do
write (u, "(1x,A)") "========================================"
end do
if (color_data%color_is_conserved) then
write (u, "(1x,A)") "Color is conserved."
else
write (u, "(1x,A)") "Fatal error: Color conversation is violated."
end if
end subroutine color_data_write
@ %def color_data_write
@
<<nlo data: nlo data: TBP>>=
procedure :: compute_sqme_real_fin => nlo_data_compute_sqme_real_fin
<<nlo data: procedures>>=
function nlo_data_compute_sqme_real_fin (nlo_data, weight, p_real) result (sqme_fin)
class(nlo_data_t), intent(inout) :: nlo_data
real(default), intent(in) :: weight
type(vector4_t), intent(inout), dimension(:), allocatable :: p_real
type(vector4_t), dimension(:), allocatable :: p_born
real(default) :: sqme_fin
integer :: emitter
if (.not. nlo_data%alpha_s_born_set) &
call msg_fatal ("Strong coupling not set for real calculation - abort")
emitter = nlo_data%get_active_emitter ()
p_born = interaction_get_momenta (nlo_data%int_born)
- sqme_fin = nlo_data%real_terms%compute (emitter, p_real, p_born, &
- nlo_data%real_kinematics, &
+ call nlo_data%real_terms%set_momenta (p_born, p_real)
+ call nlo_data%real_terms%set_fks_kinematics (nlo_data%real_kinematics)
+ sqme_fin = nlo_data%real_terms%compute (emitter, &
nlo_data%alpha_s_born)
sqme_fin = sqme_fin * weight
end function nlo_data_compute_sqme_real_fin
@ %def nlo_data_compute_sqme_real_fin
@ Check if there are massive emitters. Since the mass-structure of all
underlying Born configurations have to be the same, we just use the first
one to determine this.
<<nlo data: nlo data: TBP>>=
procedure :: has_massive_emitter => nlo_data_has_massive_emitter
<<nlo data: procedures>>=
function nlo_data_has_massive_emitter (nlo_data) result (val)
class(nlo_data_t), intent(inout) :: nlo_data
logical :: val
integer :: n_tot, i
val = .false.
n_tot = nlo_data%n_in + nlo_data%n_out_born
do i = nlo_data%n_in+1, n_tot
if (any (i == nlo_data%reg_data%emitters)) &
val = val .or. nlo_data%reg_data%flv_born(1)%massive(i)
end do
end function nlo_data_has_massive_emitter
@ %def nlo_data_has_massive_emitter
@
\subsection{Putting it together}
This data type governs the whole calculation. It contains information
about color, spin and flavor as well as the information about the Born
process.
<<nlo data: public>>=
public :: nlo_data_t
<<nlo data: types>>=
type :: nlo_data_t
type(string_t) :: nlo_type
type(region_data_t) :: reg_data
integer :: n_alr
integer :: n_in
integer :: n_out_born
integer :: n_out_real
integer :: n_allowed_born
integer :: n_flv_born
integer :: n_flv_real
integer :: active_emitter
complex(default), dimension(:), allocatable :: amp_born
type(color_data_t) :: color_data
type(real_kinematics_t) :: real_kinematics
type(virtual_t) :: virtual_terms
type(real_subtraction_t) :: real_terms
integer, dimension(:,:), allocatable :: flv_state_born
integer, dimension(:,:), allocatable :: flv_state_real
integer, dimension(:), allocatable :: flv_born
integer, dimension(:), allocatable :: hel_born
integer, dimension(:), allocatable :: col_born
real(default) :: alpha_s_born
logical :: alpha_s_born_set
real(default), dimension(:), allocatable :: sqme_born
real(default), dimension(:), allocatable :: sqme_real_non_sub
real(default), dimension(:,:,:), allocatable :: sqme_born_cc
complex(default) :: me_sc
real(default), public :: sqme_real
real(default), public :: sqme_virt
type(interaction_t), public :: int_born
type(kinematics_counter_t), public :: counter
logical, public :: counter_exists = .false.
logical :: use_internal_color_correlations = .true.
logical :: use_internal_spin_correlations = .false.
contains
<<nlo data: nlo data: TBP>>
end type nlo_data_t
@ %def nlo_data_t
@
<<nlo data: public>>=
public :: nlo_pointer_t
<<nlo data: types>>=
type :: nlo_pointer_t
type(nlo_data_t), public, pointer :: nlo_data => null ()
end type nlo_pointer_t
@ %def nlo_pointer_t
@
<<nlo data: nlo data: TBP>>=
procedure :: init => nlo_data_init
<<nlo data: procedures>>=
subroutine nlo_data_init (nlo_data, nlo_type, prc_constants, template, model)
class(nlo_data_t), intent(inout) :: nlo_data
type(string_t), intent(in) :: nlo_type
type(process_constants_t), intent(inout), dimension(2) :: prc_constants
type(fks_template_t), intent(in) :: template
type(model_data_t), intent(in) :: model
integer :: i
nlo_data%nlo_type = nlo_type
select case (char (nlo_type))
case ('Real')
nlo_data%flv_state_born = prc_constants(1)%get_flv_state ()
nlo_data%flv_state_real = prc_constants(2)%get_flv_state ()
call nlo_data%reg_data%init (model, &
nlo_data%flv_state_born, &
nlo_data%flv_state_real, &
template%mapping_type)
select type (mapping => nlo_data%reg_data%fks_mapping)
type is (fks_mapping_default_t)
call mapping%set_parameter (template%fks_dij_exp1, template%fks_dij_exp2)
end select
nlo_data%n_flv_born = nlo_data%reg_data%n_flv_born
nlo_data%n_flv_real = nlo_data%reg_data%n_flv_real
allocate (nlo_data%sqme_born (nlo_data%n_flv_born))
allocate (nlo_data%sqme_real_non_sub (nlo_data%n_flv_real))
nlo_data%sqme_born = 0.0
nlo_data%n_in = prc_constants(2)%n_in
nlo_data%n_out_born = prc_constants(1)%n_out
nlo_data%n_out_real = prc_constants(2)%n_out
allocate (nlo_data%sqme_born_cc (nlo_data%n_in + nlo_data%n_out_born, &
nlo_data%n_in + nlo_data%n_out_born, &
nlo_data%n_flv_born))
if (nlo_data%use_internal_color_correlations) &
call nlo_data%color_data%init (nlo_data%reg_data, prc_constants)
nlo_data%alpha_s_born_set = .false.
call nlo_data%init_real_kinematics
call nlo_data%real_terms%init (nlo_data%reg_data, &
nlo_data%n_in + nlo_data%n_out_born, &
nlo_data%n_in + nlo_data%n_out_born)
nlo_data%counter_exists = template%count_kinematics
if (nlo_data%counter_exists) call nlo_data%counter%init(20)
end select
end subroutine nlo_data_init
@ %def nlo_data_init
@
<<nlo data: nlo data: TBP>>=
procedure :: copy_matrix_elements => nlo_data_copy_matrix_elements
<<nlo data: procedures>>=
subroutine nlo_data_copy_matrix_elements (nlo_data)
class(nlo_data_t), intent(inout) :: nlo_data
nlo_data%real_terms%sqme_born = nlo_data%sqme_born
nlo_data%real_terms%sqme_born_cc = nlo_data%sqme_born_cc
nlo_data%real_terms%sqme_real_bare = nlo_data%sqme_real_non_sub
nlo_data%real_terms%me_sc = nlo_data%me_sc
end subroutine nlo_data_copy_matrix_elements
@ %def nlo_data_copy_matrix_elements
@
<<nlo data: nlo data: TBP>>=
procedure :: write => nlo_data_write
<<nlo data: procedures>>=
!!! Incomplete
subroutine nlo_data_write (nlo_data, unit)
class(nlo_data_t), intent(in) :: nlo_data
integer, intent(in), optional :: unit
integer :: i, u
u = given_output_unit (unit); if (u < 0) return
write (u, "(1x,A,I0)") "n_alr = ", nlo_data%n_alr
write (u, "(1x,A,I0)") "n_in = ", nlo_data%n_in
write (u, "(1x,A,I0)") "n_out_born = ", nlo_data%n_out_born
write (u, "(1x,A,I0)") "n_out_real = ", nlo_data%n_out_real
write (u, "(1x,A,I0)") "n_allowed_born = ", nlo_data%n_allowed_born
end subroutine nlo_data_write
@ %def nlo_data_write
@
<<nlo data: nlo data: TBP>>=
procedure :: init_born_amps => nlo_data_init_born_amps
<<nlo data: procedures>>=
subroutine nlo_data_init_born_amps (nlo_data, n, internal_correlations)
class(nlo_data_t), intent(inout) :: nlo_data
integer, intent(in) :: n
logical, intent(in) :: internal_correlations
nlo_data%n_allowed_born = n
allocate (nlo_data%amp_born (n))
end subroutine nlo_data_init_born_amps
@ %def nlo_data_init_born_amps
@
<<nlo data: nlo data: TBP>>=
procedure :: set_internal_procedures => nlo_data_set_internal_procedures
<<nlo data: procedures>>=
subroutine nlo_data_set_internal_procedures (nlo_data, flag_color, flag_spin)
class(nlo_data_t), intent(inout) :: nlo_data
logical, intent(in) :: flag_color, flag_spin
nlo_data%use_internal_color_correlations = flag_color
nlo_data%real_terms%sub_soft%use_internal_color_correlations = flag_color
nlo_data%virtual_terms%use_internal_color_correlations = flag_color
nlo_data%use_internal_spin_correlations = flag_spin
end subroutine nlo_data_set_internal_procedures
@ %def nlo_data_set_internal_procedures
@
<<nlo data: nlo data: TBP>>=
procedure :: init_virtual => nlo_data_init_virtual
<<nlo data: procedures>>=
subroutine nlo_data_init_virtual (nlo_data)
class(nlo_data_t), intent(inout) :: nlo_data
call nlo_data%virtual_terms%init (nlo_data%flv_state_born)
end subroutine nlo_data_init_virtual
@ %def nlo_data_init_virtual
@
<<nlo data: nlo data: TBP>>=
procedure :: set_nlo_type => nlo_data_set_nlo_type
<<nlo data: procedures>>=
subroutine nlo_data_set_nlo_type (nlo_data, nlo_type)
class(nlo_data_t), intent(inout) :: nlo_data
type(string_t), intent(in) :: nlo_type
nlo_data%nlo_type = nlo_type
end subroutine nlo_data_set_nlo_type
@ %def nlo_data_set_nlo_type
@
<<nlo data: nlo data: TBP>>=
procedure :: get_nlo_type => nlo_data_get_nlo_type
<<nlo data: procedures>>=
pure function nlo_data_get_nlo_type (nlo_data) result(nlo_type)
class(nlo_data_t), intent(in) :: nlo_data
type(string_t) :: nlo_type
nlo_type = nlo_data%nlo_type
end function nlo_data_get_nlo_type
@ %def nlo_data_get_nlo_type
@
<<nlo data: nlo data: TBP>>=
procedure :: get_emitters => nlo_data_get_emitters
<<nlo data: procedures>>=
function nlo_data_get_emitters (nlo_data) result(emitters)
class(nlo_data_t), intent(inout) :: nlo_data
integer, dimension(:), allocatable :: emitters
emitters = nlo_data%reg_data%get_emitters ()
end function nlo_data_get_emitters
@ %def nlo_data_get_emmiters
@
<<nlo data: nlo data: TBP>>=
procedure :: set_active_emitter => nlo_data_set_active_emitter
<<nlo data: procedures>>=
subroutine nlo_data_set_active_emitter (nlo_data, emitter)
class(nlo_data_t), intent(inout) :: nlo_data
integer, intent(in) :: emitter
nlo_data%active_emitter = emitter
end subroutine nlo_data_set_active_emitter
@ %def nlo_data_set_active_emitter
@
<<nlo data: nlo data: TBP>>=
procedure :: get_active_emitter => nlo_data_get_active_emitter
<<nlo data: procedures>>=
function nlo_data_get_active_emitter (nlo_data) result(emitter)
class(nlo_data_t), intent(inout) :: nlo_data
integer :: emitter
emitter = nlo_data%active_emitter
end function nlo_data_get_active_emitter
@ %def nlo_data_get_active_emitter
@
<<nlo data: nlo data: TBP>>=
procedure :: set_flv_born => nlo_data_set_flv_born
<<nlo data: procedures>>=
subroutine nlo_data_set_flv_born (nlo_data, flv_in)
class(nlo_data_t), intent(inout) :: nlo_data
integer, intent(in), dimension(:), allocatable :: flv_in
allocate (nlo_data%flv_born (size (flv_in)))
nlo_data%flv_born = flv_in
end subroutine nlo_data_set_flv_born
@ %def nlo_data_set_flv_born
@
<<nlo data: nlo data: TBP>>=
procedure :: set_hel_born => nlo_data_set_hel_born
<<nlo data: procedures>>=
subroutine nlo_data_set_hel_born (nlo_data, hel_in)
class(nlo_data_t), intent(inout) :: nlo_data
integer, intent(in), dimension(:), allocatable :: hel_in
allocate (nlo_data%hel_born (size (hel_in)))
nlo_data%hel_born = hel_in
end subroutine nlo_data_set_hel_born
@ %def nlo_data_set_hel_born
@
<<nlo data: nlo data: TBP>>=
procedure :: set_col_born => nlo_data_set_col_born
<<nlo data: procedures>>=
subroutine nlo_data_set_col_born (nlo_data, col_in)
class(nlo_data_t), intent(inout) :: nlo_data
integer, intent(in), dimension(:), allocatable :: col_in
allocate (nlo_data%col_born (size (col_in)))
nlo_data%col_born = col_in
end subroutine nlo_data_set_col_born
@ %def nlo_data_set_col_born
@
<<nlo data: nlo data: TBP>>=
procedure :: set_alpha_s_born => nlo_data_set_alpha_s_born
<<nlo data: procedures>>=
subroutine nlo_data_set_alpha_s_born (nlo_data, as_born)
class (nlo_data_t), intent(inout) :: nlo_data
real(default), intent(in) :: as_born
nlo_data%alpha_s_born = as_born
nlo_data%alpha_s_born_set = .true.
end subroutine nlo_data_set_alpha_s_born
@ %def nlo_data_set_alpha_s_born
@
<<nlo data: nlo data: TBP>>=
procedure :: init_real_kinematics => nlo_data_init_real_kinematics
<<nlo data: procedures>>=
subroutine nlo_data_init_real_kinematics (nlo_data)
class(nlo_data_t), intent(inout) :: nlo_data
integer :: n_tot
n_tot = nlo_data%n_in + nlo_data%n_out_born
allocate (nlo_data%real_kinematics%xi_max (n_tot))
allocate (nlo_data%real_kinematics%y (n_tot))
allocate (nlo_data%real_kinematics%y_soft (n_tot))
allocate (nlo_data%real_kinematics%jac_rand (n_tot))
nlo_data%real_kinematics%xi_tilde = 0
nlo_data%real_kinematics%y = 0
nlo_data%real_kinematics%xi_max = 0
nlo_data%real_kinematics%phi = 0
end subroutine nlo_data_init_real_kinematics
@ %def nlo_data_init_real_kinematics
@
<<nlo data: nlo data: TBP>>=
procedure :: set_real_kinematics => nlo_data_set_real_kinematics
<<nlo data: procedures>>=
subroutine nlo_data_set_real_kinematics (nlo_data, xi_tilde, y, phi, xi_max, &
jac, jac_rand)
class(nlo_data_t), intent(inout) :: nlo_data
real(default), dimension(:), allocatable :: xi_max, y
real(default), intent(in) :: xi_tilde
real(default), intent(in) :: phi
real(default), intent(in), dimension(3) :: jac
real(default), intent(in), dimension(:), allocatable :: jac_rand
nlo_data%real_kinematics%xi_tilde = xi_tilde
nlo_data%real_kinematics%y = y
nlo_data%real_kinematics%phi = phi
nlo_data%real_kinematics%xi_max = xi_max
nlo_data%real_kinematics%jac = jac
nlo_data%real_kinematics%jac_rand = jac_rand
end subroutine nlo_data_set_real_kinematics
@ %def nlo_data_set_real_kinematics
@
<<nlo data: nlo data: TBP>>=
procedure :: get_real_kinematics => nlo_data_get_real_kinematics
<<nlo data: procedures>>=
subroutine nlo_data_get_real_kinematics &
(nlo_data, em, xi_tilde, y, xi_max, jac, phi, jac_rand)
class(nlo_data_t), intent(in) :: nlo_data
integer, intent(in) :: em
real(default), intent(out) :: xi_tilde, y, xi_max
real(default), intent(out), dimension(3), optional :: jac
!!! For most applications, phi is not relevant. Thus, it is not
!!! always transferred as a dummy-variable
real(default), intent(out), optional :: phi
real(default), intent(out), dimension(:), optional :: jac_rand
xi_tilde = nlo_data%real_kinematics%xi_tilde
y = nlo_data%real_kinematics%y(em)
xi_max = nlo_data%real_kinematics%xi_max (em)
if (present (jac)) jac = nlo_data%real_kinematics%jac
if (present (phi)) phi = nlo_data%real_kinematics%phi
if (present (jac_rand)) jac_rand = nlo_data%real_kinematics%jac_rand
end subroutine nlo_data_get_real_kinematics
@ %def nlo_data_get_real_kinematics
@
<<nlo data: nlo data: TBP>>=
procedure :: set_real_momenta => nlo_data_set_real_momenta
<<nlo data: procedures>>=
subroutine nlo_data_set_real_momenta (nlo_data, p)
class(nlo_data_t), intent(inout) :: nlo_data
type(vector4_t), intent(in), dimension(:), allocatable :: p
nlo_data%real_kinematics%p_real = p
end subroutine nlo_data_set_real_momenta
@ %def nlo_data_set_real_momenta
<<nlo data: nlo data: TBP>>=
procedure :: get_real_momenta => nlo_data_get_real_momenta
<<nlo data: procedures>>=
function nlo_data_get_real_momenta (nlo_data) result (p)
class(nlo_data_t), intent(inout) :: nlo_data
type(vector4_t), dimension(:), allocatable :: p
p = nlo_data%real_kinematics%p_real
end function nlo_data_get_real_momenta
@ %def nlo_data_get_real_kinematics
@
<<nlo data: nlo data: TBP>>=
procedure :: set_y_soft => nlo_data_set_y_soft
<<nlo data: procedures>>=
subroutine nlo_data_set_y_soft (nlo_data, y_soft, emitter)
class(nlo_data_t), intent(inout) :: nlo_data
real(default), intent(in) :: y_soft
integer, intent(in) :: emitter
nlo_data%real_kinematics%y_soft(emitter) = y_soft
end subroutine nlo_data_set_y_soft
@ %def nlo_data_set_y_soft
@
<<nlo data: nlo data: TBP>>=
procedure :: set_jacobian => nlo_data_set_jacobian
<<nlo data: procedures>>=
subroutine nlo_data_set_jacobian (nlo_data, jac)
class(nlo_data_t), intent(inout) :: nlo_data
real(default), intent(in), dimension(3) :: jac
nlo_data%real_kinematics%jac = jac
end subroutine nlo_data_set_jacobian
@ %def nlo_data_set_jacobian
@
\subsection{xxx}
<<nlo data: nlo data: TBP>>=
procedure :: compute_virt => nlo_data_compute_virt
<<nlo data: procedures>>=
function nlo_data_compute_virt (nlo_data, i_flv, int_born) result(sqme_virt)
class(nlo_data_t), intent(inout) :: nlo_data
integer, intent(in) :: i_flv
type(interaction_t), intent(in) :: int_born
real(default) :: sqme_virt
type(vector4_t), dimension(:), allocatable :: p_born
p_born = interaction_get_momenta (int_born)
if (nlo_data%use_internal_color_correlations) then
call nlo_data%virtual_terms%evaluate (nlo_data%reg_data, &
i_flv, &
nlo_data%alpha_s_born, &
p_born, &
nlo_data%sqme_born(i_flv), &
nlo_data%color_data%beta_ij)
else
call nlo_data%virtual_terms%evaluate (nlo_data%reg_data, &
i_flv, &
nlo_data%alpha_s_born, &
p_born, &
nlo_data%sqme_born(i_flv), &
nlo_data%sqme_born_cc)
end if
sqme_virt = nlo_data%virtual_terms%sqme_virt
end function nlo_data_compute_virt
@ %def nlo_data_compute_virt
@
<<nlo data: nlo data: TBP>>=
procedure :: requires_spin_correlation => &
nlo_data_requires_spin_correlation
<<nlo data: procedures>>=
function nlo_data_requires_spin_correlation (nlo_data, i_flv) result (val)
class(nlo_data_t), intent(in) :: nlo_data
integer, intent(in) :: i_flv
logical :: val
val = nlo_data%real_terms%sc_required (i_flv)
end function nlo_data_requires_spin_correlation
@ %def nlo_data_requires_spin_correlation
@
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The FKS phase space}
To handle the FKS phase space and adapt it to multi-channel
integration, some extensions have to be made for pre-existing types.
The \texttt{phs\_tree}-type has an attribute containing the tree-code
of the emitting branch.
<<PHS trees: phs tree: TBP>>=
procedure :: get_emitter => phs_tree_get_emitter
<<PHS trees: procedures>>=
function phs_tree_get_emitter (tree) result (emitter)
class(phs_tree_t), intent(in) :: tree
integer :: emitter
emitter = tree%emitter
end function phs_tree_get_emitter
@ %def phs_tree_get_emitter
@
<<[[phs_fks.f90]]>>=
<<File header>>
module phs_fks
<<Use kinds>>
<<Use strings>>
use constants
use diagnostics
use lorentz
use flavors
use sf_mappings
use sf_base
use phs_base
use phs_wood
use process_constants
<<Standard module head>>
<<phs fks: public>>
<<phs fks: types>>
<<phs fks: interfaces>>
contains
<<phs fks: procedures>>
end module phs_fks
@ %def phs_fks
@
<<phs fks: public>>=
public :: phs_fks_config_t
<<phs fks: types>>=
type, extends (phs_wood_config_t) :: phs_fks_config_t
contains
<<phs fks: fks config: TBP>>
end type phs_fks_config_t
@ %def phs_fks_config_t
<<phs fks: fks config: TBP>>=
procedure :: final => phs_fks_config_final
<<phs fks: procedures>>=
subroutine phs_fks_config_final (object)
class(phs_fks_config_t), intent(inout) :: object
! call object%phs_wood_config_t%final ()
end subroutine phs_fks_config_final
@ %def phs_fks_config_final
@
<<phs fks: fks config: TBP>>=
procedure :: write => phs_fks_config_write
<<phs fks: procedures>>=
subroutine phs_fks_config_write (object, unit)
class(phs_fks_config_t), intent(in) :: object
integer, intent(in), optional :: unit
call object%phs_wood_config_t%write
end subroutine phs_fks_config_write
@ %def phs_fks_config_write
@
<<phs fks: fks config: TBP>>=
procedure :: configure => phs_fks_config_configure
<<phs fks: procedures>>=
subroutine phs_fks_config_configure (phs_config, sqrts, &
sqrts_fixed, cm_frame, azimuthal_dependence, rebuild, &
ignore_mismatch, nlo_type)
class(phs_fks_config_t), intent(inout) :: phs_config
real(default), intent(in) :: sqrts
logical, intent(in), optional :: sqrts_fixed
logical, intent(in), optional :: cm_frame
logical, intent(in), optional :: azimuthal_dependence
logical, intent(in), optional :: rebuild
logical, intent(in), optional :: ignore_mismatch
type(string_t), intent(inout), optional :: nlo_type
if (present (nlo_type)) then
if (nlo_type /= 'Real') &
call msg_fatal ("FKS config has to be called with nlo_type = 'Real'")
end if
phs_config%n_par = phs_config%n_par + 3
!!! Channel equivalences not accessible yet
phs_config%provides_equivalences = .false.
end subroutine phs_fks_config_configure
@ %def phs_fks_config_configure
@
<<phs fks: fks config: TBP>>=
procedure :: startup_message => phs_fks_config_startup_message
<<phs fks: procedures>>=
subroutine phs_fks_config_startup_message (phs_config, unit)
class(phs_fks_config_t), intent(in) :: phs_config
integer, intent(in), optional :: unit
call phs_config%phs_wood_config_t%startup_message
end subroutine phs_fks_config_startup_message
@ %def phs_fks_config_startup_message
@
<<phs fks: fks config: TBP>>=
procedure, nopass :: allocate_instance => phs_fks_config_allocate_instance
<<phs fks: procedures>>=
subroutine phs_fks_config_allocate_instance (phs)
class(phs_t), intent(inout), pointer :: phs
allocate (phs_fks_t :: phs)
end subroutine phs_fks_config_allocate_instance
@ %def phs_fks_config_allocate_instance
@
<<phs fks: fks config: TBP>>=
procedure :: set_born_config => phs_fks_config_set_born_config
<<phs fks: procedures>>=
subroutine phs_fks_config_set_born_config (phs_config, phs_cfg_born)
class(phs_fks_config_t), intent(inout) :: phs_config
type(phs_wood_config_t), intent(in), target :: phs_cfg_born
phs_config%forest = phs_cfg_born%forest
phs_config%n_channel = phs_cfg_born%n_channel
phs_config%n_par = phs_cfg_born%n_par
phs_config%sqrts = phs_cfg_born%sqrts
phs_config%par = phs_cfg_born%par
phs_config%sqrts_fixed = phs_cfg_born%sqrts_fixed
phs_config%azimuthal_dependence = phs_cfg_born%azimuthal_dependence
phs_config%provides_chains = phs_cfg_born%provides_chains
phs_config%chain = phs_cfg_born%chain
end subroutine phs_fks_config_set_born_config
@ %def phs_fks_config_set_born_config
@ Keep score about the real kinematics.
<<phs fks: public>>=
public :: kinematics_counter_t
<<phs fks: types>>=
type :: kinematics_counter_t
integer :: n_bins = 0
integer, dimension(:), allocatable :: histo_xi
integer, dimension(:), allocatable :: histo_xi_tilde
integer, dimension(:), allocatable :: histo_xi_max
integer, dimension(:), allocatable :: histo_y
integer, dimension(:), allocatable :: histo_phi
contains
<<phs fks: kinematics counter: TBP>>
end type kinematics_counter_t
@ %def kinematics_counter_t
@
<<phs fks: kinematics counter: TBP>>=
procedure :: init => kinematics_counter_init
<<phs fks: procedures>>=
subroutine kinematics_counter_init (counter, n_bins)
class(kinematics_counter_t), intent(inout) :: counter
integer, intent(in) :: n_bins
counter%n_bins = n_bins
allocate (counter%histo_xi (n_bins), counter%histo_xi_tilde (n_bins))
allocate (counter%histo_y (n_bins), counter%histo_phi (n_bins))
allocate (counter%histo_xi_max (n_bins))
counter%histo_xi = 0
counter%histo_xi_tilde = 0
counter%histo_xi_max = 0
counter%histo_y = 0
counter%histo_phi = 0
end subroutine kinematics_counter_init
@ %def kinematics_counter_init
@
<<phs fks: kinematics counter: TBP>>=
procedure :: record => kinematics_counter_record
<<phs fks: procedures>>=
subroutine kinematics_counter_record (counter, xi, xi_tilde, &
xi_max, y, phi)
class(kinematics_counter_t), intent(inout) :: counter
real(default), intent(in), optional :: xi, xi_tilde, xi_max
real(default), intent(in), optional :: y, phi
if (counter%n_bins > 0) then
if (present (xi)) then
call fill_histogram (counter%histo_xi, xi, &
0.0_default, 1.0_default)
end if
if (present (xi_tilde)) then
call fill_histogram (counter%histo_xi_tilde, xi_tilde, &
0.0_default, 1.0_default)
end if
if (present (xi_max)) then
call fill_histogram (counter%histo_xi_max, xi_max, &
0.0_default, 1.0_default)
end if
if (present (y)) then
call fill_histogram (counter%histo_y, y, -1.0_default, 1.0_default)
end if
if (present (phi)) then
call fill_histogram (counter%histo_phi, phi, 0.0_default, twopi)
end if
end if
contains
subroutine fill_histogram (histo, value, val_min, val_max)
integer, dimension(:), allocatable :: histo
real(default), intent(in) :: value, val_min, val_max
real(default) :: step, lo, hi
integer :: bin
step = (val_max-val_min) / counter%n_bins
do bin = 1, counter%n_bins
lo = (bin-1) * step
hi = bin * step
if (value >= lo .and. value < hi) then
histo (bin) = histo (bin) + 1
exit
end if
end do
end subroutine fill_histogram
end subroutine kinematics_counter_record
@ %def kinematics_counter_record
<<phs fks: kinematics counter: TBP>>=
procedure :: display => kinematics_counter_display
<<phs fks: procedures>>=
subroutine kinematics_counter_display (counter)
class(kinematics_counter_t), intent(in) :: counter
print *, 'xi: ', counter%histo_xi
print *, 'xi_tilde: ', counter%histo_xi_tilde
print *, 'xi_max: ', counter%histo_xi_max
print *, 'y: ', counter%histo_y
print *, 'phi: ', counter%histo_phi
end subroutine kinematics_counter_display
@ %def kinematics_counter_display
@ The fks phase space type contains the wood phase space and
separately the in- and outcoming momenta for the real process and the
corresponding Born momenta. Additionally, there are the variables
$\xi$,$\xi_{max}$, $y$ and $\phi$ which are used to create the real
phase space, as well as the jacobian and its corresponding soft and
collinear limit. Lastly, the array \texttt{ch\_to\_em} connects each
channel with an emitter.
<<phs fks: public>>=
public :: phs_fks_t
<<phs fks: types>>=
type, extends (phs_wood_t) :: phs_fks_t
type(vector4_t), dimension(:), allocatable :: p_born
type(vector4_t), dimension(:), allocatable :: q_born
type(vector4_t), dimension(:), allocatable :: p_real
type(vector4_t), dimension(:), allocatable :: q_real
type(vector4_t), dimension(:), allocatable :: p_born_tot
real(default), dimension(3) :: r_real
real(default) :: E_beam_tot
real(default) :: E_gluon
real(default) :: mrec2
real(default), dimension(:), allocatable :: m2
real(default) :: xi_tilde, phi
real(default), dimension(:), allocatable :: xi_max
real(default), dimension(:), allocatable :: y
real(default) :: xi_min, y_max
real(default) :: y_soft
real(default), dimension(:), allocatable :: jac_rand
real(default), dimension(3) :: jac
integer, dimension(:), allocatable :: emitters
type(kinematics_counter_t) :: counter
logical :: massive_phsp = .false.
logical :: perform_generation = .true.
contains
<<phs fks: phs fks: TBP>>
end type phs_fks_t
@ %def phs_fks_t
@
<<phs fks: interfaces>>=
interface compute_beta
module procedure compute_beta_massless
module procedure compute_beta_massive
end interface
interface get_xi_max_fsr
module procedure get_xi_max_fsr_massless
module procedure get_xi_max_fsr_massive
end interface
@ %def interfaces
@ Initializer for the phase space. Calls the initialization of the
corresponding Born phase space, sets up the
channel-emitter-association and allocates space for the momenta.
<<phs fks: phs fks: TBP>>=
procedure :: init => phs_fks_init
<<phs fks: procedures>>=
subroutine phs_fks_init (phs, phs_config)
class(phs_fks_t), intent(out) :: phs
class(phs_config_t), intent(in), target :: phs_config
call phs%base_init (phs_config)
select type (phs_config)
type is (phs_fks_config_t)
phs%config => phs_config
phs%forest = phs_config%forest
end select
deallocate (phs%r)
allocate (phs%r (phs%config%n_par-3, phs%config%n_channel)); phs%r = 0
select type(phs)
type is (phs_fks_t)
call phs%init_momenta (phs_config)
allocate (phs%m2 (phs_config%n_tot))
allocate (phs%xi_max (phs_config%n_tot))
allocate (phs%y (phs_config%n_tot))
allocate (phs%jac_rand (phs_config%n_tot))
phs%xi_max = 0._default
phs%y = 0._default
phs%jac_rand = 1._default
phs%jac = 1._default
phs%xi_min = 1d-5
phs%y_max = 1._default
end select
end subroutine phs_fks_init
@ %def phs_fks_init
@
<<phs fks: phs fks: TBP>>=
procedure :: final => phs_fks_final
<<phs fks: procedures>>=
subroutine phs_fks_final (object)
class(phs_fks_t), intent(inout) :: object
end subroutine phs_fks_final
@ %def phs_fks_final
@
<<phs fks: phs fks: TBP>>=
procedure :: init_momenta => phs_fks_init_momenta
<<phs fks: procedures>>=
subroutine phs_fks_init_momenta (phs, phs_config)
class(phs_fks_t), intent(inout) :: phs
class(phs_config_t), intent(in) :: phs_config
allocate (phs%p_born (phs_config%n_in))
allocate (phs%q_born (phs_config%n_out-1))
allocate (phs%p_real (phs_config%n_in))
allocate (phs%q_real (phs_config%n_out-1))
allocate (phs%p_born_tot (phs%config%n_in + phs%config%n_out-1))
end subroutine phs_fks_init_momenta
@ %def phs_fks_init_momenta
@ Evaluate selected channel. First, the subroutine calls the
evaluation procedure of the underlying Born phase space, using $n_r -
3$ random numbers. Then, the remaining three random numbers are used
to create $\xi$, $y$ and $\phi$, from which the real momenta are
calculated from the Born momenta.
<<phs fks: phs fks: TBP>>=
procedure :: evaluate_selected_channel => phs_fks_evaluate_selected_channel
<<phs fks: procedures>>=
subroutine phs_fks_evaluate_selected_channel (phs, c_in, r_in)
class(phs_fks_t), intent(inout) :: phs
integer, intent(in) :: c_in
real(default), intent(in), dimension(:) :: r_in
integer :: em
real(default), dimension(:), allocatable :: r_born
integer :: n_r_born
integer :: n_in
n_in = phs%config%n_in
if (phs%perform_generation) then
n_r_born = size(r_in) - 3
else
n_r_born = size (r_in)
end if
allocate (r_born (n_r_born))
r_born = r_in(1:n_r_born)
call phs%phs_wood_t%evaluate_selected_channel (c_in, r_born)
phs%p_born = phs%phs_wood_t%p
phs%q_born = phs%phs_wood_t%q
phs%p_born_tot (1:n_in) = phs%p_born
phs%p_born_tot (n_in+1:) = phs%q_born
phs%r(:,c_in) = r_in(1:n_r_born)
if (phs%perform_generation) call phs%generate_real_kinematics &
(r_in(n_r_born+1:n_r_born+3))
end subroutine phs_fks_evaluate_selected_channel
@ %def phs_fks_evaluate_selected_channel
@
<<phs fks: phs fks: TBP>>=
procedure :: evaluate_other_channels => phs_fks_evaluate_other_channels
<<phs fks: procedures>>=
subroutine phs_fks_evaluate_other_channels (phs, c_in)
class(phs_fks_t), intent(inout) :: phs
integer, intent(in) :: c_in
integer :: c
call phs%phs_wood_t%evaluate_other_channels (c_in)
phs%r_defined = .true.
end subroutine phs_fks_evaluate_other_channels
@ %def phs_fks_evaluate_other_channels
@
<<phs fks: phs fks: TBP>>=
procedure :: get_mcpar => phs_fks_get_mcpar
<<phs fks: procedures>>=
subroutine phs_fks_get_mcpar (phs, c, r)
class(phs_fks_t), intent(in) :: phs
integer, intent(in) :: c
real(default), dimension(:), intent(out) :: r
integer :: n_r_born
n_r_born = size (phs%r(:,c))
r(1:n_r_born) = phs%r(:,c)
if (phs%perform_generation) &
r(n_r_born+1:) = phs%r_real
end subroutine phs_fks_get_mcpar
@ %def phs_fks_get_mcpar
@
<<phs fks: phs fks: TBP>>=
procedure :: get_real_kinematics => phs_fks_get_real_kinematics
<<phs fks: procedures>>=
subroutine phs_fks_get_real_kinematics (phs, xit, y, phi, xi_max, jac, jac_rand)
class(phs_fks_t), intent(inout) :: phs
real(default), intent(out), dimension(:), allocatable :: xi_max
real(default), intent(out) :: xit
real(default), intent(out), dimension(:), allocatable :: y
real(default), intent(out) :: phi
real(default), intent(out), dimension(3) :: jac
real(default), intent(out), dimension(:), allocatable :: jac_rand
xit = phs%xi_tilde
y = phs%y
phi = phs%phi
xi_max = phs%xi_max
jac = phs%jac
jac_rand = phs%jac_rand
end subroutine phs_fks_get_real_kinematics
@ %def phs_fks_get_real_kinematics
@
<<phs fks: phs fks: TBP>>=
procedure :: enable_massive_emitters => &
phs_fks_enable_massive_emitters
<<phs fks: procedures>>=
subroutine phs_fks_enable_massive_emitters (phs_fks)
class(phs_fks_t), intent(inout) :: phs_fks
phs_fks%massive_phsp = .true.
end subroutine phs_fks_enable_massive_emitters
@ %def phs_fks_enable_massive_emitters
@
<<phs fks: phs fks: TBP>>=
procedure :: set_emitters => phs_fks_set_emitters
<<phs fks: procedures>>=
subroutine phs_fks_set_emitters (phs, emitters)
class(phs_fks_t), intent(inout) :: phs
integer, intent(in), dimension(:), allocatable :: emitters
phs%emitters = emitters
end subroutine phs_fks_set_emitters
@ %def phs_fks_set_emitters
@
<<phs fks: phs fks: TBP>>=
procedure :: get_born_momenta => phs_fks_get_born_momenta
<<phs fks: procedures>>=
subroutine phs_fks_get_born_momenta (phs, p)
class(phs_fks_t), intent(inout) :: phs
type(vector4_t), intent(out), dimension(:) :: p
p(1:phs%config%n_in) = phs%p_born
p(phs%config%n_in+1:) = phs%q_born
end subroutine phs_fks_get_born_momenta
@ %def phs_fks_get_born_momenta
@
<<phs fks: phs fks: TBP>>=
procedure :: get_outgoing_momenta => phs_fks_get_outgoing_momenta
<<phs fks: procedures>>=
subroutine phs_fks_get_outgoing_momenta (phs, q)
class(phs_fks_t), intent(in) :: phs
type(vector4_t), intent(out), dimension(:) :: q
q = phs%q_real
end subroutine phs_fks_get_outgoing_momenta
@ %def phs_fks_get_outgoing_momenta
@
<<phs fks: phs fks: TBP>>=
procedure :: get_incoming_momenta => phs_fks_get_incoming_momenta
<<phs fks: procedures>>=
subroutine phs_fks_get_incoming_momenta (phs, p)
class(phs_fks_t), intent(in) :: phs
type(vector4_t), intent(inout), dimension(:), allocatable :: p
p = phs%p_real
end subroutine phs_fks_get_incoming_momenta
@ %def phs_fks_get_incoming_momenta
@
<<phs fks: phs fks: TBP>>=
procedure :: display_kinematics => phs_fks_display_kinematics
<<phs fks: procedures>>=
subroutine phs_fks_display_kinematics (phs)
class(phs_fks_t), intent(in) :: phs
! call phs%counter%display ()
end subroutine phs_fks_display_kinematics
@ %def phs_fks_display_kinematics
@
<<phs fks: phs fks: TBP>>=
procedure :: generate_real_kinematics => phs_fks_generate_real_kinematics
<<phs fks: procedures>>=
subroutine phs_fks_generate_real_kinematics (phs, r_in)
class(phs_fks_t), intent(inout) :: phs
real(default), intent(in), dimension(:) :: r_in
integer :: em
real(default) :: beta
if (size (r_in) /= 3) call msg_fatal &
("Real kinematics need to be generated using three random numbers!")
phs%E_beam_tot = phs%p_born(1)%p(0) + phs%p_born(2)%p(0)
!!! Jacobian corresponding to the transformation rand -> (xi, y, phi)
phs%jac_rand = 1.0
!!! Produce real momentum
phs%phi = r_in(3)*twopi
phs%xi_tilde = phs%xi_min + r_in(1)*(1-phs%xi_min)
phs%jac_rand = phs%jac_rand * (1-phs%xi_min)
do em = 1, phs%config%n_tot
if (any (phs%emitters == em)) then
if (phs%massive_phsp) then
phs%m2(em) = phs%p_born_tot(em)**2
beta = beta_emitter (phs%E_beam_tot, phs%p_born_tot(em))
phs%y(em) = 1.0/beta * (1-(1+beta) * &
exp(-r_in(2)*log((1+beta)/(1-beta))))
phs%jac_rand(em) = phs%jac_rand(em) * &
(1-beta*phs%y(em))*log((1+beta)/(1-beta))/beta
phs%xi_max (em) = get_xi_max_fsr &
(phs%p_born_tot, em, phs%m2(em), phs%y(em))
else
phs%m2(em) = 0._default
phs%xi_max (em) = get_xi_max_fsr (phs%p_born_tot, em)
end if
end if
end do
if (.not. phs%massive_phsp) then
phs%y = (1-2*r_in(2))*phs%y_max
phs%jac_rand = phs%jac_rand*3*(1-phs%y(1)**2)
phs%y = 1.5_default*(phs%y - phs%y**3/3)
end if
- phs%volume = phs%volume * phs%config%sqrts**2 / (8*twopi2)
phs%r_real = r_in
contains
function beta_emitter (q0, p) result (beta)
real(default), intent(in) :: q0
type(vector4_t), intent(in) :: p
real(default) :: beta
real(default) :: m2, mrec2, k0_max
m2 = p**2
mrec2 = (q0-p%p(0))**2 - p%p(1)**2 - p%p(2)**2 - p%p(3)**2
k0_max = (q0**2-mrec2+m2)/(2*q0)
beta = sqrt(1-m2/k0_max**2)
end function beta_emitter
end subroutine phs_fks_generate_real_kinematics
@ %def phs_fks_generate_real_kinematics
@
\subsection{Creation of the real phase space - FSR}
At this point, the Born phase space has been generated, as well as the
three random variables $\xi$, $y$ and $\phi$. The question is how the
real phase space is generated for a final-state emission
configuration. We work with two different sets of momenta, the Born
configuration $\Bigl\{ \bar{k}_{\oplus}, \bar{k}_{\ominus}, \bar{k}_{1}, ...,
\bar{k}_{n} \Bigr\}$ and the real configuration $\Bigl\{ k_{\oplus},
k_{\ominus}, k_1,..., k_n, k_{n+1} \Bigr\}$. We define the momentum of
the emitter to be on the $n$-th position and the momentum of the
radiated particle to be at position $n+1$. The magnitude of the
spatial component of k is denoted by $\underline{k}$.
For final-state emissions, it is $\bar{k}_\oplus = k_\oplus$ and
$\bar{k}_\ominus = k_\ominus$. Thus, the center-of-mass systems
coincide and it is
\begin{equation}
q = \sum_{i=1}^n \bar{k}_i = \sum_{i=1}^{n+1} k_i,
\end{equation}
with $\vec{q} = 0$ and $q^2 = \left(q^0\right)^2$.
We want to construct the real phase space from the Born phase space
using three random numbers. They are defined as follows:
\begin{itemize}
\item $\xi = \frac{2k_{n+1}^0}{\sqrt{s}} \in [0, \xi_{max}]$, where
$k_{n+1}$ denotes the four-momentum of the radiated particle.
\item $y = \cos\theta = \frac{\vec{k}_n \cdot
\vec{k}_{n+1}}{\underline{k}_n \underline{k}_{n+1}}$ is the
splitting angle.
\item The angle between tho two splitting particles in the transversal
plane, $phi \in [0,2\pi]$.
\end{itemize}
Further, $k_{rec} = \sum_{i=1}^{n-1} k_i$ denotes the sum of all
recoiling momenta.
%\begin{dubious}
% Note that this calculation only works for massless particles.
%\end{dubious}
<<phs fks: phs fks: TBP>>=
procedure :: generate_fsr => phs_fks_generate_fsr
<<phs fks: procedures>>=
subroutine phs_fks_generate_fsr (phs, emitter, p_born, p_real)
!!! Important: Momenta must be input in the center-of-mass frame
class(phs_fks_t), intent(inout) :: phs
integer, intent(in) :: emitter
type(vector4_t), intent(in), dimension(:) :: p_born
type(vector4_t), intent(out), dimension(:), allocatable :: p_real
integer nlegborn, nlegreal
type(vector4_t) :: q
real(default) :: q0, q2, uk_np1, uk_n
real(default) :: uk_rec, k_rec0
type(vector3_t) :: k_n_born, k
real(default) :: uk_n_born
real(default) :: uk, k2, k0_n
real(default) :: cpsi, beta
type(vector3_t) :: vec, vec_orth
type(lorentz_transformation_t) :: rot, lambda
integer :: i
real(default) :: xi, y, phi
xi = phs%xi_tilde * phs%xi_max(emitter)
y = phs%y(emitter)
phi = phs%phi
nlegreal = phs%config%n_in + phs%config%n_out
nlegborn = nlegreal-1
if (emitter <= 2 .or. emitter > nlegborn) then
call msg_fatal ("Generate FSR phase space: Invalid emitter!")
end if
allocate (p_real (nlegreal))
p_real(1) = p_born(1)
p_real(2) = p_born(2)
q = p_born(1) + p_born(2)
q0 = phs%E_beam_tot
q2 = q**2
phs%E_gluon = q0*xi/2
uk_np1 = phs%E_gluon
k_n_born = p_born(emitter)%p(1:3)
uk_n_born = k_n_born**1
phs%mrec2 = (q-p_born(emitter))**2
if (phs%massive_phsp) then
call phs%compute_emitter_kinematics (emitter, k0_n, uk_n, uk)
else
call phs%compute_emitter_kinematics (emitter, uk_n, uk)
phs%y_soft = y
k0_n = uk_n
end if
vec = uk_n / uk_n_born * k_n_born
vec_orth = create_orthogonal (vec)
p_real(emitter)%p(0) = k0_n
p_real(emitter)%p(1:3) = vec%p(1:3)
cpsi = (uk_n**2 + uk**2 - uk_np1**2) / (2*(uk_n * uk))
!!! This is to catch the case where cpsi = 1, but numerically
!!! turns out to be slightly larger than 1.
if (cpsi > 1._default) cpsi = 1._default
rot = rotation (cpsi, -sqrt (1-cpsi**2), vec_orth)
p_real(emitter) = rot*p_real(emitter)
vec = uk_np1 / uk_n_born * k_n_born
vec_orth = create_orthogonal (vec)
p_real(nlegreal)%p(0) = uk_np1
p_real(nlegreal)%p(1:3) = vec%p(1:3)
cpsi = (uk_np1**2 + uk**2 - uk_n**2) / (2*(uk_np1 * uk))
rot = rotation (cpsi, sqrt (1-cpsi**2), vec_orth)
p_real(nlegreal) = rot*p_real(nlegreal)
@ Construction of the recoiling momenta. The reshuffling of momenta
must not change the invariant mass of the recoiling system, which
means $k_{\rm{rec}}^2 = \bar{k_{\rm{rec}}}^2$. Therefore, the momenta
are related by a boost, $\bar{k}_i = \Lambda k_i$. The boost parameter
is
\begin{equation*}
\beta = \frac{q^2 - (k_{\rm{rec}}^0 +
\underline{k}_{\rm{rec}})^2}{q^2 + (k_{\rm{rec}}^0 +
\underline{k}_{\rm{rec}})^2}
\end{equation*}
<<phs fks: procedures>>=
k_rec0 = q0 - p_real(emitter)%p(0) - p_real(nlegreal)%p(0)
uk_rec = sqrt (k_rec0**2 - phs%mrec2)
if (phs%massive_phsp) then
beta = compute_beta (q2, k_rec0, uk_rec, &
p_born(emitter)%p(0), uk_n_born)
else
beta = compute_beta (q2, k_rec0, uk_rec)
end if
k = p_real(emitter)%p(1:3) + p_real(nlegreal)%p(1:3)
vec%p(1:3) = 1/uk*k%p(1:3)
lambda = boost (beta/sqrt(1-beta**2), vec)
do i = 3, nlegborn
if (i /= emitter) then
p_real(i) = lambda * p_born(i)
end if
end do
vec%p(1:3) = p_born(emitter)%p(1:3)/uk_n_born
rot = rotation (cos(phi), sin(phi), vec)
p_real(nlegreal) = rot * p_real(nlegreal)
p_real(emitter) = rot * p_real(emitter)
@ The factor $\frac{q^2}{(4\pi)^3}$ is not included here since it is
supplied during phase space generation. Also, we already divide by
$\xi$.
<<phs fks: procedures>>=
if (phs%massive_phsp) then
phs%jac(1) = phs%jac(1)*4/q0/uk_n_born/xi
else
k2 = 2*uk_n*uk_np1*(1-y)
phs%jac(1) = uk_n**2/uk_n_born / (uk_n - k2/(2*q0))
end if
!!! Soft jacobian
phs%jac(2) = 1.0
!!! Collinear jacobian
phs%jac(3) = 1-xi/2*q0/uk_n_born
end subroutine phs_fks_generate_fsr
@ %def phs_fks_generate_fsr
<<phs fks: phs fks: TBP>>=
generic :: compute_emitter_kinematics => &
compute_emitter_kinematics_massless, &
compute_emitter_kinematics_massive
procedure :: compute_emitter_kinematics_massless => &
phs_fks_compute_emitter_kinematics_massless
procedure :: compute_emitter_kinematics_massive => &
phs_fks_compute_emitter_kinematics_massive
<<phs fks: procedures>>=
subroutine phs_fks_compute_emitter_kinematics_massless &
(phs, em, uk_em, uk)
class(phs_fks_t), intent(inout) :: phs
integer, intent(in) :: em
real(default), intent(out) :: uk_em, uk
real(default) :: y, k0_np1, q0, q2
y = phs%y(em); k0_np1 = phs%E_gluon
q0 = phs%E_beam_tot; q2 = q0**2
uk_em = (q2 - phs%mrec2 - 2*q0*k0_np1) / (2*(q0 - k0_np1*(1-y)))
uk = sqrt (uk_em**2 + k0_np1**2 + 2*uk_em*k0_np1*y)
end subroutine phs_fks_compute_emitter_kinematics_massless
subroutine phs_fks_compute_emitter_kinematics_massive &
(phs, em, k0_em, uk_em, uk)
class(phs_fks_t), intent(inout) :: phs
integer, intent(in) :: em
real(default), intent(inout) :: k0_em, uk_em, uk
real(default) :: y, k0_np1, q0, q2, mrec2, m2
real(default) :: k0_rec_max, k0_em_max, k0_rec, uk_rec
real(default) :: z, z1, z2
y = phs%y(em); k0_np1 = phs%E_gluon
q0 = phs%E_beam_tot; q2 = q0**2
mrec2 = phs%mrec2; m2 = phs%m2(em)
k0_rec_max = (q2-m2+mrec2)/(2*q0)
k0_em_max = (q2+m2-mrec2)/(2*q0)
z1 = (k0_rec_max+sqrt (k0_rec_max**2-mrec2))/q0
z2 = (k0_rec_max-sqrt (k0_rec_max**2-mrec2))/q0
z = z2 - (z2-z1)*(1+y)/2
k0_em = k0_em_max - k0_np1*z
k0_rec = q0 - k0_np1 - k0_em
uk_em = sqrt(k0_em**2-m2)
uk_rec = sqrt(k0_rec**2 - mrec2)
uk = uk_rec
phs%jac = q0*(z1-z2)/4*k0_np1
phs%y_soft = (2*q2*z-q2-mrec2+m2)/(sqrt(k0_em_max**2-m2)*q0)/2
end subroutine phs_fks_compute_emitter_kinematics_massive
@ %def phs_fks_compute_emitter_kinematics
@
<<phs fks: procedures>>=
function compute_beta_massless (q2, k0_rec, uk_rec) result (beta)
real(default), intent(in) :: q2, k0_rec, uk_rec
real(default) :: beta
beta = (q2 - (k0_rec + uk_rec)**2) / (q2 + (k0_rec + uk_rec)**2)
end function compute_beta_massless
function compute_beta_massive (q2, k0_rec, uk_rec, &
k0_em_born, uk_em_born) result (beta)
real(default), intent(in) :: q2, k0_rec, uk_rec
real(default), intent(in) :: k0_em_born, uk_em_born
real(default) :: beta
real(default) :: k0_rec_born, uk_rec_born, alpha
k0_rec_born = sqrt(q2) - k0_em_born
uk_rec_born = uk_em_born
alpha = (k0_rec+uk_rec)/(k0_rec_born+uk_rec_born)
beta = (1-alpha**2)/(1+alpha**2)
end function compute_beta_massive
@ %def compute_beta
@ The momentum of the radiated particle is computed according to
\begin{equation}
\label{eq:phs fks:compute k_n}
\underline{k}_n = \frac{q^2 - M_{\rm{rec}}^2 -
2q^0\underline{k}_{n+1}}{2(q^0 - \underline{k}_{n+1}(1-y))},
\end{equation}
with $k = k_n + k_{n+1}$ and $M_{\rm{rec}}^2 = k_{\rm{rec}}^2 =
\left(q-k\right)^2$. Because of $\boldsymbol{\bar{k}}_n \parallel
\boldsymbol{k}_n + \boldsymbol{k}_{n+1}$ we find $M_{\rm{rec}}^2 =
\left(q-\bar{k}_n\right)^2$.
Equation \ref{eq:phs fks: compute k_n} follows from the fact that
$\left(\boldsymbol{k} - \boldsymbol{k}_n\right)^2 =
\boldsymbol{k}_{n+1}^2$, which is equivalent to $\boldsymbol{k}_n
\cdot \boldsymbol{k} = \frac{1}{2} \left(\underline{k}_n^2 +
\underline{k}^2 - \underline{k}_{n+1}^2\right)$.\\
$\boldsymbol{k}_n$ and $\boldsymbol{k}_{n+1}$ are obtained by first
setting up vectors parallel to $\boldsymbol{\bar{k}}_n$,
\begin{equation*}
\boldsymbol{k}_n' = \underline{k}_n
\frac{\bar{\pmb{k}}_n}{\underline{\bar{k}}_n}, \quad \pmb{k}_{n+1}'
= \underline{k}_{n+1}\frac{\bar{\pmb{k}}_n}{\underline{\bar{k}}_n},
\end{equation*}
and then rotating these vectors by an amount of $\cos\psi_n =
\frac{\boldsymbol{k}_n\cdot\pmb{k}}{\underline{k}_n \underline{k}}$.
@ The emitted particle cannot have more momentum than the emitter has
in the Born phase space. Thus, there is an upper bound for $\xi$,
determined by the condition $k_{n+1}^0 = \underline{\bar{k}}_n$, which
is equal to
\begin{equation*}
\xi_{\rm{max}} = \frac{2}{\underline{\bar{k}}_n}{q^0}.
\end{equation*}
<<phs fks: procedures>>=
function get_xi_max_fsr_massless (p_born, emitter) result (xi_max)
type(vector4_t), intent(in), dimension(:) :: p_born
integer, intent(in) :: emitter
real(default) :: xi_max
real(default) :: uk_n_born, q0
q0 = vector4_get_component (p_born(1) + p_born(2), 0)
uk_n_born = space_part_norm (p_born(emitter))
xi_max = 2*uk_n_born / q0
end function get_xi_max_fsr_massless
@ %def get_xi_max_fsr_massless
@
<<phs fks: procedures>>=
function get_xi_max_fsr_massive (p_born, emitter, m2, y) result (xi_max)
type(vector4_t), intent(in), dimension(:) :: p_born
integer, intent(in) :: emitter
real(default), intent(in) :: m2, y
real(default) :: xi_max
real(default) :: q0, mrec2, k0_rec_max
real(default) :: z, z1, z2
real(default) :: k_np1_max
q0 = 2*p_born(1)%p(0)
associate (p => p_born(emitter)%p)
mrec2 = (q0-p(0))**2 - p(1)**2 - p(2)**2 - p(3)**2
end associate
k0_rec_max = (q0**2-m2+mrec2)/(2*q0)
z1 = (k0_rec_max+sqrt(k0_rec_max**2-mrec2))/q0
z2 = (k0_rec_max-sqrt(k0_rec_max**2-mrec2))/q0
z = z2 - (z2-z1)*(1+y)/2
k_np1_max = -(q0**2*z**2 - 2*q0*k0_rec_max*z + mrec2)/(2*q0*z*(1-z))
xi_max = 2*k_np1_max/q0
end function get_xi_max_fsr_massive
@ %def get_xi_max_fsr_massive
@
\subsection{Creation of the real phase space - ISR}
<<phs fks: phs fks: TBP>>=
procedure :: generate_isr => phs_fks_generate_isr
<<phs fks: procedures>>=
subroutine phs_fks_generate_isr &
(phs, emitter, p_born, p_real)
!!! Important: Import momenta in the lab frame
class(phs_fks_t), intent(inout) :: phs
integer, intent(in) :: emitter
type(vector4_t), intent(in) , dimension(:), allocatable :: p_born
type(vector4_t), intent(out), dimension(:), allocatable :: p_real
real(default) :: xi, y, phi
integer :: nlegborn, nlegreal
real(default) :: sqrts_beam
real(default) :: sqrts_born
real(default) :: k0_np1
type(vector3_t) :: beta_l, vec_t
real(default) :: beta_t, beta_gamma_l, beta_gamma_t
real(default) :: k_tmp, k_p, k_m, k_p0, k_m0
real(default) :: pt_rad
type(lorentz_transformation_t) :: lambda_t, lambda_l, lambda_l_inv
real(default) :: x_plus, x_minus, barx_plus, barx_minus
integer :: i
xi = phs%xi_tilde * phs%xi_max(emitter)
y = phs%y(emitter)
phi = phs%phi
nlegreal = phs%config%n_in + phs%config%n_out
nlegborn = nlegreal - 1
sqrts_beam = phs%config%sqrts
sqrts_born = sqrt ((p_born(1) + p_born(2))**2)
allocate (p_real (nlegreal))
!!! Create radiation momentum
k0_np1 = sqrts_born*xi/2
!!! There must be the real cm-energy, not the Born one!
!!! s_real = s_born / (1-xi)
!!! Build radiation momentum in the rest frame of the real momenta
p_real(nlegreal)%p(0) = k0_np1
p_real(nlegreal)%p(1) = k0_np1*sqrt(1-y**2)*sin(phi)
p_real(nlegreal)%p(2) = k0_np1*sqrt(1-y**2)*cos(phi)
p_real(nlegreal)%p(3) = k0_np1*y
!!! Boost to lab frame missing
pt_rad = transverse_part (p_real(nlegreal))
beta_t = sqrt (1 + sqrts_born**2 * (1-xi) / pt_rad**2)
beta_gamma_t = 1/sqrt(beta_t)
k_p0 = vector4_get_component (p_born(1), 0)
k_m0 = vector4_get_component (p_born(2), 0)
beta_l%p(1:3) = (p_born(1)%p(1:3) + p_born(2)%p(1:3)) / &
(p_born(1)%p(0) + p_born(2)%p(0))
beta_gamma_l = beta_l**1
beta_l = beta_l / beta_gamma_l
beta_gamma_l = beta_gamma_l / sqrt (1 - beta_gamma_l**2)
vec_t%p(1:2) = p_real(nlegreal)%p(1:2)
vec_t%p(3) = 0._default
call normalize (vec_t)
lambda_l = boost(beta_gamma_l, beta_l)
lambda_t = boost(-beta_gamma_t, vec_t)
lambda_l_inv = boost(-beta_gamma_l, beta_l)
forall (i=3:nlegborn) &
p_real(i) = lambda_l_inv * lambda_t * lambda_l * p_born(i)
!!! Now need access to the x-variables of the IS-partons
barx_plus = 2*vector4_get_component(p_born(1), 0)/sqrts_beam
barx_minus = 2*vector4_get_component(p_born(2), 0)/sqrts_beam
x_plus = barx_plus/sqrt(1-xi) * sqrt ((2-xi*(1-y)) / (2-xi*(1+y)))
x_minus = barx_minus/sqrt(1-xi) * sqrt ((2-xi*(1+y)) / (2-xi*(1-y)))
p_real(1) = x_plus/barx_plus * p_born(1)
p_real(2) = x_minus/barx_minus * p_born(2)
!!! Total nonsense
phs%jac(1) = 1
end subroutine phs_fks_generate_isr
@ %def fks_born_to_real_isr
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Wed, May 14, 10:49 AM (1 d, 1 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111261
Default Alt Text
(202 KB)
Attached To
rWHIZARDSVN whizardsvn
Event Timeline
Log In to Comment