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

f7a7a6b Merge branch...

f7a7a6b Merge branch '367-missing-data-in-events-generated-when-reading-existing-evx-file' into 'master'
parent 26b5dc75
......@@ -4,6 +4,9 @@ Use svn log to see detailed changes.
Version 3.0.1+
2021-08-20
Bug fix: correctly reading in NLO fixed order events
2021-08-06
Generalize optional partitioning of the NLO real phase space
......
......@@ -7,6 +7,7 @@
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
[closed]
* Allocate as checkpoints:
......
......@@ -112,6 +112,7 @@ event_excess* => 0.000000000000E+00
Check MD5 sum = F
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Writing to file = eio_raw_1.evx
* Re-read the event
......@@ -120,6 +121,7 @@ event_excess* => 0.000000000000E+00
Check MD5 sum = F
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Reading from file = eio_raw_1.evx
i_prc = 42
......@@ -285,6 +287,7 @@ event_excess* => 0.000000000000E+00
Check MD5 sum = F
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Writing to file = eio_raw_1.evx
* Re-read both events
......@@ -293,6 +296,7 @@ event_excess* => 0.000000000000E+00
Check MD5 sum = F
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Reading from file = eio_raw_1.evx
i_prc = 5
......
......@@ -118,6 +118,7 @@ event_excess* => 0.000000000000E+00
Alternate weights = 2
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Writing to file = eio_raw_2.evx
* Re-read the event
......@@ -127,6 +128,7 @@ event_excess* => 0.000000000000E+00
Alternate weights = 2
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Reading from file = eio_raw_2.evx
i_prc = 42
......
......@@ -113,6 +113,7 @@ event_excess* => 0.000000000000E+00
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Writing to file = event_streams_2.evx
* Reallocate raw eio stream for reading
......@@ -123,6 +124,7 @@ event_excess* => 0.000000000000E+00
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Reading from file = event_streams_2.evx
* Reread event
......
......@@ -10,6 +10,7 @@
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Writing to file = event_streams_3.evx
* Reallocate raw eio stream for reading
......@@ -20,6 +21,7 @@
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Reading from file = event_streams_3.evx
* Reread event
......@@ -31,6 +33,7 @@
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Writing to file = event_streams_3.evx
========================================================================
......@@ -144,6 +147,7 @@ event_excess* => 0.000000000000E+00
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Reading from file = event_streams_3.evx
* Reread two events and display 2nd event
......
......@@ -10,6 +10,7 @@
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Writing to file = event_streams_4.evx
* Reallocate raw eio stream for reading
......@@ -20,6 +21,7 @@
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Reading from file = event_streams_4.evx
* Reallocate modified raw eio stream for reading (fail)
......@@ -29,6 +31,7 @@
Check MD5 sum = T
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Writing to file = event_streams_4.evx
* Repeat ignoring checksum
......@@ -39,6 +42,7 @@
Check MD5 sum = F
Alpha_s from file = F
Scale from file = F
Events for fNLO = F
Reading from file = event_streams_4.evx
* Test output end: event_streams_4
......@@ -12707,9 +12707,11 @@ choose to calculate all terms at once together with the seed kinematics, use
procedure :: compute_hard_kinematics => &
term_instance_compute_hard_kinematics
<<Instances: procedures>>=
subroutine term_instance_compute_hard_kinematics (term, skip_term, success)
subroutine term_instance_compute_hard_kinematics &
(term, recover, skip_term, success)
class(term_instance_t), intent(inout) :: term
integer, intent(in), optional :: skip_term
logical, intent(in), optional :: recover
logical, intent(out) :: success
type(vector4_t), dimension(:), allocatable :: p
if (allocated (term%core_state)) &
......@@ -12718,6 +12720,9 @@ choose to calculate all terms at once together with the seed kinematics, use
if (term%config%i_term_global == skip_term) return
end if
 
