Commit 87587001 authored by Juergen Reuter's avatar Juergen Reuter
Browse files

6cf883e Merge branch...

6cf883e Merge branch '435-photon-recombination-needs-to-include-photons-in-recombined-subevent' into 'master'
parent eb1dc8da
......@@ -4,6 +4,9 @@ Use svn log to see detailed changes.
Version 3.0.1+
2021-09-30
Bug fix: keep non-recombined photons in the event record
2021-09-13
Modular NLO event generation with real partition
......
......@@ -5883,6 +5883,54 @@ test.
end subroutine cluster_p
 
@ %def cluster_p
@ Recombine photons with other particles (usually charged leptons and
maybe quarks) given in the same subevent. If [[en0]] is present,
create a mask which is true only for those particles that pass the test.
<<Eval trees: procedures>>=
subroutine photon_recombination_p (subevt, en1, en0)
type(subevt_t), intent(inout) :: subevt
type(eval_node_t), intent(in) :: en1
type(eval_node_t), intent(inout), optional :: en0
logical, dimension(:), allocatable :: mask1
type(prt_t), dimension(:), allocatable :: prt
integer :: n, i
real(default) :: reco_r0
logical :: keep_flv
reco_r0 = en1%photon_rec_r0
n = subevt_get_length (en1%pval)
allocate (prt (n))
do i = 1, n
prt(i) = subevt_get_prt (en1%pval, i)
if (.not. prt_is_recombinable (prt (i))) then
call msg_fatal ("Only charged leptons, quarks, and " //&
"photons can be included in photon recombination.")
end if
end do
if (count (prt_is_photon (prt)) > 1) &
call msg_fatal ("Photon recombination is supported " // &
"only for single photons.")
allocate (mask1 (n))
if (present (en0)) then
do i = 1, n
en0%index = i
en0%prt1 = subevt_get_prt (en1%pval, i)
call eval_node_evaluate (en0)
mask1(i) = en0%lval
end do
else
mask1 = .true.
end if
if (associated (en1%var_list)) then
keep_flv = en1%var_list%get_lval &
(var_str("?keep_flavors_when_recombining"))
else
keep_flv = .false.
end if
call subevt_recombine &
(subevt, en1%pval, mask1, reco_r0, keep_flv)
end subroutine photon_recombination_p
@ %def photon_recombination_p
@ Select all particles of the first argument. If [[en0]] is present,
create a mask which is true only for those particles that pass the
test.
......@@ -6383,64 +6431,6 @@ ignored.
end function no_pp
 
