Commit 8a8de746 authored by Juergen Reuter's avatar Juergen Reuter
Browse files

c5344c3 Merge branch '427-real-partition-for-hadron-collision-cross-sections' into 'master'

parent ab7e9b10
......@@ -4,6 +4,11 @@ Use svn log to see detailed changes.
Version 3.0.1+
2021-08-06
Generalize optional partitioning of the NLO real phase space
##################################################################
2021-07-08
RELEASE: version 3.0.1
......
......@@ -862,6 +862,7 @@ REF_OUTPUT_FILES_BASE = \
ext_tests_nlo/ref-output/nlo_eettzj.ref \
ext_tests_nlo/ref-output/nlo_eettzjj.ref \
ext_tests_nlo/ref-output/nlo_eettzz.ref \
ext_tests_nlo/ref-output/nlo_ppjj_real_partition.ref \
ext_tests_nlo/ref-output/nlo_pptttt.ref \
ext_tests_nlo/ref-output/nlo_ppw.ref \
ext_tests_nlo/ref-output/nlo_ppz.ref \
......@@ -1542,6 +1543,7 @@ EXT_NLO_SIN = \
ext_tests_nlo/nlo_eettzj.sin \
ext_tests_nlo/nlo_eettzjj.sin \
ext_tests_nlo/nlo_eettzz.sin \
ext_tests_nlo/nlo_ppjj_real_partition.sin \
ext_tests_nlo/nlo_pptttt.sin \
ext_tests_nlo/nlo_ppw.sin \
ext_tests_nlo/nlo_ppz.sin \
......
include("nlo_settings.sin")
?pacify = true
alias pr = u:U:d:D:s:S:gl
alias jet = u:U:d:D:s:S:gl
$exclude_gauge_splittings = "c:b:t"
?alphas_from_mz = false
?alphas_from_lhapdf = true
alphas_nf = 5
alphas_order = 2
alpha_power = 0
alphas_power = 2
beams = p, p => lhapdf
$lhapdf_file = "MSTW2008nlo68cl"
sqrts = 13 TeV
# Settings for MPI
$rng_method = "rng_stream"
$integration_method = "vamp2"
cuts = let subevt @clustered_jets = cluster [jet] in
let subevt @pt_selected = select if Pt > 80 GeV [@clustered_jets] in
let subevt @eta_selected = select if abs(Eta) < 4 [@pt_selected] in
let subevt @hardest_jets = select if Pt > 100 GeV [@eta_selected] in
count [@eta_selected] >= 2 and count [@hardest_jets] >= 1
scale = let int njet = count [jet] in
if njet == 2 then
( eval Pt [extract index 1 [jet]]
+ eval Pt [extract index 2 [jet]]) / 2
elsif njet == 3 then
( eval Pt [extract index 1 [jet]]
+ eval Pt [extract index 2 [jet]]
+ eval Pt [extract index 3 [jet]]) / 2
else
sqrts
endif
!process ppjj_np = pr, pr => jet, jet { nlo_calculation = real $restrictions="!W+:W-" }
?nlo_use_real_partition = true
real_partition_scale = 10 GeV
process ppjj_wp = pr, pr => jet, jet { nlo_calculation = real $restrictions="!W+:W-" }
!mult_call_real = 2
!integrate (ppjj_np) { iterations = 10:2000000:"gw",5:2000000 }
!integrate (ppjj_wp) { iterations = 10:2000000:"gw",5:2000000 }
!expect (integral(ppjj_np) == integral(ppjj_wp))
!{
! tolerance = sqrt (error(ppjj_np)**2 + error(ppjj_wp)**2)
!}
integrate (ppjj_wp) { iterations = 2:1000:"gw",2:1000 }
?openmp_logging = false
?vis_history = false
?integration_timer = false
openmp_num_threads = 1
?debug_decay = false
?debug_process = false
?debug_verbose = false
?write_raw = false
| Switching to model 'SM', scheme 'GF_MW_MZ'
$blha_ew_scheme = "alpha_qed"
SM.mZ => 9.118800000000E+01
SM.mW => 8.041900200000E+01
SM.mH => 1.250000000000E+02
SM.GF => 1.166390000000E-05
SM.wZ => 0.000000000000E+00
SM.wtop => 0.000000000000E+00
SM.wW => 0.000000000000E+00
SM.wH => 0.000000000000E+00
SM.ms => 0.000000000000E+00
SM.mc => 0.000000000000E+00
SM.mb => 0.000000000000E+00
SM.mtop => 1.732000000000E+02
SM.me => 0.000000000000E+00
SM.mmu => 0.000000000000E+00
SM.mtau => 1.777000000000E+00
SM.alphas => 1.180000000000E-01
?alphas_is_fixed = false
?alphas_from_mz = true
?alphas_from_lambda_qcd = false
alphas_nf = 5
alphas_order = 2
[user variable] jet = PDG(2, -2, 1, -1, 3, -3, 4, -4, 5, -5, 21)
$exclude_gauge_splittings = "t"
$method = "openloops"
seed = 8131
sqrts = 1.000000000000E+03
jet_algorithm = 2
jet_r = 5.000000000000E-01
$integration_method = "vamp2"
$rng_method = "rng_stream"
| End of included 'nlo_settings.sin'
?pacify = true
[user variable] pr = PDG(2, -2, 1, -1, 3, -3, 21)
[user variable] jet = PDG(2, -2, 1, -1, 3, -3, 21)
$exclude_gauge_splittings = "c:b:t"
?alphas_from_mz = false
?alphas_from_lhapdf = true
alphas_nf = 5
alphas_order = 2
alpha_power = 0
alphas_power = 2
$lhapdf_file = "MSTW2008nlo68cl"
sqrts = 1.30000E+04
$rng_method = "rng_stream"
$integration_method = "vamp2"
?nlo_use_real_partition = true
real_partition_scale = 1.00000E+01
$restrictions = "!W+:W-"
| Process library 'nlo_ppjj_real_partition_lib': recorded process 'ppjj_wp'
| Integrate: current process library needs compilation
| Process library 'nlo_ppjj_real_partition_lib': compiling ...
| Process library 'nlo_ppjj_real_partition_lib': writing makefile
| Process library 'nlo_ppjj_real_partition_lib': removing old files
| Process library 'nlo_ppjj_real_partition_lib': writing driver
| Process library 'nlo_ppjj_real_partition_lib': creating source code
| Process library 'nlo_ppjj_real_partition_lib': compiling sources
| Process library 'nlo_ppjj_real_partition_lib': linking
| Process library 'nlo_ppjj_real_partition_lib': loading
| Process library 'nlo_ppjj_real_partition_lib': ... success.
| Integrate: compilation done
| QCD alpha: using a running strong coupling
| RNG: Initializing RNG Stream random-number generator
| RNG: Setting seed for random-number generator to 8131
| Initializing integration for process ppjj_wp:
| Beam structure: p, p => lhapdf
| Beam data (collision):
| p (mass = 0.0000000E+00 GeV)
| p (mass = 0.0000000E+00 GeV)
| sqrts = 1.300000000000E+04 GeV
| Phase space: generating configuration ...
Warning: Intermediate decay of zero-width particle Z may be possible.
Warning: Intermediate decay of zero-width particle W+ may be possible.
Warning: Intermediate decay of zero-width particle W- may be possible.
| Phase space: ... success.
| Phase space: writing configuration file 'ppjj_wp.i1.phs'
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'ppjj_wp.i3.phs'
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'ppjj_wp.i6.phs'
| One-Loop-Provider: Using OpenLoops
| Loading library: [...]
| One-Loop-Provider: Using OpenLoops
| Loading library: [...]
| One-Loop-Provider: Using OpenLoops
| Loading library: [...]
| One-Loop-Provider: Using OpenLoops
| Loading library: [...]
| ------------------------------------------------------------------------
| Process [scattering]: 'ppjj_wp'
| Library name = 'nlo_ppjj_real_partition_lib'
| Process index = 1
| Process components:
| 1: 'ppjj_wp_i1': u:ubar:d:dbar:s:sbar:gl, u:ubar:d:dbar:s:sbar:gl => u:ubar:d:dbar:s:sbar:gl, u:ubar:d:dbar:s:sbar:gl [inactive]
| 2: 'ppjj_wp_i2': dbar:d:ubar:u:sbar:s:gl, dbar:d:ubar:u:sbar:s:gl => d:dbar:u:ubar:s:sbar:gl, d:dbar:u:ubar:s:sbar:gl, d:dbar:u:ubar:s:sbar:gl [openloops], [real]
| 3: 'ppjj_wp_i3': u:ubar:d:dbar:s:sbar:gl, u:ubar:d:dbar:s:sbar:gl => u:ubar:d:dbar:s:sbar:gl, u:ubar:d:dbar:s:sbar:gl [inactive], [virtual]
| 4: 'ppjj_wp_i4': u:ubar:d:dbar:s:sbar:gl, u:ubar:d:dbar:s:sbar:gl => u:ubar:d:dbar:s:sbar:gl, u:ubar:d:dbar:s:sbar:gl [inactive], [subtraction]
| 5: 'ppjj_wp_i5': u:ubar:d:dbar:s:sbar:gl, u:ubar:d:dbar:s:sbar:gl => u:ubar:d:dbar:s:sbar:gl, u:ubar:d:dbar:s:sbar:gl [inactive], [dglap]
| 6: 'ppjj_wp_i6': dbar:d:ubar:u:sbar:s:gl, dbar:d:ubar:u:sbar:s:gl => d:dbar:u:ubar:s:sbar:gl, d:dbar:u:ubar:s:sbar:gl, d:dbar:u:ubar:s:sbar:gl [openloops], [real]
| ------------------------------------------------------------------------
| Phase space: 5 channels, 2 dimensions
| Phase space: found 5 channels, collected in 2 groves.
| Phase space: Using 18 equivalences between channels.
| Phase space: wood
| Phase space: 5 channels, 5 dimensions
| Phase space: found 5 channels, collected in 2 groves.
| Phase space: Using 18 equivalences between channels.
| Phase space: wood
| Phase space: 5 channels, 2 dimensions
| Phase space: found 5 channels, collected in 2 groves.
| Phase space: Using 18 equivalences between channels.
| Phase space: wood
| Phase space: 5 channels, 3 dimensions
| Phase space: found 5 channels, collected in 2 groves.
| Phase space: Using 18 equivalences between channels.
| Phase space: wood
| Phase space: 57 channels, 5 dimensions
| Phase space: found 57 channels, collected in 8 groves.
| Phase space: Using 170 equivalences between channels.
| Phase space: wood
| Beam structure: lhapdf, none => none, lhapdf
| Beam structure: 1 channels, 2 dimensions
| Applying user-defined cuts.
| Using user-defined general scale.
| Starting integration for process 'ppjj_wp' part 'real'
| Integrate: iterations = 2:1000:"gw", 2:1000
| Integrator: 2 chains, 5 channels, 7 dimensions
| Integrator: Using VAMP2 channel equivalences
| Integrator: Write grid header and grids to 'ppjj_wp.m2.vg2'
| Integrator: Grid checkpoint after each iteration
| Integrator: 1000 initial calls, 20 max. bins, stratified = T
| Integrator: VAMP2
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
| VAMP2: Initialize new grids and write to file 'ppjj_wp.m2.vg2'.
| VAMP2: set chain: use chained weights.
1 932 -3.010E+09 1.68E+09 55.73 17.01 0.6
2 980 -5.270E+09 1.61E+09 30.63 9.59 1.3
|-----------------------------------------------------------------------------|
2 1912 -4.183E+09 1.16E+09 27.81 12.16 1.3 0.94 2
|-----------------------------------------------------------------------------|
3 980 -6.109E+09 7.80E+08 12.77 4.00 2.5
4 980 -5.506E+09 5.68E+08 10.32 3.23 3.9
|-----------------------------------------------------------------------------|
4 1960 -5.715E+09 4.59E+08 8.04 3.56 3.9 0.39 2
|=============================================================================|
| Starting integration for process 'ppjj_wp' part 'real'
| Integrate: iterations = 2:1000:"gw", 2:1000
| Integrator: 8 chains, 57 channels, 7 dimensions
| Integrator: Using VAMP2 channel equivalences
| Integrator: Write grid header and grids to 'ppjj_wp.m5.vg2'
| Integrator: Grid checkpoint after each iteration
| Integrator: 1000 initial calls, 20 max. bins, stratified = T
| Integrator: VAMP2
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
| VAMP2: Initialize new grids and write to file 'ppjj_wp.m5.vg2'.
| VAMP2: set chain: use chained weights.
1 1006 9.749E+09 6.29E+09 64.55 20.47 5.7
2 1000 4.877E+09 1.16E+09 23.82 7.53 6.0
|-----------------------------------------------------------------------------|
2 2006 5.037E+09 1.14E+09 22.68 10.16 6.0 0.58 2
|-----------------------------------------------------------------------------|
3 1000 6.069E+09 1.71E+09 28.15 8.90 5.6
4 1000 4.310E+09 9.16E+08 21.26 6.72 6.3
|-----------------------------------------------------------------------------|
4 2000 4.703E+09 8.07E+08 17.17 7.68 6.3 0.82 2
|=============================================================================|
| Integrate: sum of all components
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
1 0 -1.012E+09 9.29E+08 91.78 0.00 -1.3
|=============================================================================|
| There were no errors and 3 warning(s).
| WHIZARD run finished.
|=============================================================================|
......@@ -3645,7 +3645,9 @@ resonances are possible for all ftuples
end subroutine region_data_init_phs_identifiers
@ %def region_data_init_phs_identifiers
@
@ Gathers all ftuples from all ALRs.
There are at most $n \cdot (n-1)$ ftuples with $i$ and $j$ in the final state
and up to $3n$ ftuples with $i$ in the initial state.
<<fks regions: reg data: TBP>>=
procedure :: get_all_ftuples => region_data_get_all_ftuples
<<fks regions: procedures>>=
......@@ -3653,14 +3655,14 @@ resonances are possible for all ftuples
class(region_data_t), intent(in) :: reg_data
type(ftuple_t), intent(inout), dimension(:), allocatable :: ftuples
type(ftuple_t), dimension(:), allocatable :: ftuple_tmp
integer :: i, j, alr
!!! Can have at most n * (n-1) ftuples
integer :: i, j, alr, n_fs
j = 0
allocate (ftuple_tmp (reg_data%n_legs_real * (reg_data%n_legs_real - 1)))
n_fs = reg_data%n_legs_real - reg_data%n_in
allocate (ftuple_tmp (n_fs * (n_fs - 1) + 3 * n_fs))
do i = 1, reg_data%n_regions
associate (region => reg_data%regions(i))
do alr = 1, region%nregions
if (.not. any (region%ftuples(alr) == ftuple_tmp)) then
if (.not. any (ftuple_equal_ireg (region%ftuples(alr), ftuple_tmp))) then
j = j + 1
ftuple_tmp(j) = region%ftuples(alr)
end if
......@@ -8986,7 +8988,11 @@ extracted higher up.
end subroutine real_partition_fixed_order_write
@ %def real_partition_fixed_order_write
@
@ Implements the suppression factor
\begin{equation}
F(\Phi_{n+1}) = 1 - \prod_{(i,j) \in P_{FKS}} \theta \left[ m_i + m_j + h - \sqrt{(p_j + p_j)^2} \right]
\end{equation}
to split the real matrix element into singular and finite part.
<<real subtraction: real partition fixed order: TBP>>=
procedure :: get_f => real_partition_fixed_order_get_f
<<real subtraction: procedures>>=
......@@ -8994,16 +9000,25 @@ extracted higher up.
real(default) :: f
class(real_partition_fixed_order_t), intent(in) :: partition
type(vector4_t), intent(in), dimension(:) :: p
integer :: i
integer :: i, em
f = zero
do i = 1, size (partition%fks_pairs)
PAIRS: do i = 1, size (partition%fks_pairs)
associate (ii => partition%fks_pairs(i)%ireg)
if ((p(ii(1)) + p(ii(2)))**1 < p(ii(1))**1 + p(ii(2))**1 + partition%scale) then
f = one
exit
if (ii(1) == 0) then
IS: do em = 1, 2
if ((p(em) + p(ii(2)))**1 < p(em)**1 + p(ii(2))**1 + partition%scale) then
f = one
exit PAIRS
end if
end do IS
else
if ((p(ii(1)) + p(ii(2)))**1 < p(ii(1))**1 + p(ii(2))**1 + partition%scale) then
f = one
exit PAIRS
end if
end if
end associate
end do
end do PAIRS
end function real_partition_fixed_order_get_f
@ %def real_partition_fixed_order_get_f
......
......@@ -4819,7 +4819,7 @@ correlations that are lost in the individual particle entries.
public :: particle_set_t
<<Particles: types>>=
type :: particle_set_t
! private !!!
! private !!!
integer :: n_beam = 0
integer :: n_in = 0
integer :: n_vir = 0
......@@ -5162,7 +5162,7 @@ virtual.
end select
end do
end subroutine particle_set_transfer
@ %def particle_set_transfer
@ Insert a new particle as an intermediate into a previously empty position.
......@@ -5228,7 +5228,7 @@ children.
call pset%prt(c)%set_parents ([i])
end do
end subroutine particle_set_insert
@ %def particle_set_insert
@ This should be done after completing all insertions: recover color
assignments for the inserted particles, working backwards from
......@@ -5347,9 +5347,9 @@ with invalid color assignments. If the array size is zero, all is fine.
index = pack ([(i, i = 1, n)], mask)
if (present (prt)) prt = pack (particle_set%prt, mask)
end subroutine particle_set_find_prt_invalid_color
@ %def particle_set_find_prt_invalid_color
@
@
<<Particles: particle set: TBP>>=
generic :: get_momenta => get_momenta_all
generic :: get_momenta => get_momenta_indices
......@@ -5442,7 +5442,7 @@ is possible. The reconstructed momentum is not projected on-shell.
p = particle_set%get_momenta (index)
call particle_set%set_momentum (i, sum (p))
end subroutine particle_set_recover_momentum
@ %def particle_set_recover_momentum
@
<<Particles: particle set: TBP>>=
......@@ -5481,7 +5481,7 @@ is possible. The reconstructed momentum is not projected on-shell.
end subroutine particle_set_replace_outgoing_momenta
@ %def particle_set_replace_outgoing_momenta
@
@
<<Particles: particle set: TBP>>=
procedure :: get_outgoing_momenta => particle_set_get_outgoing_momenta
<<Particles: procedures>>=
......@@ -6720,7 +6720,7 @@ To be on the safe side, allocate four times the original size.
@ %def particle_set_to_hepevt_form
@ This procedure aims at reconstructing the momenta of an interaction,
given a particle set. The main task is to find the original hard process, by
following the event history.
following the event history.
In-state: take those particles which are flagged as [[PRT_INCOMING]]
......@@ -6864,7 +6864,7 @@ parent.
p_child(c) = 0
end if
end do
if (n_shift < 0) then
if (n_shift < 0) then
if (present (success)) then
success = .false.
return
......@@ -6916,14 +6916,14 @@ parent.
call state_flv%match (pdg_i, match, map_i)
if (present (success)) then
success = match
end if
end if
if (.not. match) then
if (present (success)) then
return
else
if (present (success)) then
return
else
call err_mismatch
end if
end if
end if
allocate (mask_p (pset%n_tot), source = .false.)
mask_p (p_in) = .true.
mask_p (p_out) = .true.
......@@ -8522,6 +8522,3 @@ remnants first.
end subroutine particles_9
@ %def particles_9
......@@ -12798,7 +12798,9 @@ stored in the kinematics subobject to their source. This is a side effect.
end subroutine term_instance_return_beam_momenta
 
@ %def term_instance_return_beam_momenta
@
@ Applies the real partition by computing the real partition function $F(\Phi)$
and multiplying either $\mathcal{R}_\text{sin} = \mathcal{R} \cdot F$ or
$\mathcal{R}_\text{fin} = \mathcal{R} \cdot (1-F)$.
<<Instances: term instance: TBP>>=
procedure :: apply_real_partition => term_instance_apply_real_partition
<<Instances: procedures>>=
......@@ -12807,7 +12809,7 @@ stored in the kinematics subobject to their source. This is a side effect.
type(process_t), intent(in) :: process
real(default) :: f, sqme
integer :: i_component
integer :: i_amp, n_amps
integer :: i_amp, n_amps, qn_index
logical :: is_subtraction
i_component = term%config%i_component
if (process%component_is_selected (i_component) .and. &
......@@ -12815,14 +12817,18 @@ stored in the kinematics subobject to their source. This is a side effect.
is_subtraction = process%get_component_type (i_component) == COMP_REAL_SING &
.and. term%k_term%emitter < 0
if (is_subtraction) return
select case (process%get_component_type (i_component))
case (COMP_REAL_FIN)
call term%connected%trace%set_duplicate_flv_zero()
end select
select type (pcm => process%get_pcm_ptr ())
type is (pcm_nlo_t)
f = pcm%real_partition%get_f (term%p_hard)
end select
n_amps = term%connected%trace%get_n_matrix_elements ()
do i_amp = 1, n_amps
sqme = real (term%connected%trace%get_matrix_element ( &
term%connected%trace%get_qn_index (i_amp, i_sub = 0)))
qn_index = term%connected%trace%get_qn_index (i_amp, i_sub = 0)
sqme = real (term%connected%trace%get_matrix_element (qn_index))
if (debug_on) call msg_debug2 (D_PROCESS_INTEGRATION, "term_instance_apply_real_partition")
select type (pcm => term%pcm_instance%config)
type is (pcm_nlo_t)
......@@ -12839,7 +12845,7 @@ stored in the kinematics subobject to their source. This is a side effect.
end select
end select
if (debug_on) call msg_debug2 (D_PROCESS_INTEGRATION, "apply_damping: sqme", sqme)
call term%connected%trace%set_matrix_element (i_amp, cmplx (sqme, zero, default))
call term%connected%trace%set_matrix_element (qn_index, cmplx (sqme, zero, default))
end do
end if
end subroutine term_instance_apply_real_partition
......
......@@ -6883,6 +6883,7 @@ with its density matrix in flavor, color, and helicity space.
module state_matrices
 
<<Use kinds>>
use constants, only: zero
use io_units
use format_utils, only: pac_fmt
use format_defs, only: FMT_17, FMT_19
......@@ -7958,6 +7959,54 @@ and one wants to reintroduce the previous order of the matrix elements.
end subroutine state_matrix_reorder_me
 
@ %def state_matrix_order_by_flavors
@ Sets all matrix elements whose flavor structure is a duplicate
of another flavor structure to zero. We need this for the real finite to
ignore duplicate flavor structures while keeping the indices identical to the
singular real component.
When comparing the flavor structures, we take into account permutations of final-
state particles. To do this properly, we keep only the non-hard flavors and the
initial-state flavors, i.e. the first two hard flavors fixed.
<<State matrices: state matrix: TBP>>=
procedure :: set_duplicate_flv_zero => state_matrix_set_duplicate_flv_zero
<<State matrices: procedures>>=
subroutine state_matrix_set_duplicate_flv_zero (state)
class(state_matrix_t), intent(inout), target :: state
type(quantum_numbers_t), dimension(state%depth) :: qn
type(flavor_t) :: flv
class(state_flv_content_t), allocatable :: state_flv
logical, dimension(:), allocatable :: hard_mask, sort_mask, duplicate_mask
integer :: i, j, n_in, n_flvs
n_flvs = state%get_depth ()
n_in = 2
!!! TODO (PS-28-07-21) n_in should not be hard coded to work for decays
!!! This assumes that the positions of the non-hard flavors are the same for all flavor structures.
qn = state%get_quantum_number(1)
allocate (hard_mask(n_flvs))
do i = 1, n_flvs
flv = qn(i)%get_flavor()
hard_mask(i) = flv%is_hard_process ()
end do
allocate (sort_mask(n_flvs))
sort_mask = hard_mask
j = 0
do i = 1, n_flvs
if (j == n_in) exit
if (sort_mask(i)) then
sort_mask(i) = .false.
j = j + 1
end if
end do
allocate (state_flv)
call state_flv%fill (state, sort_mask)
call state_flv%find_duplicates (duplicate_mask)
do i = 1, state%get_n_matrix_elements ()
if (duplicate_mask(i)) then
call state%set_matrix_element_single(i, cmplx(zero, zero, default))
end if
end do
end subroutine state_matrix_set_duplicate_flv_zero
@ %def state_matrix_set_duplicate_flv_zero
@ This subroutine sets up the matrix-element array. The leaf nodes
aquire the index values that point to the appropriate matrix-element
entry.
......@@ -9195,7 +9244,7 @@ from that.
integer, dimension(:), allocatable :: pdg, pdg_subset
integer, dimension(:), allocatable :: idx, map_subset, idx_subset, map
type(quantum_numbers_t), dimension(:), allocatable :: qn
integer :: n, d, c, i, j
integer :: n, d, c, i
call state_tmp%init ()
d = state_full%get_depth ()
allocate (flv (d), qn (d), pdg (d), idx (d), map (d))
......@@ -9305,6 +9354,44 @@ must be equal to the length $d$ of the individual flavor-state entries.
end function pacify_complex
 
@ %def pacify_complex
@ Looks for flavor structures that only differ by a permutation
of the masked flavors.
The result is returned in form of a mask which is [[.true.]] at the
positions of a duplicate flavor structure from the second encounter on.
This routine implements the naive approach: We go through all flavor
structures and compare each one with each preceeding one. This works
but is $\mathcal{O}(n^2)$ in the number of flavor structures. Using
a table to remember which flavor structure has already been encountered,
if would be possible to find the duplicates in $\mathcal{O}(n)$.
<<State matrices: state flv content: TBP>>=
procedure :: find_duplicates => state_flv_content_find_duplicates
<<State matrices: procedures>>=
subroutine state_flv_content_find_duplicates (state_flv, duplicate_mask)
class(state_flv_content_t), intent(in) :: state_flv
logical, dimension(:), allocatable, intent(out) :: duplicate_mask
integer, dimension(:), allocatable :: flvst
integer :: i1, i2, n_flvsts
logical :: found_once
n_flvsts = size (state_flv%pdg, 2)
allocate (duplicate_mask (n_flvsts))
duplicate_mask = .false.
do i1 = 1, n_flvsts
found_once = .false.
flvst = state_flv%pdg(:,i1)
do i2 = 1, i1
if (all(flvst == state_flv%pdg(:,i2))) then
if (found_once) then
duplicate_mask(i1) = .true.
exit
else
found_once = .true.
end if
end if
end do
end do
end subroutine state_flv_content_find_duplicates
@ %def state_flv_content_find_duplicates
@
\subsection{Unit tests}
Test module, followed by the corresponding implementation module.
......@@ -11099,6 +11186,16 @@ After this, the value array has to be updated.
end subroutine interaction_add_state
 
@ %def interaction_add_state
@
<<Interactions: interaction: TBP>>=
procedure :: set_duplicate_flv_zero => interaction_set_duplicate_flv_zero
<<Interactions: procedures>>=
subroutine interaction_set_duplicate_flv_zero (int)
class(interaction_t), intent(inout) :: int
call int%state_matrix%set_duplicate_flv_zero ()
end subroutine interaction_set_duplicate_flv_zero
@ %def interaction_set_duplicate_flv_zero
@ Freeze the quantum state: First collapse the quantum state, i.e.,
remove quantum numbers if any mask has changed, then fix the array of
value pointers.
......@@ -11644,7 +11741,7 @@ tagged as part of the hard process.
call it%advance ()
end do
end subroutine interaction_retag_hard_process