if (present (recover)) then
if (recover) return
end if
if (term%nlo_type == NLO_REAL .and. term%k_term%emitter >= 0) then
call term%k_term%evaluate_radiation (term%p_seed, p, success)
select type (config => term%pcm_instance%config)
......@@ -12759,14 +12764,19 @@ left undefined. We remedy this by calling [[receive_kinematics]] once.
procedure :: recover_seed_kinematics => &
term_instance_recover_seed_kinematics
<<Instances: procedures>>=
subroutine term_instance_recover_seed_kinematics (term)
subroutine term_instance_recover_seed_kinematics (term, p_seed_ref)
class(term_instance_t), intent(inout) :: term
integer :: n_in
type(vector4_t), dimension(:), intent(in), optional :: p_seed_ref
n_in = term%k_term%n_in
call term%k_term%get_incoming_momenta (term%p_seed(1:n_in))
associate (int_eff => term%isolated%int_eff)
call int_eff%set_momenta (term%p_seed(1:n_in), outgoing = .false.)
term%p_seed(n_in + 1 : ) = int_eff%get_momenta (outgoing = .true.)
if (present (p_seed_ref)) then
term%p_seed(n_in + 1 : ) = p_seed_ref
else
term%p_seed(n_in + 1 : ) = int_eff%get_momenta (outgoing = .true.)
end if
end associate
call term%isolated%receive_kinematics ()
end subroutine term_instance_recover_seed_kinematics
......@@ -15040,8 +15050,10 @@ recovered, and we are just computing the rest.
procedure :: compute_seed_kinematics => &
process_instance_compute_seed_kinematics
<<Instances: procedures>>=
subroutine process_instance_compute_seed_kinematics (instance, skip_term)
subroutine process_instance_compute_seed_kinematics &
(instance, recover, skip_term)
class(process_instance_t), intent(inout) :: instance
logical, intent(in), optional :: recover
integer, intent(in), optional :: skip_term
integer :: channel, skip_component, i, j
logical :: success
......@@ -15055,6 +15067,9 @@ recovered, and we are just computing the rest.
else
skip_component = 0
end if
if (present (recover)) then
if (recover) return
end if
if (instance%evaluation_status >= STAT_ACTIVATED) then
success = .true.
do i = 1, instance%process%get_n_components ()
......@@ -15126,7 +15141,7 @@ a channel, reconstruct the MC parameter set.
subroutine process_instance_recover_mcpar (instance, i_term)
class(process_instance_t), intent(inout) :: instance
integer, intent(in) :: i_term
integer :: channel
integer :: channel, i
if (instance%evaluation_status >= STAT_EFF_KINEMATICS) then
channel = instance%selected_channel
if (channel == 0) then
......@@ -15134,6 +15149,16 @@ a channel, reconstruct the MC parameter set.
end if
call instance%term(i_term)%recover_mcpar &
(instance%mci_work(instance%i_mci), channel)
if (instance%term(i_term)%nlo_type == NLO_REAL) then
do i = 1, size (instance%term)
if (i /= i_term .and. instance%term(i)%nlo_type == NLO_REAL) then
if (instance%term(i)%active) then
call instance%term(i)%recover_mcpar &
(instance%mci_work(instance%i_mci), channel)
end if
end if
end do
end if
end if
end subroutine process_instance_recover_mcpar
 
......@@ -15164,21 +15189,24 @@ components, from the seed kinematics.
procedure :: compute_hard_kinematics => &
process_instance_compute_hard_kinematics
<<Instances: procedures>>=
subroutine process_instance_compute_hard_kinematics (instance, skip_term)
subroutine process_instance_compute_hard_kinematics (instance, recover, skip_term)
class(process_instance_t), intent(inout) :: instance
integer, intent(in), optional :: skip_term
logical, intent(in), optional :: recover
integer :: i
logical :: success
success = .true.
if (instance%evaluation_status >= STAT_SEED_KINEMATICS) then
do i = 1, size (instance%term)
if (instance%term(i)%active) then
call instance%term(i)%compute_hard_kinematics (skip_term, success)
call instance%term(i)%compute_hard_kinematics &
(recover, skip_term, success)
if (.not. success) exit
!!! Ren scale is zero when this is commented out! Understand!
if (instance%term(i)%nlo_type == NLO_REAL) &
call instance%term(i)%redo_sf_chain (instance%mci_work(instance%i_mci), &
instance%selected_channel)
call instance%term(i)%redo_sf_chain &
(instance%mci_work(instance%i_mci), &
instance%selected_channel)
end if
end do
if (success) then
......@@ -15200,8 +15228,23 @@ for one specific term.
subroutine process_instance_recover_seed_kinematics (instance, i_term)
class(process_instance_t), intent(inout) :: instance
integer, intent(in) :: i_term
if (instance%evaluation_status >= STAT_EFF_KINEMATICS) &
type(vector4_t), dimension(:), allocatable :: p_seed_ref
integer :: i
if (instance%evaluation_status >= STAT_EFF_KINEMATICS) then
call instance%term(i_term)%recover_seed_kinematics ()
if (instance%term(i_term)%nlo_type == NLO_REAL) then
allocate (p_seed_ref (instance%term(i_term)%isolated%int_eff%get_n_out ()))
p_seed_ref = instance%term(i_term)%isolated%int_eff%get_momenta &
(outgoing = .true.)
do i = 1, size (instance%term)
if (i /= i_term .and. instance%term(i)%nlo_type == NLO_REAL) then
if (instance%term(i)%active) then
call instance%term(i)%recover_seed_kinematics (p_seed_ref)
end if
end if
end do
end if
end if
end subroutine process_instance_recover_seed_kinematics
 
@ %def process_instance_recover_seed_kinematics
......@@ -15403,8 +15446,9 @@ terms, iterating over all active process components.
<<Instances: process instance: TBP>>=
procedure :: evaluate_trace => process_instance_evaluate_trace
<<Instances: procedures>>=
subroutine process_instance_evaluate_trace (instance)
subroutine process_instance_evaluate_trace (instance, recover)
class(process_instance_t), intent(inout) :: instance
logical, intent(in), optional :: recover
class(prc_core_t), pointer :: core => null ()
integer :: i, i_real_fin, i_core
real(default) :: alpha_s, alpha_qed
......@@ -15467,6 +15511,9 @@ terms, iterating over all active process components.
end if
end select
end if
if (present (recover)) then
if (recover) return
end if
select case (term%nlo_type)
case (NLO_REAL)
call term%apply_fks (alpha_s, alpha_qed)
......@@ -15836,17 +15883,20 @@ Before calling this, we should call [[choose_mci]].
logical, intent(in) :: update_sqme
logical, intent(in) :: recover_phs
real(default), intent(in), allocatable, optional :: scale_forced
logical :: skip_phs
logical :: skip_phs, recover
call instance%activate ()
instance%evaluation_status = STAT_EFF_KINEMATICS
call instance%recover_hard_kinematics (i_term)
call instance%recover_seed_kinematics (i_term)
call instance%select_channel (channel)
recover = instance%pcm%config%is_nlo ()
if (recover_phs) then
call instance%recover_mcpar (i_term)
call instance%recover_beam_momenta (i_term)
call instance%compute_seed_kinematics (i_term)
call instance%compute_hard_kinematics (i_term)
call instance%compute_seed_kinematics &
(recover = recover, skip_term = i_term)
call instance%compute_hard_kinematics &
(recover = recover, skip_term = i_term)
call instance%compute_eff_kinematics (i_term)
call instance%compute_other_channels (i_term)
else
......@@ -15855,7 +15905,7 @@ Before calling this, we should call [[choose_mci]].
call instance%evaluate_expressions (scale_forced)
if (update_sqme) then
call instance%reset_core_kinematics ()
call instance%evaluate_trace ()
call instance%evaluate_trace (recover)
end if
end subroutine process_instance_recover
 
......
......@@ -10052,14 +10052,17 @@ object in the [[pcm_instance]].
has already been called for one generation point. If this variable,
[[i_evaluation]], is zero, this means that [[evt_nlo_generate_weighted]] is called
for the first time, so that the generation of an $N$-particle event is required.
In all other cases, emission events are generated.\\
In all other cases, emission events are generated.
During a separate integration of the real component, the first event of
each event group will become the counterevent. In this case, we return
the sum of all subtraction matrix elements. \\
the sum of all subtraction matrix elements.
During a combined integration, the first event will be a combination of all
Born-like events. To get the sum of their matrix elements, we subtract the
sum of all real emissions from the sum of all matrix elements as the real
contribution is the only non-Born contribution. \\
contribution is the only non-Born contribution.
Note that the argument named [[probablity]], the use of the routine
[[generate_weighted]] and the procedure we use to generate NLO events via an
event transformation is an abuse of the interface which should be refactored.
......@@ -11240,9 +11243,12 @@ For testing purpose, we allow for providing them explicitly, as an option.
(event%selected_i_mci, event%selected_i_term)
evt => evt%next
end do
if (debug_on) call msg_debug (D_TRANSFORMS, "Before event transformations")
if (debug_on) call msg_debug (D_TRANSFORMS, "event%weight_prc", event%weight_prc)
if (debug_on) call msg_debug (D_TRANSFORMS, "event%sqme_prc", event%sqme_prc)
if (debug_on) call msg_debug &
(D_TRANSFORMS, "Before event transformations")
if (debug_on) call msg_debug &
(D_TRANSFORMS, "event%weight_prc", event%weight_prc)
if (debug_on) call msg_debug &
(D_TRANSFORMS, "event%sqme_prc", event%sqme_prc)
evt => event%transform_first
do while (associated (evt))
call print_transform_name_if_debug ()
......@@ -11300,10 +11306,14 @@ For testing purpose, we allow for providing them explicitly, as an option.
call event%link_particle_set (evt%particle_set)
end if
end if
if (debug_on) call msg_debug (D_TRANSFORMS, "After event transformations")
if (debug_on) call msg_debug (D_TRANSFORMS, "event%weight_prc", event%weight_prc)
if (debug_on) call msg_debug (D_TRANSFORMS, "event%sqme_prc", event%sqme_prc)
if (debug_on) call msg_debug (D_TRANSFORMS, "evt%particle_set_exists", evt%particle_set_exists)
if (debug_on) call msg_debug &
(D_TRANSFORMS, "After event transformations")
if (debug_on) call msg_debug &
(D_TRANSFORMS, "event%weight_prc", event%weight_prc)
if (debug_on) call msg_debug &
(D_TRANSFORMS, "event%sqme_prc", event%sqme_prc)
if (debug_on) call msg_debug &
(D_TRANSFORMS, "evt%particle_set_exists", evt%particle_set_exists)
end if
contains
subroutine print_transform_name_if_debug ()
......@@ -12690,6 +12700,7 @@ initialization, which should enforce some compatibility mode.
logical :: check = .false.
logical :: use_alphas_from_file = .false.
logical :: use_scale_from_file = .false.
logical :: fixed_order_nlo = .false.
integer :: file_version = CURRENT_FILE_VERSION
contains
<<EIO raw: eio raw: TBP>>
......@@ -12715,6 +12726,8 @@ of the current object status.
object%use_alphas_from_file
write (u, "(3x,A,L1)") "Scale from file = ", &
object%use_scale_from_file
write (u, "(3x,A,L1)") "Events for fNLO = ", &
object%fixed_order_nlo
if (object%reading) then
write (u, "(3x,A,A)") "Reading from file = ", char (object%filename)
else if (object%writing) then
......@@ -12747,10 +12760,10 @@ of the current object status.
procedure :: set_parameters => eio_raw_set_parameters
<<EIO raw: procedures>>=
subroutine eio_raw_set_parameters (eio, check, use_alphas_from_file, &
use_scale_from_file, version_string, extension)
use_scale_from_file, fixed_order_nlo, version_string, extension)
class(eio_raw_t), intent(inout) :: eio
logical, intent(in), optional :: check, use_alphas_from_file, &
use_scale_from_file
use_scale_from_file, fixed_order_nlo
type(string_t), intent(in), optional :: version_string
type(string_t), intent(in), optional :: extension
if (present (check)) eio%check = check
......@@ -12758,6 +12771,8 @@ of the current object status.
use_alphas_from_file
if (present (use_scale_from_file)) eio%use_scale_from_file = &
use_scale_from_file
if (present (fixed_order_nlo)) eio%fixed_order_nlo = &
fixed_order_nlo
if (present (version_string)) then
select case (char (version_string))
case ("", "2.2.4")
......@@ -13089,6 +13104,12 @@ a different version of [[event%set_hard_particle_set]].
call event%set_hard_particle_set (pset)
if (eio%use_alphas_from_file .or. eio%use_scale_from_file) then
call event%recalculate (update_sqme = .true.)
if (eio%fixed_order_nlo) then
if (event%weight_prc /= event%weight_ref .and. &
event%weight_prc == 0) then
event%weight_prc = event%weight_ref
end if
end if
end if
call pset%final ()
deallocate (pset)
......@@ -13849,6 +13870,7 @@ which is defined here.
class(event_callback_t), allocatable, intent(in) :: event_callback
logical :: check, keep_beams, keep_remnants, recover_beams
logical :: use_alphas_from_file, use_scale_from_file
logical :: fixed_order_nlo_events
logical :: write_sqme_prc, write_sqme_ref, write_sqme_alt
logical :: output_cross_section, ensure_order
type(string_t) :: lhef_version, lhef_extension, raw_version
......@@ -13876,6 +13898,8 @@ which is defined here.
var_list%get_lval (var_str ("?use_alphas_from_file"))
use_scale_from_file = &
var_list%get_lval (var_str ("?use_scale_from_file"))
fixed_order_nlo_events = &
var_list%get_lval (var_str ("?fixed_order_nlo_events"))
select case (char (method))
case ("raw")
allocate (eio_raw_t :: eio)
......@@ -13888,7 +13912,8 @@ which is defined here.
extension_raw = &
var_list%get_sval (var_str ("$extension_raw"))
call eio%set_parameters (check, use_alphas_from_file, &
use_scale_from_file, raw_version, extension_raw)
use_scale_from_file, fixed_order_nlo_events, &
raw_version, extension_raw)
end select
case ("checkpoint")
allocate (eio_checkpoints_t :: eio)
......
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