@ %def all_pp any_pp no_pp
@ Recombine photons with the particles of the first argument. If [[en0]] is present,
create a mask which is true only for those particles that pass the test.
<<Eval trees: procedures>>=
subroutine photon_recombination_pp (subevt, en1, en2, en0)
type(subevt_t), intent(inout) :: subevt
type(eval_node_t), intent(in) :: en1
type(eval_node_t), intent(in) :: en2
type(eval_node_t), intent(inout), optional :: en0
logical, dimension(:), allocatable :: mask1
type(prt_t), dimension(:), allocatable :: prt_gam0
type(prt_t) :: prt
integer :: n1, i
real(default) :: reco_r0
logical :: keep_flv
reco_r0 = en1%photon_rec_r0
n1 = subevt_get_length (en1%pval)
if (n1 == 0) then
subevt = en2%pval
else
if (n1 > 1) then
call msg_fatal ("Photon recombination is supported " // &
"only for single photons.")
end if
allocate (prt_gam0 (n1), mask1 (n1))
do i = 1, n1
prt_gam0(i) = subevt_get_prt (en1%pval, i)
end do
if (present (en0)) then
do i = 1, n1
en0%index = i
en0%prt1 = prt_gam0(i)
if (.not. prt_is_photon (en0%prt1)) &
call msg_fatal ("Photon recombination can only " // &
"be applied to photons.")
call eval_node_evaluate (en0)
mask1(i) = en0%lval
end do
else
mask1 = .true.
end if
do i = 1, subevt_get_length (en2%pval)
if (.not. prt_is_recombinable (subevt_get_prt (en2%pval, i))) then
call msg_fatal ("Only charged leptons and QCD partons " //&
"can be recombined with photons")
end if
end do
if (associated (en1%var_list)) then
keep_flv = en1%var_list%get_lval &
(var_str("?keep_flavors_when_recombining"))
else
keep_flv = .false.
end if
call subevt_recombine &
(subevt, en2%pval, prt_gam0(1), mask1(1), reco_r0, keep_flv)
end if
end subroutine photon_recombination_pp
@ %def photon_recombination_pp
The conditional restriction encoded in the [[eval_node_t]] [[en_0]] is
applied only to the photons from [[en1]], not to the objects being
isolated from in [[en2]].
......@@ -9277,6 +9267,11 @@ appropriate (unary/binary) observables.
call var_list%get_rptr (var_str ("jet_p"), en1%jet_p)
call var_list%get_rptr (var_str ("jet_ycut"), en1%jet_ycut)
call var_list%get_rptr (var_str ("jet_dcut"), en1%jet_dcut)
case ("photon_recombination")
en1%var_list => var_list
call eval_node_init_prt_fun_unary &
(en, en1, key, photon_recombination_p)
call var_list%get_rptr (var_str ("photon_rec_r0"), en1%photon_rec_r0)
case ("select")
call eval_node_init_prt_fun_unary (en, en1, key, select_p)
case ("extract")
......@@ -9308,11 +9303,6 @@ appropriate (unary/binary) observables.
call eval_node_init_prt_fun_binary (en, en1, en2, key, select_pp)
case ("sort")
call eval_node_init_prt_fun_binary (en, en1, en2, key, sort_pp)
case ("photon_recombination")
en1%var_list => var_list
call eval_node_init_prt_fun_binary &
(en, en1, en2, key, photon_recombination_pp)
call var_list%get_rptr (var_str ("photon_rec_r0"), en1%photon_rec_r0)
case default
call msg_bug (" Binary particle function '" // char (key) // &
"' undefined")
......
......@@ -3167,7 +3167,8 @@ is not set.
type(prt_t), intent(in) :: prt
flag = prt_is_parton (prt) .or. &
abs(prt%pdg) == TOP_Q .or. &
prt_is_lepton (prt)
prt_is_lepton (prt) .or. &
prt_is_photon (prt)
end function prt_is_recombinable
@ %def prt_is_recombinable
......@@ -4299,46 +4300,83 @@ that have been constructed by FastJet and its interface.
end subroutine subevt_cluster
@ %def subevt_cluster
@ Do recombination.
@ Do recombination. The incoming subevent [[pl]] is left unchanged if
it either does not contain photons at all, or consists just of a
single photon and nothing else or the photon does have a larger $R>R_0$
distance to the nearest other particle or does not fulfill the
[[mask1]] condition. Otherwise, the subevent is one entry shorter and
contains a single recombined particle whose original flavor is kept
depending on the setting [[keep_flv]]. When this subroutine is called,
it is explicitly assumed that there is only one photon. For the
moment, we take here the first photon from the subevent to possibly
recombine and leave this open for generalization.
<<Subevents: public>>=
public :: subevt_recombine
<<Subevents: procedures>>=
subroutine subevt_recombine (subevt, pl, prt, mask1, reco_r0, keep_flv)
subroutine subevt_recombine (subevt, pl, mask1, reco_r0, keep_flv)
type(subevt_t), intent(inout) :: subevt
type(subevt_t), intent(in) :: pl
type(prt_t), intent(in) :: prt
logical, intent(in) :: mask1, keep_flv
type(prt_t), dimension(:), allocatable :: prt_rec
logical, dimension(:), intent(in) :: mask1
logical, intent(in) :: keep_flv
real(default), intent(in) :: reco_r0
real(default), dimension(:), allocatable :: del_rij
integer, dimension(:), allocatable :: i_sortr
type(prt_t) :: prt_comb
logical :: ok
integer :: i, n, pdg_orig
type(prt_t) :: prt_gam, prt_comb
logical :: recombine, ok
integer :: i, n, i_gam, n_gam, n_rec, pdg_orig
allocate (prt_rec (0))
n = subevt_get_length (pl)
allocate (del_rij (n), i_sortr (n))
do i = 1, n
del_rij(i) = eta_phi_distance(prt_get_momentum (prt), &
prt_get_momentum (pl%prt(i)))
end do
i_sortr = order (del_rij)
call subevt_reset (subevt, pl%n_active)
do i = 1, n
if (i == i_sortr(1)) then
if (del_rij (i_sortr (1)) <= reco_r0 .and. mask1) then
pdg_orig = prt_get_pdg (pl%prt (i_sortr (1)))
call prt_combine (prt_comb, prt, pl%prt(i_sortr (1)), ok)
if (ok) then
subevt%prt(i_sortr (1)) = prt_comb
if (keep_flv) call prt_set_pdg &
(subevt%prt(i_sortr (1)), pdg_orig)
n_gam = 0
FIND_FIRST_PHOTON: do i = 1, n
if (prt_is_photon (pl%prt (i))) then
n_gam = n_gam + 1
prt_gam = pl%prt (i)
i_gam = i
exit FIND_FIRST_PHOTON
end if
end do FIND_FIRST_PHOTON
n_rec = n - n_gam
if (n_gam == 0) then
subevt = pl
else
if (n_rec > 0) then
allocate (prt_rec (n_rec))
do i = 1, n_rec
if (i == i_gam) cycle
if (i < i_gam) then
prt_rec(i) = pl%prt(i)
else
prt_rec(i) = pl%prt(i+n_gam)
end if
end do
allocate (del_rij (n_rec), i_sortr (n_rec))
del_rij(1:n_rec) = eta_phi_distance(prt_get_momentum (prt_gam), &
prt_get_momentum (prt_rec(1:n_rec)))
i_sortr = order (del_rij)
recombine = del_rij (i_sortr (1)) <= reco_r0 .and. mask1(i_gam)
if (recombine) then
call subevt_reset (subevt, pl%n_active-n_gam)
do i = 1, n_rec
if (i == i_sortr(1)) then
pdg_orig = prt_get_pdg (prt_rec(i_sortr (1)))
call prt_combine (prt_comb, prt_gam, prt_rec(i_sortr (1)), ok)
if (ok) then
subevt%prt(i_sortr (1)) = prt_comb
if (keep_flv) call prt_set_pdg &
(subevt%prt(i_sortr (1)), pdg_orig)
end if
else
subevt%prt(i) = prt_rec(i)
end if
end do
else
subevt%prt(i) = pl%prt(i)
subevt = pl
end if
else
subevt%prt(i) = pl%prt(i)
subevt = pl
end if
end do
end if
end subroutine subevt_recombine
@ %def subevt_recombine
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment