Index: trunk/src/types/types.nw =================================================================== --- trunk/src/types/types.nw (revision 8747) +++ trunk/src/types/types.nw (revision 8748) @@ -1,8034 +1,8033 @@ %% -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*- % WHIZARD code as NOWEB source: common types and objects %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Sindarin Built-In Types} \includemodulegraph{types} Here, we define a couple of types and objects which are useful both internally for \whizard, and visible to the user, so they correspond to Sindarin types. \begin{description} \item[particle\_specifiers] Expressions for particles and particle alternatives, involving particle names. \item[pdg\_arrays] Integer (PDG) codes for particles. Useful for particle aliases (e.g., 'quark' for $u,d,s$ etc.). \item[jets] Define (pseudo)jets as objects. Functional only if the [[fastjet]] library is linked. (This may change in the future.) \item[subevents] Particle collections built from event records, for use in analysis and other Sindarin expressions \item[analysis] Observables, histograms, and plots. \end{description} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Particle Specifiers} In this module we introduce a type for specifying a particle or particle alternative. In addition to the particle specifiers (strings separated by colons), the type contains an optional flag [[polarized]] and a string [[decay]]. If the [[polarized]] flag is set, particle polarization information should be kept when generating events for this process. If the [[decay]] string is set, it is the ID of a decay process which should be applied to this particle when generating events. In input/output form, the [[polarized]] flag is indicated by an asterisk [[(*)]] in brackets, and the [[decay]] is indicated by its ID in brackets. The [[read]] and [[write]] procedures in this module are not type-bound but generic procedures which handle scalar and array arguments. <<[[particle_specifiers.f90]]>>= <> module particle_specifiers <> use io_units use diagnostics <> <> <> <> contains <> end module particle_specifiers @ %def particle_specifiers @ \subsection{Base type} This is an abstract type which can hold a single particle or an expression. <>= type, abstract :: prt_spec_expr_t contains <> end type prt_spec_expr_t @ %def prt_expr_t @ Output, as a string. <>= procedure (prt_spec_expr_to_string), deferred :: to_string <>= abstract interface function prt_spec_expr_to_string (object) result (string) import class(prt_spec_expr_t), intent(in) :: object type(string_t) :: string end function prt_spec_expr_to_string end interface @ %def prt_spec_expr_to_string @ Call an [[expand]] method for all enclosed subexpressions (before handling the current expression). <>= procedure (prt_spec_expr_expand_sub), deferred :: expand_sub <>= abstract interface subroutine prt_spec_expr_expand_sub (object) import class(prt_spec_expr_t), intent(inout) :: object end subroutine prt_spec_expr_expand_sub end interface @ %def prt_spec_expr_expand_sub @ \subsection{Wrapper type} This wrapper can hold a particle expression of any kind. We need it so we can make variadic arrays. <>= public :: prt_expr_t <>= type :: prt_expr_t class(prt_spec_expr_t), allocatable :: x contains <> end type prt_expr_t @ %def prt_expr_t @ Output as a string: delegate. <>= procedure :: to_string => prt_expr_to_string <>= recursive function prt_expr_to_string (object) result (string) class(prt_expr_t), intent(in) :: object type(string_t) :: string if (allocated (object%x)) then string = object%x%to_string () else string = "" end if end function prt_expr_to_string @ %def prt_expr_to_string @ Allocate the expression as a particle specifier and copy the value. <>= procedure :: init_spec => prt_expr_init_spec <>= subroutine prt_expr_init_spec (object, spec) class(prt_expr_t), intent(out) :: object type(prt_spec_t), intent(in) :: spec allocate (prt_spec_t :: object%x) select type (x => object%x) type is (prt_spec_t) x = spec end select end subroutine prt_expr_init_spec @ %def prt_expr_init_spec @ Allocate as a list/sum and allocate for a given length <>= procedure :: init_list => prt_expr_init_list procedure :: init_sum => prt_expr_init_sum <>= subroutine prt_expr_init_list (object, n) class(prt_expr_t), intent(out) :: object integer, intent(in) :: n allocate (prt_spec_list_t :: object%x) select type (x => object%x) type is (prt_spec_list_t) allocate (x%expr (n)) end select end subroutine prt_expr_init_list subroutine prt_expr_init_sum (object, n) class(prt_expr_t), intent(out) :: object integer, intent(in) :: n allocate (prt_spec_sum_t :: object%x) select type (x => object%x) type is (prt_spec_sum_t) allocate (x%expr (n)) end select end subroutine prt_expr_init_sum @ %def prt_expr_init_list @ %def prt_expr_init_sum @ Return the number of terms. This is unity, except if the expression is a sum. <>= procedure :: get_n_terms => prt_expr_get_n_terms <>= function prt_expr_get_n_terms (object) result (n) class(prt_expr_t), intent(in) :: object integer :: n if (allocated (object%x)) then select type (x => object%x) type is (prt_spec_sum_t) n = size (x%expr) class default n = 1 end select else n = 0 end if end function prt_expr_get_n_terms @ %def prt_expr_get_n_terms @ Transform one of the terms, as returned by the previous method, to an array of particle specifiers. The array has more than one entry if the selected term is a list. This makes sense only if the expression has been completely expanded, so the list contains only atoms. <>= procedure :: term_to_array => prt_expr_term_to_array <>= recursive subroutine prt_expr_term_to_array (object, array, i) class(prt_expr_t), intent(in) :: object type(prt_spec_t), dimension(:), intent(inout), allocatable :: array integer, intent(in) :: i integer :: j if (allocated (array)) deallocate (array) select type (x => object%x) type is (prt_spec_t) allocate (array (1)) array(1) = x type is (prt_spec_list_t) allocate (array (size (x%expr))) do j = 1, size (array) select type (y => x%expr(j)%x) type is (prt_spec_t) array(j) = y end select end do type is (prt_spec_sum_t) call x%expr(i)%term_to_array (array, 1) end select end subroutine prt_expr_term_to_array @ %def prt_expr_term_to_array @ \subsection{The atomic type} The trivial case is a single particle, including optional decay and polarization attributes. \subsubsection{Definition} The particle is unstable if the [[decay]] array is allocated. The [[polarized]] flag and decays may not be set simultaneously. <>= public :: prt_spec_t <>= type, extends (prt_spec_expr_t) :: prt_spec_t private type(string_t) :: name logical :: polarized = .false. type(string_t), dimension(:), allocatable :: decay contains <> end type prt_spec_t @ %def prt_spec_t @ \subsubsection{I/O} Output. Old-style subroutines. <>= public :: prt_spec_write <>= interface prt_spec_write module procedure prt_spec_write1 module procedure prt_spec_write2 end interface prt_spec_write <>= subroutine prt_spec_write1 (object, unit, advance) type(prt_spec_t), intent(in) :: object integer, intent(in), optional :: unit character(len=*), intent(in), optional :: advance character(3) :: adv integer :: u u = given_output_unit (unit) adv = "yes"; if (present (advance)) adv = advance write (u, "(A)", advance = adv) char (object%to_string ()) end subroutine prt_spec_write1 @ %def prt_spec_write1 @ Write an array as a list of particle specifiers. <>= subroutine prt_spec_write2 (prt_spec, unit, advance) type(prt_spec_t), dimension(:), intent(in) :: prt_spec integer, intent(in), optional :: unit character(len=*), intent(in), optional :: advance character(3) :: adv integer :: u, i u = given_output_unit (unit) adv = "yes"; if (present (advance)) adv = advance do i = 1, size (prt_spec) if (i > 1) write (u, "(A)", advance="no") ", " call prt_spec_write (prt_spec(i), u, advance="no") end do write (u, "(A)", advance = adv) end subroutine prt_spec_write2 @ %def prt_spec_write2 @ Read. Input may be string or array of strings. <>= public :: prt_spec_read <>= interface prt_spec_read module procedure prt_spec_read1 module procedure prt_spec_read2 end interface prt_spec_read @ Read a single particle specifier <>= pure subroutine prt_spec_read1 (prt_spec, string) type(prt_spec_t), intent(out) :: prt_spec type(string_t), intent(in) :: string type(string_t) :: arg, buffer integer :: b1, b2, c, n, i b1 = scan (string, "(") b2 = scan (string, ")") if (b1 == 0) then prt_spec%name = trim (adjustl (string)) else prt_spec%name = trim (adjustl (extract (string, 1, b1-1))) arg = trim (adjustl (extract (string, b1+1, b2-1))) if (arg == "*") then prt_spec%polarized = .true. else n = 0 buffer = arg do if (verify (buffer, " ") == 0) exit n = n + 1 c = scan (buffer, "+") if (c == 0) exit buffer = extract (buffer, c+1) end do allocate (prt_spec%decay (n)) buffer = arg do i = 1, n c = scan (buffer, "+") if (c == 0) c = len (buffer) + 1 prt_spec%decay(i) = trim (adjustl (extract (buffer, 1, c-1))) buffer = extract (buffer, c+1) end do end if end if end subroutine prt_spec_read1 @ %def prt_spec_read1 @ Read a particle specifier array, given as a single string. The array is allocated to the correct size. <>= pure subroutine prt_spec_read2 (prt_spec, string) type(prt_spec_t), dimension(:), intent(out), allocatable :: prt_spec type(string_t), intent(in) :: string type(string_t) :: buffer integer :: c, i, n n = 0 buffer = string do n = n + 1 c = scan (buffer, ",") if (c == 0) exit buffer = extract (buffer, c+1) end do allocate (prt_spec (n)) buffer = string do i = 1, size (prt_spec) c = scan (buffer, ",") if (c == 0) c = len (buffer) + 1 call prt_spec_read (prt_spec(i), & trim (adjustl (extract (buffer, 1, c-1)))) buffer = extract (buffer, c+1) end do end subroutine prt_spec_read2 @ %def prt_spec_read2 @ \subsubsection{Constructor} Initialize a particle specifier. <>= public :: new_prt_spec <>= interface new_prt_spec module procedure new_prt_spec module procedure new_prt_spec_polarized module procedure new_prt_spec_unstable end interface new_prt_spec <>= elemental function new_prt_spec (name) result (prt_spec) type(string_t), intent(in) :: name type(prt_spec_t) :: prt_spec prt_spec%name = name end function new_prt_spec elemental function new_prt_spec_polarized (name, polarized) result (prt_spec) type(string_t), intent(in) :: name logical, intent(in) :: polarized type(prt_spec_t) :: prt_spec prt_spec%name = name prt_spec%polarized = polarized end function new_prt_spec_polarized pure function new_prt_spec_unstable (name, decay) result (prt_spec) type(string_t), intent(in) :: name type(string_t), dimension(:), intent(in) :: decay type(prt_spec_t) :: prt_spec prt_spec%name = name allocate (prt_spec%decay (size (decay))) prt_spec%decay = decay end function new_prt_spec_unstable @ %def new_prt_spec @ \subsubsection{Access Methods} Return the particle name without qualifiers <>= procedure :: get_name => prt_spec_get_name <>= elemental function prt_spec_get_name (prt_spec) result (name) class(prt_spec_t), intent(in) :: prt_spec type(string_t) :: name name = prt_spec%name end function prt_spec_get_name @ %def prt_spec_get_name @ Return the name with qualifiers <>= procedure :: to_string => prt_spec_to_string <>= function prt_spec_to_string (object) result (string) class(prt_spec_t), intent(in) :: object type(string_t) :: string integer :: i string = object%name if (allocated (object%decay)) then string = string // "(" do i = 1, size (object%decay) if (i > 1) string = string // " + " string = string // object%decay(i) end do string = string // ")" else if (object%polarized) then string = string // "(*)" end if end function prt_spec_to_string @ %def prt_spec_to_string @ Return the polarization flag <>= procedure :: is_polarized => prt_spec_is_polarized <>= elemental function prt_spec_is_polarized (prt_spec) result (flag) class(prt_spec_t), intent(in) :: prt_spec logical :: flag flag = prt_spec%polarized end function prt_spec_is_polarized @ %def prt_spec_is_polarized @ The particle is unstable if there is a decay array. <>= procedure :: is_unstable => prt_spec_is_unstable <>= elemental function prt_spec_is_unstable (prt_spec) result (flag) class(prt_spec_t), intent(in) :: prt_spec logical :: flag flag = allocated (prt_spec%decay) end function prt_spec_is_unstable @ %def prt_spec_is_unstable @ Return the number of decay channels <>= procedure :: get_n_decays => prt_spec_get_n_decays <>= elemental function prt_spec_get_n_decays (prt_spec) result (n) class(prt_spec_t), intent(in) :: prt_spec integer :: n if (allocated (prt_spec%decay)) then n = size (prt_spec%decay) else n = 0 end if end function prt_spec_get_n_decays @ %def prt_spec_get_n_decays @ Return the decay channels <>= procedure :: get_decays => prt_spec_get_decays <>= subroutine prt_spec_get_decays (prt_spec, decay) class(prt_spec_t), intent(in) :: prt_spec type(string_t), dimension(:), allocatable, intent(out) :: decay if (allocated (prt_spec%decay)) then allocate (decay (size (prt_spec%decay))) decay = prt_spec%decay else allocate (decay (0)) end if end subroutine prt_spec_get_decays @ %def prt_spec_get_decays @ \subsubsection{Miscellaneous} There is nothing to expand here: <>= procedure :: expand_sub => prt_spec_expand_sub <>= subroutine prt_spec_expand_sub (object) class(prt_spec_t), intent(inout) :: object end subroutine prt_spec_expand_sub @ %def prt_spec_expand_sub @ \subsection{List} A list of particle specifiers, indicating, e.g., the final state of a process. <>= public :: prt_spec_list_t <>= type, extends (prt_spec_expr_t) :: prt_spec_list_t type(prt_expr_t), dimension(:), allocatable :: expr contains <> end type prt_spec_list_t @ %def prt_spec_list_t @ Output: Concatenate the components. Insert brackets if the component is also a list. The components of the [[expr]] array, if any, should all be filled. <>= procedure :: to_string => prt_spec_list_to_string <>= recursive function prt_spec_list_to_string (object) result (string) class(prt_spec_list_t), intent(in) :: object type(string_t) :: string integer :: i string = "" if (allocated (object%expr)) then do i = 1, size (object%expr) if (i > 1) string = string // ", " select type (x => object%expr(i)%x) type is (prt_spec_list_t) string = string // "(" // x%to_string () // ")" class default string = string // x%to_string () end select end do end if end function prt_spec_list_to_string @ %def prt_spec_list_to_string @ Flatten: if there is a subexpression which is also a list, include the components as direct members of the current list. <>= procedure :: flatten => prt_spec_list_flatten <>= subroutine prt_spec_list_flatten (object) class(prt_spec_list_t), intent(inout) :: object type(prt_expr_t), dimension(:), allocatable :: tmp_expr integer :: i, n_flat, i_flat n_flat = 0 do i = 1, size (object%expr) select type (y => object%expr(i)%x) type is (prt_spec_list_t) n_flat = n_flat + size (y%expr) class default n_flat = n_flat + 1 end select end do if (n_flat > size (object%expr)) then allocate (tmp_expr (n_flat)) i_flat = 0 do i = 1, size (object%expr) select type (y => object%expr(i)%x) type is (prt_spec_list_t) tmp_expr (i_flat + 1 : i_flat + size (y%expr)) = y%expr i_flat = i_flat + size (y%expr) class default tmp_expr (i_flat + 1) = object%expr(i) i_flat = i_flat + 1 end select end do end if if (allocated (tmp_expr)) & call move_alloc (from = tmp_expr, to = object%expr) end subroutine prt_spec_list_flatten @ %def prt_spec_list_flatten @ Convert a list of sums into a sum of lists. (Subexpressions which are not sums are left untouched.) <>= subroutine distribute_prt_spec_list (object) class(prt_spec_expr_t), intent(inout), allocatable :: object class(prt_spec_expr_t), allocatable :: new_object integer, dimension(:), allocatable :: n, ii integer :: k, n_expr, n_terms, i_term select type (object) type is (prt_spec_list_t) n_expr = size (object%expr) allocate (n (n_expr), source = 1) allocate (ii (n_expr), source = 1) do k = 1, size (object%expr) select type (y => object%expr(k)%x) type is (prt_spec_sum_t) n(k) = size (y%expr) end select end do n_terms = product (n) if (n_terms > 1) then allocate (prt_spec_sum_t :: new_object) select type (new_object) type is (prt_spec_sum_t) allocate (new_object%expr (n_terms)) do i_term = 1, n_terms allocate (prt_spec_list_t :: new_object%expr(i_term)%x) select type (x => new_object%expr(i_term)%x) type is (prt_spec_list_t) allocate (x%expr (n_expr)) do k = 1, n_expr select type (y => object%expr(k)%x) type is (prt_spec_sum_t) x%expr(k) = y%expr(ii(k)) class default x%expr(k) = object%expr(k) end select end do end select INCR_INDEX: do k = n_expr, 1, -1 if (ii(k) < n(k)) then ii(k) = ii(k) + 1 exit INCR_INDEX else ii(k) = 1 end if end do INCR_INDEX end do end select end if end select if (allocated (new_object)) call move_alloc (from = new_object, to = object) end subroutine distribute_prt_spec_list @ %def distribute_prt_spec_list @ Apply [[expand]] to all components of the list. <>= procedure :: expand_sub => prt_spec_list_expand_sub <>= recursive subroutine prt_spec_list_expand_sub (object) class(prt_spec_list_t), intent(inout) :: object integer :: i if (allocated (object%expr)) then do i = 1, size (object%expr) call object%expr(i)%expand () end do end if end subroutine prt_spec_list_expand_sub @ %def prt_spec_list_expand_sub @ \subsection{Sum} A sum of particle specifiers, indicating, e.g., a sum of final states. <>= public :: prt_spec_sum_t <>= type, extends (prt_spec_expr_t) :: prt_spec_sum_t type(prt_expr_t), dimension(:), allocatable :: expr contains <> end type prt_spec_sum_t @ %def prt_spec_sum_t @ Output: Concatenate the components. Insert brackets if the component is a list or also a sum. The components of the [[expr]] array, if any, should all be filled. <>= procedure :: to_string => prt_spec_sum_to_string <>= recursive function prt_spec_sum_to_string (object) result (string) class(prt_spec_sum_t), intent(in) :: object type(string_t) :: string integer :: i string = "" if (allocated (object%expr)) then do i = 1, size (object%expr) if (i > 1) string = string // " + " select type (x => object%expr(i)%x) type is (prt_spec_list_t) string = string // "(" // x%to_string () // ")" type is (prt_spec_sum_t) string = string // "(" // x%to_string () // ")" class default string = string // x%to_string () end select end do end if end function prt_spec_sum_to_string @ %def prt_spec_sum_to_string @ Flatten: if there is a subexpression which is also a sum, include the components as direct members of the current sum. This is identical to [[prt_spec_list_flatten]] above, except for the type. <>= procedure :: flatten => prt_spec_sum_flatten <>= subroutine prt_spec_sum_flatten (object) class(prt_spec_sum_t), intent(inout) :: object type(prt_expr_t), dimension(:), allocatable :: tmp_expr integer :: i, n_flat, i_flat n_flat = 0 do i = 1, size (object%expr) select type (y => object%expr(i)%x) type is (prt_spec_sum_t) n_flat = n_flat + size (y%expr) class default n_flat = n_flat + 1 end select end do if (n_flat > size (object%expr)) then allocate (tmp_expr (n_flat)) i_flat = 0 do i = 1, size (object%expr) select type (y => object%expr(i)%x) type is (prt_spec_sum_t) tmp_expr (i_flat + 1 : i_flat + size (y%expr)) = y%expr i_flat = i_flat + size (y%expr) class default tmp_expr (i_flat + 1) = object%expr(i) i_flat = i_flat + 1 end select end do end if if (allocated (tmp_expr)) & call move_alloc (from = tmp_expr, to = object%expr) end subroutine prt_spec_sum_flatten @ %def prt_spec_sum_flatten @ Apply [[expand]] to all terms in the sum. <>= procedure :: expand_sub => prt_spec_sum_expand_sub <>= recursive subroutine prt_spec_sum_expand_sub (object) class(prt_spec_sum_t), intent(inout) :: object integer :: i if (allocated (object%expr)) then do i = 1, size (object%expr) call object%expr(i)%expand () end do end if end subroutine prt_spec_sum_expand_sub @ %def prt_spec_sum_expand_sub @ \subsection{Expression Expansion} The [[expand]] method transforms each particle specifier expression into a sum of lists, according to the rules \begin{align} a, (b, c) &\to a, b, c \\ a + (b + c) &\to a + b + c \\ a, b + c &\to (a, b) + (a, c) \end{align} Note that the precedence of comma and plus are opposite to this expansion, so the parentheses in the final expression are necessary. We assume that subexpressions are filled, i.e., arrays are allocated. <>= procedure :: expand => prt_expr_expand <>= recursive subroutine prt_expr_expand (expr) class(prt_expr_t), intent(inout) :: expr if (allocated (expr%x)) then call distribute_prt_spec_list (expr%x) call expr%x%expand_sub () select type (x => expr%x) type is (prt_spec_list_t) call x%flatten () type is (prt_spec_sum_t) call x%flatten () end select end if end subroutine prt_expr_expand @ %def prt_expr_expand @ \subsection{Unit Tests} Test module, followed by the corresponding implementation module. <<[[particle_specifiers_ut.f90]]>>= <> module particle_specifiers_ut use unit_tests use particle_specifiers_uti <> <> contains <> end module particle_specifiers_ut @ %def particle_specifiers_ut @ <<[[particle_specifiers_uti.f90]]>>= <> module particle_specifiers_uti <> use particle_specifiers <> <> contains <> end module particle_specifiers_uti @ %def particle_specifiers_ut @ API: driver for the unit tests below. <>= public :: particle_specifiers_test <>= subroutine particle_specifiers_test (u, results) integer, intent(in) :: u type(test_results_t), intent(inout) :: results <> end subroutine particle_specifiers_test @ %def particle_specifiers_test @ \subsubsection{Particle specifier array} Define, read and write an array of particle specifiers. <>= call test (particle_specifiers_1, "particle_specifiers_1", & "Handle particle specifiers", & u, results) <>= public :: particle_specifiers_1 <>= subroutine particle_specifiers_1 (u) integer, intent(in) :: u type(prt_spec_t), dimension(:), allocatable :: prt_spec type(string_t), dimension(:), allocatable :: decay type(string_t), dimension(0) :: no_decay integer :: i, j write (u, "(A)") "* Test output: particle_specifiers_1" write (u, "(A)") "* Purpose: Read and write a particle specifier array" write (u, "(A)") allocate (prt_spec (5)) prt_spec = [ & new_prt_spec (var_str ("a")), & new_prt_spec (var_str ("b"), .true.), & new_prt_spec (var_str ("c"), [var_str ("dec1")]), & new_prt_spec (var_str ("d"), [var_str ("dec1"), var_str ("dec2")]), & new_prt_spec (var_str ("e"), no_decay) & ] do i = 1, size (prt_spec) write (u, "(A)") char (prt_spec(i)%to_string ()) end do write (u, "(A)") call prt_spec_read (prt_spec, & var_str (" a, b( *), c( dec1), d (dec1 + dec2 ), e()")) call prt_spec_write (prt_spec, u) do i = 1, size (prt_spec) write (u, "(A)") write (u, "(A,A)") char (prt_spec(i)%get_name ()), ":" write (u, "(A,L1)") "polarized = ", prt_spec(i)%is_polarized () write (u, "(A,L1)") "unstable = ", prt_spec(i)%is_unstable () write (u, "(A,I0)") "n_decays = ", prt_spec(i)%get_n_decays () call prt_spec(i)%get_decays (decay) write (u, "(A)", advance="no") "decays =" do j = 1, size (decay) write (u, "(1x,A)", advance="no") char (decay(j)) end do write (u, "(A)") end do write (u, "(A)") write (u, "(A)") "* Test output end: particle_specifiers_1" end subroutine particle_specifiers_1 @ %def particle_specifiers_1 @ \subsubsection{Particle specifier expressions} Nested expressions (only basic particles, no decay specs). <>= call test (particle_specifiers_2, "particle_specifiers_2", & "Particle specifier expressions", & u, results) <>= public :: particle_specifiers_2 <>= subroutine particle_specifiers_2 (u) integer, intent(in) :: u type(prt_spec_t) :: a, b, c, d, e, f type(prt_expr_t) :: pe1, pe2, pe3 type(prt_expr_t) :: pe4, pe5, pe6, pe7, pe8, pe9 integer :: i type(prt_spec_t), dimension(:), allocatable :: pa write (u, "(A)") "* Test output: particle_specifiers_2" write (u, "(A)") "* Purpose: Create and display particle expressions" write (u, "(A)") write (u, "(A)") "* Basic expressions" write (u, *) a = new_prt_spec (var_str ("a")) b = new_prt_spec (var_str ("b")) c = new_prt_spec (var_str ("c")) d = new_prt_spec (var_str ("d")) e = new_prt_spec (var_str ("e")) f = new_prt_spec (var_str ("f")) call pe1%init_spec (a) write (u, "(A)") char (pe1%to_string ()) call pe2%init_sum (2) select type (x => pe2%x) type is (prt_spec_sum_t) call x%expr(1)%init_spec (a) call x%expr(2)%init_spec (b) end select write (u, "(A)") char (pe2%to_string ()) call pe3%init_list (2) select type (x => pe3%x) type is (prt_spec_list_t) call x%expr(1)%init_spec (a) call x%expr(2)%init_spec (b) end select write (u, "(A)") char (pe3%to_string ()) write (u, *) write (u, "(A)") "* Nested expressions" write (u, *) call pe4%init_list (2) select type (x => pe4%x) type is (prt_spec_list_t) call x%expr(1)%init_sum (2) select type (y => x%expr(1)%x) type is (prt_spec_sum_t) call y%expr(1)%init_spec (a) call y%expr(2)%init_spec (b) end select call x%expr(2)%init_spec (c) end select write (u, "(A)") char (pe4%to_string ()) call pe5%init_list (2) select type (x => pe5%x) type is (prt_spec_list_t) call x%expr(1)%init_list (2) select type (y => x%expr(1)%x) type is (prt_spec_list_t) call y%expr(1)%init_spec (a) call y%expr(2)%init_spec (b) end select call x%expr(2)%init_spec (c) end select write (u, "(A)") char (pe5%to_string ()) call pe6%init_sum (2) select type (x => pe6%x) type is (prt_spec_sum_t) call x%expr(1)%init_spec (a) call x%expr(2)%init_sum (2) select type (y => x%expr(2)%x) type is (prt_spec_sum_t) call y%expr(1)%init_spec (b) call y%expr(2)%init_spec (c) end select end select write (u, "(A)") char (pe6%to_string ()) call pe7%init_list (2) select type (x => pe7%x) type is (prt_spec_list_t) call x%expr(1)%init_sum (2) select type (y => x%expr(1)%x) type is (prt_spec_sum_t) call y%expr(1)%init_spec (a) call y%expr(2)%init_list (2) select type (z => y%expr(2)%x) type is (prt_spec_list_t) call z%expr(1)%init_spec (b) call z%expr(2)%init_spec (c) end select end select call x%expr(2)%init_spec (d) end select write (u, "(A)") char (pe7%to_string ()) call pe8%init_sum (2) select type (x => pe8%x) type is (prt_spec_sum_t) call x%expr(1)%init_list (2) select type (y => x%expr(1)%x) type is (prt_spec_list_t) call y%expr(1)%init_spec (a) call y%expr(2)%init_spec (b) end select call x%expr(2)%init_list (2) select type (y => x%expr(2)%x) type is (prt_spec_list_t) call y%expr(1)%init_spec (c) call y%expr(2)%init_spec (d) end select end select write (u, "(A)") char (pe8%to_string ()) call pe9%init_list (3) select type (x => pe9%x) type is (prt_spec_list_t) call x%expr(1)%init_sum (2) select type (y => x%expr(1)%x) type is (prt_spec_sum_t) call y%expr(1)%init_spec (a) call y%expr(2)%init_spec (b) end select call x%expr(2)%init_spec (c) call x%expr(3)%init_sum (3) select type (y => x%expr(3)%x) type is (prt_spec_sum_t) call y%expr(1)%init_spec (d) call y%expr(2)%init_spec (e) call y%expr(3)%init_spec (f) end select end select write (u, "(A)") char (pe9%to_string ()) write (u, *) write (u, "(A)") "* Expand as sum" write (u, *) call pe1%expand () write (u, "(A)") char (pe1%to_string ()) call pe4%expand () write (u, "(A)") char (pe4%to_string ()) call pe5%expand () write (u, "(A)") char (pe5%to_string ()) call pe6%expand () write (u, "(A)") char (pe6%to_string ()) call pe7%expand () write (u, "(A)") char (pe7%to_string ()) call pe8%expand () write (u, "(A)") char (pe8%to_string ()) call pe9%expand () write (u, "(A)") char (pe9%to_string ()) write (u, *) write (u, "(A)") "* Transform to arrays:" write (u, "(A)") "* Atomic specifier" do i = 1, pe1%get_n_terms () call pe1%term_to_array (pa, i) call prt_spec_write (pa, u) end do write (u, *) write (u, "(A)") "* List" do i = 1, pe5%get_n_terms () call pe5%term_to_array (pa, i) call prt_spec_write (pa, u) end do write (u, *) write (u, "(A)") "* Sum of atoms" do i = 1, pe6%get_n_terms () call pe6%term_to_array (pa, i) call prt_spec_write (pa, u) end do write (u, *) write (u, "(A)") "* Sum of lists" do i = 1, pe9%get_n_terms () call pe9%term_to_array (pa, i) call prt_spec_write (pa, u) end do write (u, "(A)") write (u, "(A)") "* Test output end: particle_specifiers_2" end subroutine particle_specifiers_2 @ %def particle_specifiers_2 @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{PDG arrays} For defining aliases, we introduce a special type which holds a set of (integer) PDG codes. <<[[pdg_arrays.f90]]>>= <> module pdg_arrays use io_units use sorting use physics_defs <> <> <> <> contains <> end module pdg_arrays @ %def pdg_arrays @ \subsection{Type definition} Using an allocatable array eliminates the need for initializer and/or finalizer. <>= public :: pdg_array_t <>= type :: pdg_array_t private integer, dimension(:), allocatable :: pdg contains <> end type pdg_array_t @ %def pdg_array_t @ Output <>= public :: pdg_array_write <>= procedure :: write => pdg_array_write <>= subroutine pdg_array_write (aval, unit) class(pdg_array_t), intent(in) :: aval integer, intent(in), optional :: unit integer :: u, i u = given_output_unit (unit); if (u < 0) return write (u, "(A)", advance="no") "PDG(" if (allocated (aval%pdg)) then do i = 1, size (aval%pdg) if (i > 1) write (u, "(A)", advance="no") ", " write (u, "(I0)", advance="no") aval%pdg(i) end do end if write (u, "(A)", advance="no") ")" end subroutine pdg_array_write @ %def pdg_array_write @ <>= public :: pdg_array_write_set <>= subroutine pdg_array_write_set (aval, unit) type(pdg_array_t), intent(in), dimension(:) :: aval integer, intent(in), optional :: unit integer :: i do i = 1, size (aval) call aval(i)%write (unit) print *, '' end do end subroutine pdg_array_write_set @ %def pdg_array_write_set @ \subsection{Basic operations} Assignment. We define assignment from and to an integer array. Note that the integer array, if it is the l.h.s., must be declared allocatable by the caller. <>= public :: assignment(=) <>= interface assignment(=) module procedure pdg_array_from_int_array module procedure pdg_array_from_int module procedure int_array_from_pdg_array end interface <>= subroutine pdg_array_from_int_array (aval, iarray) type(pdg_array_t), intent(out) :: aval integer, dimension(:), intent(in) :: iarray allocate (aval%pdg (size (iarray))) aval%pdg = iarray end subroutine pdg_array_from_int_array elemental subroutine pdg_array_from_int (aval, int) type(pdg_array_t), intent(out) :: aval integer, intent(in) :: int allocate (aval%pdg (1)) aval%pdg = int end subroutine pdg_array_from_int subroutine int_array_from_pdg_array (iarray, aval) integer, dimension(:), allocatable, intent(out) :: iarray type(pdg_array_t), intent(in) :: aval if (allocated (aval%pdg)) then allocate (iarray (size (aval%pdg))) iarray = aval%pdg else allocate (iarray (0)) end if end subroutine int_array_from_pdg_array @ %def pdg_array_from_int_array pdg_array_from_int int_array_from_pdg_array @ Allocate space for a PDG array <>= public :: pdg_array_init <>= subroutine pdg_array_init (aval, n_elements) type(pdg_array_t), intent(inout) :: aval integer, intent(in) :: n_elements allocate(aval%pdg(n_elements)) end subroutine pdg_array_init @ %def pdg_array_init @ Deallocate a previously allocated pdg array <>= public :: pdg_array_delete <>= subroutine pdg_array_delete (aval) type(pdg_array_t), intent(inout) :: aval if (allocated (aval%pdg)) deallocate (aval%pdg) end subroutine pdg_array_delete @ %def pdg_array_delete @ Merge two pdg arrays, i.e. append a particle string to another leaving out doublettes <>= public :: pdg_array_merge <>= subroutine pdg_array_merge (aval1, aval2) type(pdg_array_t), intent(inout) :: aval1 type(pdg_array_t), intent(in) :: aval2 type(pdg_array_t) :: aval if (allocated (aval1%pdg) .and. allocated (aval2%pdg)) then if (.not. any (aval1%pdg == aval2%pdg)) aval = aval1 // aval2 else if (allocated (aval1%pdg)) then aval = aval1 else if (allocated (aval2%pdg)) then aval = aval2 end if call pdg_array_delete (aval1) aval1 = aval%pdg end subroutine pdg_array_merge @ %def pdg_array_merge @ Length of the array. <>= public :: pdg_array_get_length <>= procedure :: get_length => pdg_array_get_length <>= elemental function pdg_array_get_length (aval) result (n) class(pdg_array_t), intent(in) :: aval integer :: n if (allocated (aval%pdg)) then n = size (aval%pdg) else n = 0 end if end function pdg_array_get_length @ %def pdg_array_get_length @ Return the element with index i. <>= public :: pdg_array_get <>= procedure :: get => pdg_array_get <>= elemental function pdg_array_get (aval, i) result (pdg) class(pdg_array_t), intent(in) :: aval integer, intent(in), optional :: i integer :: pdg if (present (i)) then pdg = aval%pdg(i) else pdg = aval%pdg(1) end if end function pdg_array_get @ %def pdg_array_get @ Explicitly set the element with index i. <>= procedure :: set => pdg_array_set <>= subroutine pdg_array_set (aval, i, pdg) class(pdg_array_t), intent(inout) :: aval integer, intent(in) :: i integer, intent(in) :: pdg aval%pdg(i) = pdg end subroutine pdg_array_set @ %def pdg_array_set @ <>= procedure :: add => pdg_array_add <>= function pdg_array_add (aval, aval_add) result (aval_out) type(pdg_array_t) :: aval_out class(pdg_array_t), intent(in) :: aval type(pdg_array_t), intent(in) :: aval_add integer :: n, n_add, i n = size (aval%pdg) n_add = size (aval_add%pdg) allocate (aval_out%pdg (n + n_add)) aval_out%pdg(1:n) = aval%pdg do i = 1, n_add aval_out%pdg(n+i) = aval_add%pdg(i) end do end function pdg_array_add @ %def pdg_array_add @ Replace element with index [[i]] by a new array of elements. <>= public :: pdg_array_replace <>= procedure :: replace => pdg_array_replace <>= function pdg_array_replace (aval, i, pdg_new) result (aval_new) class(pdg_array_t), intent(in) :: aval integer, intent(in) :: i integer, dimension(:), intent(in) :: pdg_new type(pdg_array_t) :: aval_new integer :: n, l n = size (aval%pdg) l = size (pdg_new) allocate (aval_new%pdg (n + l - 1)) aval_new%pdg(:i-1) = aval%pdg(:i-1) aval_new%pdg(i:i+l-1) = pdg_new aval_new%pdg(i+l:) = aval%pdg(i+1:) end function pdg_array_replace @ %def pdg_array_replace @ Concatenate two PDG arrays <>= public :: operator(//) <>= interface operator(//) module procedure concat_pdg_arrays end interface <>= function concat_pdg_arrays (aval1, aval2) result (aval) type(pdg_array_t) :: aval type(pdg_array_t), intent(in) :: aval1, aval2 integer :: n1, n2 if (allocated (aval1%pdg) .and. allocated (aval2%pdg)) then n1 = size (aval1%pdg) n2 = size (aval2%pdg) allocate (aval%pdg (n1 + n2)) aval%pdg(:n1) = aval1%pdg aval%pdg(n1+1:) = aval2%pdg else if (allocated (aval1%pdg)) then aval = aval1 else if (allocated (aval2%pdg)) then aval = aval2 end if end function concat_pdg_arrays @ %def concat_pdg_arrays @ \subsection{Matching} A PDG array matches a given PDG code if the code is present within the array. If either one is zero (UNDEFINED), the match also succeeds. <>= public :: operator(.match.) <>= interface operator(.match.) module procedure pdg_array_match_integer module procedure pdg_array_match_pdg_array end interface @ %def .match. @ Match a single code against the array. <>= elemental function pdg_array_match_integer (aval, pdg) result (flag) logical :: flag type(pdg_array_t), intent(in) :: aval integer, intent(in) :: pdg if (allocated (aval%pdg)) then flag = pdg == UNDEFINED & .or. any (aval%pdg == UNDEFINED) & .or. any (aval%pdg == pdg) else flag = .false. end if end function pdg_array_match_integer @ %def pdg_array_match_integer @ Check if the pdg-number corresponds to a quark <>= public :: is_quark <>= elemental function is_quark (pdg_nr) logical :: is_quark integer, intent(in) :: pdg_nr if (abs (pdg_nr) >= 1 .and. abs (pdg_nr) <= 6) then is_quark = .true. else is_quark = .false. end if end function is_quark @ %def is_quark @ Check if pdg-number corresponds to a gluon <>= public :: is_gluon <>= elemental function is_gluon (pdg_nr) logical :: is_gluon integer, intent(in) :: pdg_nr if (pdg_nr == GLUON) then is_gluon = .true. else is_gluon = .false. end if end function is_gluon @ %def is_gluon @ Check if pdg-number corresponds to a photon <>= public :: is_photon <>= elemental function is_photon (pdg_nr) logical :: is_photon integer, intent(in) :: pdg_nr if (pdg_nr == PHOTON) then is_photon = .true. else is_photon = .false. end if end function is_photon @ %def is_photon @ Check if pdg-number corresponds to a colored particle <>= public :: is_colored <>= elemental function is_colored (pdg_nr) logical :: is_colored integer, intent(in) :: pdg_nr is_colored = is_quark (pdg_nr) .or. is_gluon (pdg_nr) end function is_colored @ %def is_colored @ Check if the pdg-number corresponds to a lepton <>= public :: is_lepton <>= elemental function is_lepton (pdg_nr) logical :: is_lepton integer, intent(in) :: pdg_nr if (abs (pdg_nr) >= ELECTRON .and. & abs (pdg_nr) <= TAU_NEUTRINO) then is_lepton = .true. else is_lepton = .false. end if end function is_lepton @ %def is_lepton @ <>= public :: is_fermion <>= elemental function is_fermion (pdg_nr) logical :: is_fermion integer, intent(in) :: pdg_nr is_fermion = is_lepton(pdg_nr) .or. is_quark(pdg_nr) end function is_fermion @ %def is_fermion @ Check if the pdg-number corresponds to a massless vector boson <>= public :: is_massless_vector <>= elemental function is_massless_vector (pdg_nr) integer, intent(in) :: pdg_nr logical :: is_massless_vector if (pdg_nr == GLUON .or. pdg_nr == PHOTON) then is_massless_vector = .true. else is_massless_vector = .false. end if end function is_massless_vector @ %def is_massless_vector @ Check if pdg-number corresponds to a massive vector boson <>= public :: is_massive_vector <>= elemental function is_massive_vector (pdg_nr) integer, intent(in) :: pdg_nr logical :: is_massive_vector if (abs (pdg_nr) == Z_BOSON .or. abs (pdg_nr) == W_BOSON) then is_massive_vector = .true. else is_massive_vector = .false. end if end function is_massive_vector @ %def is massive_vector @ Check if pdg-number corresponds to a vector boson <>= public :: is_vector <>= elemental function is_vector (pdg_nr) integer, intent(in) :: pdg_nr logical :: is_vector if (is_massless_vector (pdg_nr) .or. is_massive_vector (pdg_nr)) then is_vector = .true. else is_vector = .false. end if end function is_vector @ %def is vector @ Check if particle is elementary. <>= public :: is_elementary <>= elemental function is_elementary (pdg_nr) integer, intent(in) :: pdg_nr logical :: is_elementary if (is_vector (pdg_nr) .or. is_fermion (pdg_nr) .or. pdg_nr == 25) then is_elementary = .true. else is_elementary = .false. end if end function is_elementary @ %def is_elementary @ Check if particle is strongly interacting <>= procedure :: has_colored_particles => pdg_array_has_colored_particles <>= function pdg_array_has_colored_particles (pdg) result (colored) class(pdg_array_t), intent(in) :: pdg logical :: colored integer :: i, pdg_nr colored = .false. do i = 1, size (pdg%pdg) pdg_nr = pdg%pdg(i) if (is_quark (pdg_nr) .or. is_gluon (pdg_nr)) then colored = .true. exit end if end do end function pdg_array_has_colored_particles @ %def pdg_array_has_colored_particles @ Match two arrays. Succeeds if any pair of entries matches. <>= function pdg_array_match_pdg_array (aval1, aval2) result (flag) logical :: flag type(pdg_array_t), intent(in) :: aval1, aval2 if (allocated (aval1%pdg) .and. allocated (aval2%pdg)) then flag = any (aval1 .match. aval2%pdg) else flag = .false. end if end function pdg_array_match_pdg_array @ %def pdg_array_match_pdg_array @ Comparison. Here, we take the PDG arrays as-is, assuming that they are sorted. The ordering is a bit odd: first, we look only at the absolute values of the PDG codes. If they all match, the particle comes before the antiparticle, scanning from left to right. <>= public :: operator(<) public :: operator(>) public :: operator(<=) public :: operator(>=) public :: operator(==) public :: operator(/=) <>= interface operator(<) module procedure pdg_array_lt end interface interface operator(>) module procedure pdg_array_gt end interface interface operator(<=) module procedure pdg_array_le end interface interface operator(>=) module procedure pdg_array_ge end interface interface operator(==) module procedure pdg_array_eq end interface interface operator(/=) module procedure pdg_array_ne end interface <>= elemental function pdg_array_lt (aval1, aval2) result (flag) type(pdg_array_t), intent(in) :: aval1, aval2 logical :: flag integer :: i if (size (aval1%pdg) /= size (aval2%pdg)) then flag = size (aval1%pdg) < size (aval2%pdg) else do i = 1, size (aval1%pdg) if (abs (aval1%pdg(i)) /= abs (aval2%pdg(i))) then flag = abs (aval1%pdg(i)) < abs (aval2%pdg(i)) return end if end do do i = 1, size (aval1%pdg) if (aval1%pdg(i) /= aval2%pdg(i)) then flag = aval1%pdg(i) > aval2%pdg(i) return end if end do flag = .false. end if end function pdg_array_lt elemental function pdg_array_gt (aval1, aval2) result (flag) type(pdg_array_t), intent(in) :: aval1, aval2 logical :: flag flag = .not. (aval1 < aval2 .or. aval1 == aval2) end function pdg_array_gt elemental function pdg_array_le (aval1, aval2) result (flag) type(pdg_array_t), intent(in) :: aval1, aval2 logical :: flag flag = aval1 < aval2 .or. aval1 == aval2 end function pdg_array_le elemental function pdg_array_ge (aval1, aval2) result (flag) type(pdg_array_t), intent(in) :: aval1, aval2 logical :: flag flag = .not. (aval1 < aval2) end function pdg_array_ge elemental function pdg_array_eq (aval1, aval2) result (flag) type(pdg_array_t), intent(in) :: aval1, aval2 logical :: flag if (size (aval1%pdg) /= size (aval2%pdg)) then flag = .false. else flag = all (aval1%pdg == aval2%pdg) end if end function pdg_array_eq elemental function pdg_array_ne (aval1, aval2) result (flag) type(pdg_array_t), intent(in) :: aval1, aval2 logical :: flag flag = .not. (aval1 == aval2) end function pdg_array_ne @ Equivalence. Two PDG arrays are equivalent if either one contains [[UNDEFINED]] or if each element of array 1 is present in array 2, and vice versa. <>= public :: operator(.eqv.) public :: operator(.neqv.) <>= interface operator(.eqv.) module procedure pdg_array_equivalent end interface interface operator(.neqv.) module procedure pdg_array_inequivalent end interface <>= elemental function pdg_array_equivalent (aval1, aval2) result (eq) logical :: eq type(pdg_array_t), intent(in) :: aval1, aval2 logical, dimension(:), allocatable :: match1, match2 integer :: i if (allocated (aval1%pdg) .and. allocated (aval2%pdg)) then eq = any (aval1%pdg == UNDEFINED) & .or. any (aval2%pdg == UNDEFINED) if (.not. eq) then allocate (match1 (size (aval1%pdg))) allocate (match2 (size (aval2%pdg))) match1 = .false. match2 = .false. do i = 1, size (aval1%pdg) match2 = match2 .or. aval1%pdg(i) == aval2%pdg end do do i = 1, size (aval2%pdg) match1 = match1 .or. aval2%pdg(i) == aval1%pdg end do eq = all (match1) .and. all (match2) end if else eq = .false. end if end function pdg_array_equivalent elemental function pdg_array_inequivalent (aval1, aval2) result (neq) logical :: neq type(pdg_array_t), intent(in) :: aval1, aval2 neq = .not. pdg_array_equivalent (aval1, aval2) end function pdg_array_inequivalent @ %def pdg_array_equivalent @ \subsection{Sorting} Sort a PDG array by absolute value, particle before antiparticle. After sorting, we eliminate double entries. <>= public :: sort_abs <>= interface sort_abs module procedure pdg_array_sort_abs end interface <>= procedure :: sort_abs => pdg_array_sort_abs <>= function pdg_array_sort_abs (aval1, unique) result (aval2) class(pdg_array_t), intent(in) :: aval1 logical, intent(in), optional :: unique type(pdg_array_t) :: aval2 integer, dimension(:), allocatable :: tmp logical, dimension(:), allocatable :: mask integer :: i, n logical :: uni uni = .false.; if (present (unique)) uni = unique n = size (aval1%pdg) if (uni) then allocate (tmp (n), mask(n)) tmp = sort_abs (aval1%pdg) mask(1) = .true. do i = 2, n mask(i) = tmp(i) /= tmp(i-1) end do allocate (aval2%pdg (count (mask))) aval2%pdg = pack (tmp, mask) else allocate (aval2%pdg (n)) aval2%pdg = sort_abs (aval1%pdg) end if end function pdg_array_sort_abs @ %def sort_abs @ <>= procedure :: intersect => pdg_array_intersect <>= function pdg_array_intersect (aval1, match) result (aval2) class(pdg_array_t), intent(in) :: aval1 integer, dimension(:) :: match type(pdg_array_t) :: aval2 integer, dimension(:), allocatable :: isec integer :: i isec = pack (aval1%pdg, [(any(aval1%pdg(i) == match), i=1,size(aval1%pdg))]) aval2 = isec end function pdg_array_intersect @ %def pdg_array_intersect @ <>= procedure :: search_for_particle => pdg_array_search_for_particle <>= elemental function pdg_array_search_for_particle (pdg, i_part) result (found) class(pdg_array_t), intent(in) :: pdg integer, intent(in) :: i_part logical :: found found = any (pdg%pdg == i_part) end function pdg_array_search_for_particle @ %def pdg_array_search_for_particle @ <>= procedure :: invert => pdg_array_invert <>= function pdg_array_invert (pdg) result (pdg_inverse) class(pdg_array_t), intent(in) :: pdg type(pdg_array_t) :: pdg_inverse integer :: i, n n = size (pdg%pdg) allocate (pdg_inverse%pdg (n)) do i = 1, n select case (pdg%pdg(i)) case (GLUON, PHOTON, Z_BOSON, 25) pdg_inverse%pdg(i) = pdg%pdg(i) case default pdg_inverse%pdg(i) = -pdg%pdg(i) end select end do end function pdg_array_invert @ %def pdg_array_invert @ \subsection{PDG array list} A PDG array list, or PDG list, is an array of PDG-array objects with some convenience methods. <>= public :: pdg_list_t <>= type :: pdg_list_t type(pdg_array_t), dimension(:), allocatable :: a contains <> end type pdg_list_t @ %def pdg_list_t @ Output, as a comma-separated list without advancing I/O. <>= procedure :: write => pdg_list_write <>= subroutine pdg_list_write (object, unit) class(pdg_list_t), intent(in) :: object integer, intent(in), optional :: unit integer :: u, i u = given_output_unit (unit) if (allocated (object%a)) then do i = 1, size (object%a) if (i > 1) write (u, "(A)", advance="no") ", " call object%a(i)%write (u) end do end if end subroutine pdg_list_write @ %def pdg_list_write @ Initialize for a certain size. The entries are initially empty PDG arrays. <>= generic :: init => pdg_list_init_size procedure, private :: pdg_list_init_size <>= subroutine pdg_list_init_size (pl, n) class(pdg_list_t), intent(out) :: pl integer, intent(in) :: n allocate (pl%a (n)) end subroutine pdg_list_init_size @ %def pdg_list_init_size @ Initialize with a definite array of PDG codes. That is, each entry in the list becomes a single-particle PDG array. <>= generic :: init => pdg_list_init_int_array procedure, private :: pdg_list_init_int_array <>= subroutine pdg_list_init_int_array (pl, pdg) class(pdg_list_t), intent(out) :: pl integer, dimension(:), intent(in) :: pdg integer :: i allocate (pl%a (size (pdg))) do i = 1, size (pdg) pl%a(i) = pdg(i) end do end subroutine pdg_list_init_int_array @ %def pdg_list_init_array @ Set one of the entries. No bounds-check. <>= generic :: set => pdg_list_set_int generic :: set => pdg_list_set_int_array generic :: set => pdg_list_set_pdg_array procedure, private :: pdg_list_set_int procedure, private :: pdg_list_set_int_array procedure, private :: pdg_list_set_pdg_array <>= subroutine pdg_list_set_int (pl, i, pdg) class(pdg_list_t), intent(inout) :: pl integer, intent(in) :: i integer, intent(in) :: pdg pl%a(i) = pdg end subroutine pdg_list_set_int subroutine pdg_list_set_int_array (pl, i, pdg) class(pdg_list_t), intent(inout) :: pl integer, intent(in) :: i integer, dimension(:), intent(in) :: pdg pl%a(i) = pdg end subroutine pdg_list_set_int_array subroutine pdg_list_set_pdg_array (pl, i, pa) class(pdg_list_t), intent(inout) :: pl integer, intent(in) :: i type(pdg_array_t), intent(in) :: pa pl%a(i) = pa end subroutine pdg_list_set_pdg_array @ %def pdg_list_set @ Array size, not the length of individual entries <>= procedure :: get_size => pdg_list_get_size <>= function pdg_list_get_size (pl) result (n) class(pdg_list_t), intent(in) :: pl integer :: n if (allocated (pl%a)) then n = size (pl%a) else n = 0 end if end function pdg_list_get_size @ %def pdg_list_get_size @ Return an entry, as a PDG array. <>= procedure :: get => pdg_list_get <>= function pdg_list_get (pl, i) result (pa) type(pdg_array_t) :: pa class(pdg_list_t), intent(in) :: pl integer, intent(in) :: i pa = pl%a(i) end function pdg_list_get @ %def pdg_list_get @ Check if the list entries are all either mutually disjoint or identical. The individual entries (PDG arrays) should already be sorted, so we can test for equality. <>= procedure :: is_regular => pdg_list_is_regular <>= function pdg_list_is_regular (pl) result (flag) class(pdg_list_t), intent(in) :: pl logical :: flag integer :: i, j, s s = pl%get_size () flag = .true. do i = 1, s do j = i + 1, s if (pl%a(i) .match. pl%a(j)) then if (pl%a(i) /= pl%a(j)) then flag = .false. return end if end if end do end do end function pdg_list_is_regular @ %def pdg_list_is_regular @ Sort the list. First, each entry gets sorted, including elimination of doublers. Then, we sort the list, using the first member of each PDG array as the marker. No removal of doublers at this stage. If [[n_in]] is supplied, we do not reorder the first [[n_in]] particle entries. <>= procedure :: sort_abs => pdg_list_sort_abs <>= function pdg_list_sort_abs (pl, n_in) result (pl_sorted) class(pdg_list_t), intent(in) :: pl integer, intent(in), optional :: n_in type(pdg_list_t) :: pl_sorted type(pdg_array_t), dimension(:), allocatable :: pa integer, dimension(:), allocatable :: pdg, map integer :: i, n0 call pl_sorted%init (pl%get_size ()) if (allocated (pl%a)) then allocate (pa (size (pl%a))) do i = 1, size (pl%a) pa(i) = pl%a(i)%sort_abs (unique = .true.) end do allocate (pdg (size (pa)), source = 0) do i = 1, size (pa) if (allocated (pa(i)%pdg)) then if (size (pa(i)%pdg) > 0) then pdg(i) = pa(i)%pdg(1) end if end if end do if (present (n_in)) then n0 = n_in else n0 = 0 end if allocate (map (size (pdg))) map(:n0) = [(i, i = 1, n0)] map(n0+1:) = n0 + order_abs (pdg(n0+1:)) do i = 1, size (pa) call pl_sorted%set (i, pa(map(i))) end do end if end function pdg_list_sort_abs @ %def pdg_list_sort_abs @ Compare sorted lists: equality. The result is undefined if some entries are not allocated. <>= generic :: operator (==) => pdg_list_eq procedure, private :: pdg_list_eq <>= function pdg_list_eq (pl1, pl2) result (flag) class(pdg_list_t), intent(in) :: pl1, pl2 logical :: flag integer :: i flag = .false. if (allocated (pl1%a) .and. allocated (pl2%a)) then if (size (pl1%a) == size (pl2%a)) then do i = 1, size (pl1%a) associate (a1 => pl1%a(i), a2 => pl2%a(i)) if (allocated (a1%pdg) .and. allocated (a2%pdg)) then if (size (a1%pdg) == size (a2%pdg)) then if (size (a1%pdg) > 0) then if (a1%pdg(1) /= a2%pdg(1)) return end if else return end if else return end if end associate end do flag = .true. end if end if end function pdg_list_eq @ %def pdg_list_eq @ Compare sorted lists. The result is undefined if some entries are not allocated. The ordering is quite complicated. First, a shorter list comes before a longer list. Comparing entry by entry, a shorter entry comes first. Next, we check the first PDG code within corresponding entries. This is compared by absolute value. If equal, particle comes before antiparticle. Finally, if all is equal, the result is false. <>= generic :: operator (<) => pdg_list_lt procedure, private :: pdg_list_lt <>= function pdg_list_lt (pl1, pl2) result (flag) class(pdg_list_t), intent(in) :: pl1, pl2 logical :: flag integer :: i flag = .false. if (allocated (pl1%a) .and. allocated (pl2%a)) then if (size (pl1%a) < size (pl2%a)) then flag = .true.; return else if (size (pl1%a) > size (pl2%a)) then return else do i = 1, size (pl1%a) associate (a1 => pl1%a(i), a2 => pl2%a(i)) if (allocated (a1%pdg) .and. allocated (a2%pdg)) then if (size (a1%pdg) < size (a2%pdg)) then flag = .true.; return else if (size (a1%pdg) > size (a2%pdg)) then return else if (size (a1%pdg) > 0) then if (abs (a1%pdg(1)) < abs (a2%pdg(1))) then flag = .true.; return else if (abs (a1%pdg(1)) > abs (a2%pdg(1))) then return else if (a1%pdg(1) > 0 .and. a2%pdg(1) < 0) then flag = .true.; return else if (a1%pdg(1) < 0 .and. a2%pdg(1) > 0) then return end if end if end if else return end if end associate end do flag = .false. end if end if end function pdg_list_lt @ %def pdg_list_lt @ Replace an entry. In the result, the entry [[#i]] is replaced by the contents of the second argument. The result is not sorted. If [[n_in]] is also set and [[i]] is less or equal to [[n_in]], replace [[#i]] only by the first entry of [[pl_insert]], and insert the remainder after entry [[n_in]]. <>= procedure :: replace => pdg_list_replace <>= function pdg_list_replace (pl, i, pl_insert, n_in) result (pl_out) type(pdg_list_t) :: pl_out class(pdg_list_t), intent(in) :: pl integer, intent(in) :: i class(pdg_list_t), intent(in) :: pl_insert integer, intent(in), optional :: n_in integer :: n, n_insert, n_out, k n = pl%get_size () n_insert = pl_insert%get_size () n_out = n + n_insert - 1 call pl_out%init (n_out) ! if (allocated (pl%a)) then do k = 1, i - 1 pl_out%a(k) = pl%a(k) end do ! end if if (present (n_in)) then pl_out%a(i) = pl_insert%a(1) do k = i + 1, n_in pl_out%a(k) = pl%a(k) end do do k = 1, n_insert - 1 pl_out%a(n_in+k) = pl_insert%a(1+k) end do do k = 1, n - n_in pl_out%a(n_in+k+n_insert-1) = pl%a(n_in+k) end do else ! if (allocated (pl_insert%a)) then do k = 1, n_insert pl_out%a(i-1+k) = pl_insert%a(k) end do ! end if ! if (allocated (pl%a)) then do k = 1, n - i pl_out%a(i+n_insert-1+k) = pl%a(i+k) end do end if ! end if end function pdg_list_replace @ %def pdg_list_replace @ <>= procedure :: fusion => pdg_list_fusion <>= function pdg_list_fusion (pl, pl_insert, i, check_if_existing) result (pl_out) type(pdg_list_t) :: pl_out class(pdg_list_t), intent(in) :: pl type(pdg_list_t), intent(in) :: pl_insert integer, intent(in) :: i logical, intent(in) :: check_if_existing integer :: n, n_insert, k, n_out logical :: new_pdg n = pl%get_size () n_insert = pl_insert%get_size () new_pdg = .not. check_if_existing .or. & (.not. any (pl%search_for_particle (pl_insert%a(1)%pdg))) call pl_out%init (n + n_insert - 1) do k = 1, n if (new_pdg .and. k == i) then pl_out%a(k) = pl%a(k)%add (pl_insert%a(1)) else pl_out%a(k) = pl%a(k) end if end do do k = n + 1, n + n_insert - 1 pl_out%a(k) = pl_insert%a(k-n) end do end function pdg_list_fusion @ %def pdg_list_fusion @ <>= procedure :: get_pdg_sizes => pdg_list_get_pdg_sizes <>= function pdg_list_get_pdg_sizes (pl) result (i_size) integer, dimension(:), allocatable :: i_size class(pdg_list_t), intent(in) :: pl integer :: i, n n = pl%get_size () allocate (i_size (n)) do i = 1, n i_size(i) = size (pl%a(i)%pdg) end do end function pdg_list_get_pdg_sizes @ %def pdg_list_get_pdg_sizes @ Replace the entries of [[pl]] by the matching entries of [[pl_match]], one by one. This is done in-place. If there is no match, return failure. <>= procedure :: match_replace => pdg_list_match_replace <>= subroutine pdg_list_match_replace (pl, pl_match, success) class(pdg_list_t), intent(inout) :: pl class(pdg_list_t), intent(in) :: pl_match logical, intent(out) :: success integer :: i, j success = .true. SCAN_ENTRIES: do i = 1, size (pl%a) do j = 1, size (pl_match%a) if (pl%a(i) .match. pl_match%a(j)) then pl%a(i) = pl_match%a(j) cycle SCAN_ENTRIES end if end do success = .false. return end do SCAN_ENTRIES end subroutine pdg_list_match_replace @ %def pdg_list_match_replace @ Just check if a PDG array matches any entry in the PDG list. The second version returns the position of the match within the list. An optional mask indicates the list elements that should be checked. <>= generic :: operator (.match.) => pdg_list_match_pdg_array procedure, private :: pdg_list_match_pdg_array procedure :: find_match => pdg_list_find_match_pdg_array <>= function pdg_list_match_pdg_array (pl, pa) result (flag) class(pdg_list_t), intent(in) :: pl type(pdg_array_t), intent(in) :: pa logical :: flag flag = pl%find_match (pa) /= 0 end function pdg_list_match_pdg_array function pdg_list_find_match_pdg_array (pl, pa, mask) result (i) class(pdg_list_t), intent(in) :: pl type(pdg_array_t), intent(in) :: pa logical, dimension(:), intent(in), optional :: mask integer :: i do i = 1, size (pl%a) if (present (mask)) then if (.not. mask(i)) cycle end if if (pl%a(i) .match. pa) return end do i = 0 end function pdg_list_find_match_pdg_array @ %def pdg_list_match_pdg_array @ %def pdg_list_find_match_pdg_array @ Some old compilers have problems with allocatable arrays as intent(out) or as function result, so be conservative here: <>= procedure :: create_pdg_array => pdg_list_create_pdg_array <>= subroutine pdg_list_create_pdg_array (pl, pdg) class(pdg_list_t), intent(in) :: pl type(pdg_array_t), dimension(:), intent(inout), allocatable :: pdg integer :: n_elements integer :: i associate (a => pl%a) n_elements = size (a) if (allocated (pdg)) deallocate (pdg) allocate (pdg (n_elements)) do i = 1, n_elements pdg(i) = a(i) end do end associate end subroutine pdg_list_create_pdg_array @ %def pdg_list_create_pdg_array @ <>= procedure :: create_antiparticles => pdg_list_create_antiparticles <>= subroutine pdg_list_create_antiparticles (pl, pl_anti, n_new_particles) class(pdg_list_t), intent(in) :: pl type(pdg_list_t), intent(out) :: pl_anti integer, intent(out) :: n_new_particles type(pdg_list_t) :: pl_inverse integer :: i, n integer :: n_identical logical, dimension(:), allocatable :: collect n = pl%get_size (); n_identical = 0 allocate (collect (n)); collect = .true. call pl_inverse%init (n) do i = 1, n pl_inverse%a(i) = pl%a(i)%invert() end do do i = 1, n if (any (pl_inverse%a(i) == pl%a)) then collect(i) = .false. n_identical = n_identical + 1 end if end do n_new_particles = n - n_identical if (n_new_particles > 0) then call pl_anti%init (n_new_particles) do i = 1, n if (collect (i)) pl_anti%a(i) = pl_inverse%a(i) end do end if end subroutine pdg_list_create_antiparticles @ %def pdg_list_create_antiparticles @ <>= procedure :: search_for_particle => pdg_list_search_for_particle <>= elemental function pdg_list_search_for_particle (pl, i_part) result (found) logical :: found class(pdg_list_t), intent(in) :: pl integer, intent(in) :: i_part integer :: i_pl do i_pl = 1, size (pl%a) found = pl%a(i_pl)%search_for_particle (i_part) if (found) return end do end function pdg_list_search_for_particle @ %def pdg_list_search_for_particle @ <>= procedure :: contains_colored_particles => pdg_list_contains_colored_particles <>= function pdg_list_contains_colored_particles (pl) result (colored) class(pdg_list_t), intent(in) :: pl logical :: colored integer :: i colored = .false. do i = 1, size (pl%a) if (pl%a(i)%has_colored_particles()) then colored = .true. exit end if end do end function pdg_list_contains_colored_particles @ %def pdg_list_contains_colored_particles @ \subsection{Unit tests} Test module, followed by the corresponding implementation module. <<[[pdg_arrays_ut.f90]]>>= <> module pdg_arrays_ut use unit_tests use pdg_arrays_uti <> <> contains <> end module pdg_arrays_ut @ %def pdg_arrays_ut @ <<[[pdg_arrays_uti.f90]]>>= <> module pdg_arrays_uti use pdg_arrays <> <> contains <> end module pdg_arrays_uti @ %def pdg_arrays_ut @ API: driver for the unit tests below. <>= public :: pdg_arrays_test <>= subroutine pdg_arrays_test (u, results) integer, intent(in) :: u type (test_results_t), intent(inout) :: results <> end subroutine pdg_arrays_test @ %def pdg_arrays_test @ Basic functionality. <>= call test (pdg_arrays_1, "pdg_arrays_1", & "create and sort PDG array", & u, results) <>= public :: pdg_arrays_1 <>= subroutine pdg_arrays_1 (u) integer, intent(in) :: u type(pdg_array_t) :: pa, pa1, pa2, pa3, pa4, pa5, pa6 integer, dimension(:), allocatable :: pdg write (u, "(A)") "* Test output: pdg_arrays_1" write (u, "(A)") "* Purpose: create and sort PDG arrays" write (u, "(A)") write (u, "(A)") "* Assignment" write (u, "(A)") call pa%write (u) write (u, *) write (u, "(A,I0)") "length = ", pa%get_length () pdg = pa write (u, "(A,3(1x,I0))") "contents = ", pdg write (u, *) pa = 1 call pa%write (u) write (u, *) write (u, "(A,I0)") "length = ", pa%get_length () pdg = pa write (u, "(A,3(1x,I0))") "contents = ", pdg write (u, *) pa = [1, 2, 3] call pa%write (u) write (u, *) write (u, "(A,I0)") "length = ", pa%get_length () pdg = pa write (u, "(A,3(1x,I0))") "contents = ", pdg write (u, "(A,I0)") "element #2 = ", pa%get (2) write (u, *) write (u, "(A)") "* Replace" write (u, *) pa = pa%replace (2, [-5, 5, -7]) call pa%write (u) write (u, *) write (u, *) write (u, "(A)") "* Sort" write (u, *) pa = [1, -7, 3, -5, 5, 3] call pa%write (u) write (u, *) pa1 = pa%sort_abs () pa2 = pa%sort_abs (unique = .true.) call pa1%write (u) write (u, *) call pa2%write (u) write (u, *) write (u, *) write (u, "(A)") "* Compare" write (u, *) pa1 = [1, 3] pa2 = [1, 2, -2] pa3 = [1, 2, 4] pa4 = [1, 2, 4] pa5 = [1, 2, -4] pa6 = [1, 2, -3] write (u, "(A,6(1x,L1))") "< ", & pa1 < pa2, pa2 < pa3, pa3 < pa4, pa4 < pa5, pa5 < pa6, pa6 < pa1 write (u, "(A,6(1x,L1))") "> ", & pa1 > pa2, pa2 > pa3, pa3 > pa4, pa4 > pa5, pa5 > pa6, pa6 > pa1 write (u, "(A,6(1x,L1))") "<=", & pa1 <= pa2, pa2 <= pa3, pa3 <= pa4, pa4 <= pa5, pa5 <= pa6, pa6 <= pa1 write (u, "(A,6(1x,L1))") ">=", & pa1 >= pa2, pa2 >= pa3, pa3 >= pa4, pa4 >= pa5, pa5 >= pa6, pa6 >= pa1 write (u, "(A,6(1x,L1))") "==", & pa1 == pa2, pa2 == pa3, pa3 == pa4, pa4 == pa5, pa5 == pa6, pa6 == pa1 write (u, "(A,6(1x,L1))") "/=", & pa1 /= pa2, pa2 /= pa3, pa3 /= pa4, pa4 /= pa5, pa5 /= pa6, pa6 /= pa1 write (u, *) pa1 = [0] pa2 = [1, 2] pa3 = [1, -2] write (u, "(A,6(1x,L1))") "eqv ", & pa1 .eqv. pa1, pa1 .eqv. pa2, & pa2 .eqv. pa2, pa2 .eqv. pa3 write (u, "(A,6(1x,L1))") "neqv", & pa1 .neqv. pa1, pa1 .neqv. pa2, & pa2 .neqv. pa2, pa2 .neqv. pa3 write (u, *) write (u, "(A,6(1x,L1))") "match", & pa1 .match. 0, pa1 .match. 1, & pa2 .match. 0, pa2 .match. 1, pa2 .match. 3 write (u, "(A)") write (u, "(A)") "* Test output end: pdg_arrays_1" end subroutine pdg_arrays_1 @ %def pdg_arrays_1 @ PDG array list, i.e., arrays of arrays. <>= call test (pdg_arrays_2, "pdg_arrays_2", & "create and sort PDG lists", & u, results) <>= public :: pdg_arrays_2 <>= subroutine pdg_arrays_2 (u) integer, intent(in) :: u type(pdg_array_t) :: pa type(pdg_list_t) :: pl, pl1 write (u, "(A)") "* Test output: pdg_arrays_2" write (u, "(A)") "* Purpose: create and sort PDG lists" write (u, "(A)") write (u, "(A)") "* Assignment" write (u, "(A)") call pl%init (3) call pl%set (1, 42) call pl%set (2, [3, 2]) pa = [5, -5] call pl%set (3, pa) call pl%write (u) write (u, *) write (u, "(A,I0)") "size = ", pl%get_size () write (u, "(A)") write (u, "(A)") "* Sort" write (u, "(A)") pl = pl%sort_abs () call pl%write (u) write (u, *) write (u, "(A)") write (u, "(A)") "* Extract item #3" write (u, "(A)") pa = pl%get (3) call pa%write (u) write (u, *) write (u, "(A)") write (u, "(A)") "* Replace item #3" write (u, "(A)") call pl1%init (2) call pl1%set (1, [2, 4]) call pl1%set (2, -7) pl = pl%replace (3, pl1) call pl%write (u) write (u, *) write (u, "(A)") write (u, "(A)") "* Test output end: pdg_arrays_2" end subroutine pdg_arrays_2 @ %def pdg_arrays_2 @ Check if a (sorted) PDG array lists is regular. The entries (PDG arrays) must not overlap, unless they are identical. <>= call test (pdg_arrays_3, "pdg_arrays_3", & "check PDG lists", & u, results) <>= public :: pdg_arrays_3 <>= subroutine pdg_arrays_3 (u) integer, intent(in) :: u type(pdg_list_t) :: pl write (u, "(A)") "* Test output: pdg_arrays_3" write (u, "(A)") "* Purpose: check for regular PDG lists" write (u, "(A)") write (u, "(A)") "* Regular list" write (u, "(A)") call pl%init (4) call pl%set (1, [1, 2]) call pl%set (2, [1, 2]) call pl%set (3, [5, -5]) call pl%set (4, 42) call pl%write (u) write (u, *) write (u, "(L1)") pl%is_regular () write (u, "(A)") write (u, "(A)") "* Irregular list" write (u, "(A)") call pl%init (4) call pl%set (1, [1, 2]) call pl%set (2, [1, 2]) call pl%set (3, [2, 5, -5]) call pl%set (4, 42) call pl%write (u) write (u, *) write (u, "(L1)") pl%is_regular () write (u, "(A)") write (u, "(A)") "* Test output end: pdg_arrays_3" end subroutine pdg_arrays_3 @ %def pdg_arrays_3 @ Compare PDG array lists. The lists must be regular, i.e., sorted and with non-overlapping (or identical) entries. <>= call test (pdg_arrays_4, "pdg_arrays_4", & "compare PDG lists", & u, results) <>= public :: pdg_arrays_4 <>= subroutine pdg_arrays_4 (u) integer, intent(in) :: u type(pdg_list_t) :: pl1, pl2, pl3 write (u, "(A)") "* Test output: pdg_arrays_4" write (u, "(A)") "* Purpose: check for regular PDG lists" write (u, "(A)") write (u, "(A)") "* Create lists" write (u, "(A)") call pl1%init (4) call pl1%set (1, [1, 2]) call pl1%set (2, [1, 2]) call pl1%set (3, [5, -5]) call pl1%set (4, 42) write (u, "(I1,1x)", advance = "no") 1 call pl1%write (u) write (u, *) call pl2%init (2) call pl2%set (1, 3) call pl2%set (2, [5, -5]) write (u, "(I1,1x)", advance = "no") 2 call pl2%write (u) write (u, *) call pl3%init (2) call pl3%set (1, 4) call pl3%set (2, [5, -5]) write (u, "(I1,1x)", advance = "no") 3 call pl3%write (u) write (u, *) write (u, "(A)") write (u, "(A)") "* a == b" write (u, "(A)") write (u, "(2x,A)") "123" write (u, *) write (u, "(I1,1x,4L1)") 1, pl1 == pl1, pl1 == pl2, pl1 == pl3 write (u, "(I1,1x,4L1)") 2, pl2 == pl1, pl2 == pl2, pl2 == pl3 write (u, "(I1,1x,4L1)") 3, pl3 == pl1, pl3 == pl2, pl3 == pl3 write (u, "(A)") write (u, "(A)") "* a < b" write (u, "(A)") write (u, "(2x,A)") "123" write (u, *) write (u, "(I1,1x,4L1)") 1, pl1 < pl1, pl1 < pl2, pl1 < pl3 write (u, "(I1,1x,4L1)") 2, pl2 < pl1, pl2 < pl2, pl2 < pl3 write (u, "(I1,1x,4L1)") 3, pl3 < pl1, pl3 < pl2, pl3 < pl3 write (u, "(A)") write (u, "(A)") "* Test output end: pdg_arrays_4" end subroutine pdg_arrays_4 @ %def pdg_arrays_4 @ Match-replace: translate all entries in the first list into the matching entries of the second list, if there is a match. <>= call test (pdg_arrays_5, "pdg_arrays_5", & "match PDG lists", & u, results) <>= public :: pdg_arrays_5 <>= subroutine pdg_arrays_5 (u) integer, intent(in) :: u type(pdg_list_t) :: pl1, pl2, pl3 logical :: success write (u, "(A)") "* Test output: pdg_arrays_5" write (u, "(A)") "* Purpose: match-replace" write (u, "(A)") write (u, "(A)") "* Create lists" write (u, "(A)") call pl1%init (2) call pl1%set (1, [1, 2]) call pl1%set (2, 42) call pl1%write (u) write (u, *) call pl3%init (2) call pl3%set (1, [42, -42]) call pl3%set (2, [1, 2, 3, 4]) call pl1%match_replace (pl3, success) call pl3%write (u) write (u, "(1x,A,1x,L1,':',1x)", advance="no") "=>", success call pl1%write (u) write (u, *) write (u, *) call pl2%init (2) call pl2%set (1, 9) call pl2%set (2, 42) call pl2%write (u) write (u, *) call pl2%match_replace (pl3, success) call pl3%write (u) write (u, "(1x,A,1x,L1,':',1x)", advance="no") "=>", success call pl2%write (u) write (u, *) write (u, "(A)") write (u, "(A)") "* Test output end: pdg_arrays_5" end subroutine pdg_arrays_5 @ %def pdg_arrays_5 @ \clearpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Jets} The FastJet library is linked externally, if available. The wrapper code is also in a separate directory. Here, we define \whizard-specific procedures and tests. <<[[jets.f90]]>>= <> module jets use fastjet !NODEP! <> <> contains <> end module jets @ %def jets @ \subsection{Re-exported symbols} We use this module as a proxy for the FastJet interface, therefore we re-export some symbols. <>= public :: fastjet_available public :: fastjet_init public :: jet_definition_t public :: pseudojet_t public :: pseudojet_vector_t public :: cluster_sequence_t public :: assignment (=) @ %def jet_definition_t pseudojet_t pseudojet_vector_t cluster_sequence_t @ The initialization routine prints the banner. <>= subroutine fastjet_init () call print_banner () end subroutine fastjet_init @ %def fastjet_init @ The jet algorithm codes (actually, integers) <>= public :: kt_algorithm public :: cambridge_algorithm public :: antikt_algorithm public :: genkt_algorithm public :: cambridge_for_passive_algorithm public :: genkt_for_passive_algorithm public :: ee_kt_algorithm public :: ee_genkt_algorithm public :: plugin_algorithm public :: undefined_jet_algorithm @ \subsection{Unit tests} Test module, followed by the corresponding implementation module. <<[[jets_ut.f90]]>>= <> module jets_ut use unit_tests use jets_uti <> <> contains <> end module jets_ut @ %def jets_ut @ <<[[jets_uti.f90]]>>= <> module jets_uti <> use fastjet !NODEP! use jets <> <> contains <> end module jets_uti @ %def jets_ut @ API: driver for the unit tests below. <>= public :: jets_test <>= subroutine jets_test (u, results) integer, intent(in) :: u type (test_results_t), intent(inout) :: results <> end subroutine jets_test @ %def jets_test @ This test is actually the minimal example from the FastJet manual, translated to Fortran. Note that FastJet creates pseudojet vectors, which we mirror in the [[pseudojet_vector_t]], but immediately assign to pseudojet arrays. Without automatic finalization available in the compilers, we should avoid this in actual code and rather introduce intermediate variables for those objects, which we can finalize explicitly. <>= call test (jets_1, "jets_1", & "basic FastJet functionality", & u, results) <>= public :: jets_1 <>= subroutine jets_1 (u) integer, intent(in) :: u type(pseudojet_t), dimension(:), allocatable :: prt, jets, constituents type(jet_definition_t) :: jet_def type(cluster_sequence_t) :: cs integer, parameter :: dp = default integer :: i, j write (u, "(A)") "* Test output: jets_1" write (u, "(A)") "* Purpose: test basic FastJet functionality" write (u, "(A)") write (u, "(A)") "* Print banner" call print_banner () write (u, *) write (u, "(A)") "* Prepare input particles" allocate (prt (3)) call prt(1)%init ( 99._dp, 0.1_dp, 0._dp, 100._dp) call prt(2)%init ( 4._dp,-0.1_dp, 0._dp, 5._dp) call prt(3)%init (-99._dp, 0._dp, 0._dp, 99._dp) write (u, *) write (u, "(A)") "* Define jet algorithm" call jet_def%init (antikt_algorithm, 0.7_dp) write (u, *) write (u, "(A)") "* Cluster particles according to jet algorithm" write (u, *) write (u, "(A,A)") "Clustering with ", jet_def%description () call cs%init (pseudojet_vector (prt), jet_def) write (u, *) write (u, "(A)") "* Sort output jets" jets = sorted_by_pt (cs%inclusive_jets ()) write (u, *) write (u, "(A)") "* Print jet observables and constituents" write (u, *) write (u, "(4x,3(7x,A3))") "pt", "y", "phi" do i = 1, size (jets) write (u, "(A,1x,I0,A,3(1x,F9.5))") & "jet", i, ":", jets(i)%perp (), jets(i)%rap (), jets(i)%phi () constituents = jets(i)%constituents () do j = 1, size (constituents) write (u, "(4x,A,1x,I0,A,F9.5)") & "constituent", j, "'s pt:", constituents(j)%perp () end do do j = 1, size (constituents) call constituents(j)%final () end do end do write (u, *) write (u, "(A)") "* Cleanup" do i = 1, size (prt) call prt(i)%final () end do do i = 1, size (jets) call jets(i)%final () end do call jet_def%final () call cs%final () write (u, "(A)") write (u, "(A)") "* Test output end: jets_1" end subroutine jets_1 @ %def jets_1 @ \clearpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Subevents} The purpose of subevents is to store the relevant part of the physical event (either partonic or hadronic), and to hold particle selections and combinations which are constructed in cut or analysis expressions. <<[[subevents.f90]]>>= <> module subevents use, intrinsic :: iso_c_binding !NODEP! <> use io_units use format_defs, only: FMT_14, FMT_19 use format_utils, only: pac_fmt use physics_defs use sorting use c_particles use lorentz use pdg_arrays use jets <> <> <> <> <> contains <> end module subevents @ %def subevents @ \subsection{Particles} For the purpose of this module, a particle has a type which can indicate a beam, incoming, outgoing, or composite particle, flavor and helicity codes (integer, undefined for composite), four-momentum and invariant mass squared. (Other particles types are used in extended event types, but also defined here.) Furthermore, each particle has an allocatable array of ancestors -- particle indices which indicate the building blocks of a composite particle. For an incoming/outgoing particle, the array contains only the index of the particle itself. For incoming particles, the momentum is inverted before storing it in the particle object. <>= integer, parameter, public :: PRT_UNDEFINED = 0 integer, parameter, public :: PRT_BEAM = -9 integer, parameter, public :: PRT_INCOMING = 1 integer, parameter, public :: PRT_OUTGOING = 2 integer, parameter, public :: PRT_COMPOSITE = 3 integer, parameter, public :: PRT_VIRTUAL = 4 integer, parameter, public :: PRT_RESONANT = 5 integer, parameter, public :: PRT_BEAM_REMNANT = 9 @ %def PRT_UNDEFINED PRT_BEAM @ %def PRT_INCOMING PRT_OUTGOING PRT_COMPOSITE @ %def PRT_COMPOSITE PRT_VIRTUAL PRT_RESONANT @ %def PRT_BEAM_REMNANT @ \subsubsection{The type} We initialize only the type here and mark as unpolarized. The initializers below do the rest. The logicals [[is_b_jet]] and [[is_c_jet]] are true only if [[prt_t]] comes out of the [[subevt_cluster]] routine and fulfils the correct flavor content. <>= public :: prt_t <>= type :: prt_t private integer :: type = PRT_UNDEFINED integer :: pdg logical :: polarized = .false. logical :: colorized = .false. logical :: clustered = .false. logical :: is_b_jet = .false. logical :: is_c_jet = .false. integer :: h type(vector4_t) :: p real(default) :: p2 integer, dimension(:), allocatable :: src integer, dimension(:), allocatable :: col integer, dimension(:), allocatable :: acl end type prt_t @ %def prt_t @ Initializers. Polarization is set separately. Finalizers are not needed. <>= subroutine prt_init_beam (prt, pdg, p, p2, src) type(prt_t), intent(out) :: prt integer, intent(in) :: pdg type(vector4_t), intent(in) :: p real(default), intent(in) :: p2 integer, dimension(:), intent(in) :: src prt%type = PRT_BEAM call prt_set (prt, pdg, - p, p2, src) end subroutine prt_init_beam subroutine prt_init_incoming (prt, pdg, p, p2, src) type(prt_t), intent(out) :: prt integer, intent(in) :: pdg type(vector4_t), intent(in) :: p real(default), intent(in) :: p2 integer, dimension(:), intent(in) :: src prt%type = PRT_INCOMING call prt_set (prt, pdg, - p, p2, src) end subroutine prt_init_incoming subroutine prt_init_outgoing (prt, pdg, p, p2, src) type(prt_t), intent(out) :: prt integer, intent(in) :: pdg type(vector4_t), intent(in) :: p real(default), intent(in) :: p2 integer, dimension(:), intent(in) :: src prt%type = PRT_OUTGOING call prt_set (prt, pdg, p, p2, src) end subroutine prt_init_outgoing subroutine prt_init_composite (prt, p, src) type(prt_t), intent(out) :: prt type(vector4_t), intent(in) :: p integer, dimension(:), intent(in) :: src prt%type = PRT_COMPOSITE call prt_set (prt, 0, p, p**2, src) end subroutine prt_init_composite @ %def prt_init_beam prt_init_incoming prt_init_outgoing prt_init_composite @ This version is for temporary particle objects, so the [[src]] array is not set. <>= public :: prt_init_combine <>= subroutine prt_init_combine (prt, prt1, prt2) type(prt_t), intent(out) :: prt type(prt_t), intent(in) :: prt1, prt2 type(vector4_t) :: p integer, dimension(0) :: src prt%type = PRT_COMPOSITE p = prt1%p + prt2%p call prt_set (prt, 0, p, p**2, src) end subroutine prt_init_combine @ %def prt_init_combine @ Init from a pseudojet object. <>= subroutine prt_init_pseudojet (prt, jet, src, pdg, is_b_jet, is_c_jet) type(prt_t), intent(out) :: prt type(pseudojet_t), intent(in) :: jet integer, dimension(:), intent(in) :: src integer, intent(in) :: pdg logical, intent(in) :: is_b_jet, is_c_jet type(vector4_t) :: p prt%type = PRT_COMPOSITE p = vector4_moving (jet%e(), & vector3_moving ([jet%px(), jet%py(), jet%pz()])) call prt_set (prt, pdg, p, p**2, src) prt%is_b_jet = is_b_jet prt%is_c_jet = is_c_jet prt%clustered = .true. end subroutine prt_init_pseudojet @ %def prt_init_pseudojet @ \subsubsection{Accessing contents} <>= public :: prt_get_pdg <>= elemental function prt_get_pdg (prt) result (pdg) integer :: pdg type(prt_t), intent(in) :: prt pdg = prt%pdg end function prt_get_pdg @ %def prt_get_pdg <>= public :: prt_get_momentum <>= elemental function prt_get_momentum (prt) result (p) type(vector4_t) :: p type(prt_t), intent(in) :: prt p = prt%p end function prt_get_momentum @ %def prt_get_momentum <>= public :: prt_get_msq <>= elemental function prt_get_msq (prt) result (msq) real(default) :: msq type(prt_t), intent(in) :: prt msq = prt%p2 end function prt_get_msq @ %def prt_get_msq <>= public :: prt_is_polarized <>= elemental function prt_is_polarized (prt) result (flag) logical :: flag type(prt_t), intent(in) :: prt flag = prt%polarized end function prt_is_polarized @ %def prt_is_polarized <>= public :: prt_get_helicity <>= elemental function prt_get_helicity (prt) result (h) integer :: h type(prt_t), intent(in) :: prt h = prt%h end function prt_get_helicity @ %def prt_get_helicity <>= public :: prt_is_colorized <>= elemental function prt_is_colorized (prt) result (flag) logical :: flag type(prt_t), intent(in) :: prt flag = prt%colorized end function prt_is_colorized @ %def prt_is_colorized <>= public :: prt_is_clustered <>= elemental function prt_is_clustered (prt) result (flag) logical :: flag type(prt_t), intent(in) :: prt flag = prt%clustered end function prt_is_clustered @ %def prt_is_clustered <>= public :: prt_is_recombinable <>= elemental function prt_is_recombinable (prt) result (flag) logical :: flag type(prt_t), intent(in) :: prt flag = prt_is_parton (prt) .or. & abs(prt%pdg) == TOP_Q .or. & prt_is_lepton (prt) .or. & prt_is_photon (prt) end function prt_is_recombinable @ %def prt_is_recombinable <>= public :: prt_is_photon <>= elemental function prt_is_photon (prt) result (flag) logical :: flag type(prt_t), intent(in) :: prt flag = prt%pdg == PHOTON end function prt_is_photon @ %def prt_is_photon We do not take the top quark into account here. <>= public :: prt_is_parton <>= elemental function prt_is_parton (prt) result (flag) logical :: flag type(prt_t), intent(in) :: prt flag = abs(prt%pdg) == DOWN_Q .or. & abs(prt%pdg) == UP_Q .or. & abs(prt%pdg) == STRANGE_Q .or. & abs(prt%pdg) == CHARM_Q .or. & abs(prt%pdg) == BOTTOM_Q .or. & prt%pdg == GLUON end function prt_is_parton @ %def prt_is_parton <>= public :: prt_is_lepton <>= elemental function prt_is_lepton (prt) result (flag) logical :: flag type(prt_t), intent(in) :: prt flag = abs(prt%pdg) == ELECTRON .or. & abs(prt%pdg) == MUON .or. & abs(prt%pdg) == TAU end function prt_is_lepton @ %def prt_is_lepton <>= public :: prt_is_b_jet <>= elemental function prt_is_b_jet (prt) result (flag) logical :: flag type(prt_t), intent(in) :: prt flag = prt%is_b_jet end function prt_is_b_jet @ %def prt_is_b_jet <>= public :: prt_is_c_jet <>= elemental function prt_is_c_jet (prt) result (flag) logical :: flag type(prt_t), intent(in) :: prt flag = prt%is_c_jet end function prt_is_c_jet @ %def prt_is_c_jet @ The number of open color (anticolor) lines. We inspect the list of color (anticolor) lines and count the entries that do not appear in the list of anticolors (colors). (There is no check against duplicates; we assume that color line indices are unique.) <>= public :: prt_get_n_col public :: prt_get_n_acl <>= elemental function prt_get_n_col (prt) result (n) integer :: n type(prt_t), intent(in) :: prt integer, dimension(:), allocatable :: col, acl integer :: i n = 0 if (prt%colorized) then do i = 1, size (prt%col) if (all (prt%col(i) /= prt%acl)) n = n + 1 end do end if end function prt_get_n_col elemental function prt_get_n_acl (prt) result (n) integer :: n type(prt_t), intent(in) :: prt integer, dimension(:), allocatable :: col, acl integer :: i n = 0 if (prt%colorized) then do i = 1, size (prt%acl) if (all (prt%acl(i) /= prt%col)) n = n + 1 end do end if end function prt_get_n_acl @ %def prt_get_n_col @ %def prt_get_n_acl @ Return the color and anticolor-flow line indices explicitly. <>= public :: prt_get_color_indices <>= subroutine prt_get_color_indices (prt, col, acl) type(prt_t), intent(in) :: prt integer, dimension(:), allocatable, intent(out) :: col, acl if (prt%colorized) then col = prt%col acl = prt%acl else col = [integer::] acl = [integer::] end if end subroutine prt_get_color_indices @ %def prt_get_color_indices @ \subsubsection{Setting data} Set the PDG, momentum and momentum squared, and ancestors. If allocate-on-assignment is available, this can be simplified. <>= subroutine prt_set (prt, pdg, p, p2, src) type(prt_t), intent(inout) :: prt integer, intent(in) :: pdg type(vector4_t), intent(in) :: p real(default), intent(in) :: p2 integer, dimension(:), intent(in) :: src prt%pdg = pdg prt%p = p prt%p2 = p2 if (allocated (prt%src)) then if (size (src) /= size (prt%src)) then deallocate (prt%src) allocate (prt%src (size (src))) end if else allocate (prt%src (size (src))) end if prt%src = src end subroutine prt_set @ %def prt_set @ Set the particle PDG code separately. <>= elemental subroutine prt_set_pdg (prt, pdg) type(prt_t), intent(inout) :: prt integer, intent(in) :: pdg prt%pdg = pdg end subroutine prt_set_pdg @ %def prt_set_pdg @ Set the momentum separately. <>= elemental subroutine prt_set_p (prt, p) type(prt_t), intent(inout) :: prt type(vector4_t), intent(in) :: p prt%p = p end subroutine prt_set_p @ %def prt_set_p @ Set the squared invariant mass separately. <>= elemental subroutine prt_set_p2 (prt, p2) type(prt_t), intent(inout) :: prt real(default), intent(in) :: p2 prt%p2 = p2 end subroutine prt_set_p2 @ %def prt_set_p2 @ Set helicity (optional). <>= subroutine prt_polarize (prt, h) type(prt_t), intent(inout) :: prt integer, intent(in) :: h prt%polarized = .true. prt%h = h end subroutine prt_polarize @ %def prt_polarize @ Set color-flow indices (optional). <>= subroutine prt_colorize (prt, col, acl) type(prt_t), intent(inout) :: prt integer, dimension(:), intent(in) :: col, acl prt%colorized = .true. prt%col = col prt%acl = acl end subroutine prt_colorize @ %def prt_colorize @ \subsubsection{Conversion} Transform a [[prt_t]] object into a [[c_prt_t]] object. <>= public :: c_prt <>= interface c_prt module procedure c_prt_from_prt end interface @ %def c_prt <>= elemental function c_prt_from_prt (prt) result (c_prt) type(c_prt_t) :: c_prt type(prt_t), intent(in) :: prt c_prt = prt%p c_prt%type = prt%type c_prt%pdg = prt%pdg if (prt%polarized) then c_prt%polarized = 1 else c_prt%polarized = 0 end if c_prt%h = prt%h end function c_prt_from_prt @ %def c_prt_from_prt @ \subsubsection{Output} <>= public :: prt_write <>= subroutine prt_write (prt, unit, testflag) type(prt_t), intent(in) :: prt integer, intent(in), optional :: unit logical, intent(in), optional :: testflag logical :: pacified type(prt_t) :: tmp character(len=7) :: fmt integer :: u, i call pac_fmt (fmt, FMT_19, FMT_14, testflag) u = given_output_unit (unit); if (u < 0) return pacified = .false. ; if (present (testflag)) pacified = testflag tmp = prt if (pacified) call pacify (tmp) write (u, "(1x,A)", advance="no") "prt(" select case (prt%type) case (PRT_UNDEFINED); write (u, "('?')", advance="no") case (PRT_BEAM); write (u, "('b:')", advance="no") case (PRT_INCOMING); write (u, "('i:')", advance="no") case (PRT_OUTGOING); write (u, "('o:')", advance="no") case (PRT_COMPOSITE); write (u, "('c:')", advance="no") end select select case (prt%type) case (PRT_BEAM, PRT_INCOMING, PRT_OUTGOING) if (prt%polarized) then write (u, "(I0,'/',I0,'|')", advance="no") prt%pdg, prt%h else write (u, "(I0,'|')", advance="no") prt%pdg end if end select select case (prt%type) case (PRT_BEAM, PRT_INCOMING, PRT_OUTGOING, PRT_COMPOSITE) if (prt%colorized) then write (u, "(*(I0,:,','))", advance="no") prt%col write (u, "('/')", advance="no") write (u, "(*(I0,:,','))", advance="no") prt%acl write (u, "('|')", advance="no") end if end select select case (prt%type) case (PRT_BEAM, PRT_INCOMING, PRT_OUTGOING, PRT_COMPOSITE) write (u, "(" // FMT_14 // ",';'," // FMT_14 // ",','," // & FMT_14 // ",','," // FMT_14 // ")", advance="no") tmp%p write (u, "('|'," // fmt // ")", advance="no") tmp%p2 end select if (allocated (prt%src)) then write (u, "('|')", advance="no") do i = 1, size (prt%src) write (u, "(1x,I0)", advance="no") prt%src(i) end do end if if (prt%is_b_jet) then write (u, "('|b jet')", advance="no") end if if (prt%is_c_jet) then write (u, "('|c jet')", advance="no") end if write (u, "(A)") ")" end subroutine prt_write @ %def prt_write @ \subsubsection{Tools} Two particles match if their [[src]] arrays are the same. <>= public :: operator(.match.) <>= interface operator(.match.) module procedure prt_match end interface @ %def .match. <>= elemental function prt_match (prt1, prt2) result (match) logical :: match type(prt_t), intent(in) :: prt1, prt2 if (size (prt1%src) == size (prt2%src)) then match = all (prt1%src == prt2%src) else match = .false. end if end function prt_match @ %def prt_match @ The combine operation makes a pseudoparticle whose momentum is the result of adding (the momenta of) the pair of input particles. We trace the particles from which a particle is built by storing a [[src]] array. Each particle entry in the [[src]] list contains a list of indices which indicates its building blocks. The indices refer to an original list of particles. Index lists are sorted, and they contain no element more than once. We thus require that in a given pseudoparticle, each original particle occurs at most once. The result is intent(inout), so it will not be initialized when the subroutine is entered. If the particles carry color, we recall that the combined particle is a composite which is understood as outgoing. If one of the arguments is an incoming particle, is color entries must be reversed. <>= subroutine prt_combine (prt, prt_in1, prt_in2, ok) type(prt_t), intent(inout) :: prt type(prt_t), intent(in) :: prt_in1, prt_in2 logical :: ok integer, dimension(:), allocatable :: src integer, dimension(:), allocatable :: col1, acl1, col2, acl2 call combine_index_lists (src, prt_in1%src, prt_in2%src) ok = allocated (src) if (ok) then call prt_init_composite (prt, prt_in1%p + prt_in2%p, src) if (prt_in1%colorized .or. prt_in2%colorized) then select case (prt_in1%type) case default call prt_get_color_indices (prt_in1, col1, acl1) case (PRT_BEAM, PRT_INCOMING) call prt_get_color_indices (prt_in1, acl1, col1) end select select case (prt_in2%type) case default call prt_get_color_indices (prt_in2, col2, acl2) case (PRT_BEAM, PRT_INCOMING) call prt_get_color_indices (prt_in2, acl2, col2) end select call prt_colorize (prt, [col1, col2], [acl1, acl2]) end if end if end subroutine prt_combine @ %def prt_combine @ This variant does not produce the combined particle, it just checks whether the combination is valid (no common [[src]] entry). <>= public :: are_disjoint <>= function are_disjoint (prt_in1, prt_in2) result (flag) logical :: flag type(prt_t), intent(in) :: prt_in1, prt_in2 flag = index_lists_are_disjoint (prt_in1%src, prt_in2%src) end function are_disjoint @ %def are_disjoint @ [[src]] Lists with length $>1$ are built by a [[combine]] operation which merges the lists in a sorted manner. If the result would have a duplicate entry, it is discarded, and the result is unallocated. <>= subroutine combine_index_lists (res, src1, src2) integer, dimension(:), intent(in) :: src1, src2 integer, dimension(:), allocatable :: res integer :: i1, i2, i allocate (res (size (src1) + size (src2))) if (size (src1) == 0) then res = src2 return else if (size (src2) == 0) then res = src1 return end if i1 = 1 i2 = 1 LOOP: do i = 1, size (res) if (src1(i1) < src2(i2)) then res(i) = src1(i1); i1 = i1 + 1 if (i1 > size (src1)) then res(i+1:) = src2(i2:) exit LOOP end if else if (src1(i1) > src2(i2)) then res(i) = src2(i2); i2 = i2 + 1 if (i2 > size (src2)) then res(i+1:) = src1(i1:) exit LOOP end if else deallocate (res) exit LOOP end if end do LOOP end subroutine combine_index_lists @ %def combine_index_lists @ This function is similar, but it does not actually merge the list, it just checks whether they are disjoint (no common [[src]] entry). <>= function index_lists_are_disjoint (src1, src2) result (flag) logical :: flag integer, dimension(:), intent(in) :: src1, src2 integer :: i1, i2, i flag = .true. i1 = 1 i2 = 1 LOOP: do i = 1, size (src1) + size (src2) if (src1(i1) < src2(i2)) then i1 = i1 + 1 if (i1 > size (src1)) then exit LOOP end if else if (src1(i1) > src2(i2)) then i2 = i2 + 1 if (i2 > size (src2)) then exit LOOP end if else flag = .false. exit LOOP end if end do LOOP end function index_lists_are_disjoint @ %def index_lists_are_disjoint @ \subsection{subevents} Particles are collected in subevents. This type is implemented as a dynamically allocated array, which need not be completely filled. The value [[n_active]] determines the number of meaningful entries. \subsubsection{Type definition} <>= public :: subevt_t <>= type :: subevt_t private integer :: n_active = 0 type(prt_t), dimension(:), allocatable :: prt contains <> end type subevt_t @ %def subevt_t @ Initialize, allocating with size zero (default) or given size. The number of contained particles is set equal to the size. <>= public :: subevt_init <>= subroutine subevt_init (subevt, n_active) type(subevt_t), intent(out) :: subevt integer, intent(in), optional :: n_active if (present (n_active)) subevt%n_active = n_active allocate (subevt%prt (subevt%n_active)) end subroutine subevt_init @ %def subevt_init @ (Re-)allocate the subevent with some given size. If the size is greater than the previous one, do a real reallocation. Otherwise, just reset the recorded size. Contents are untouched, but become invalid. <>= public :: subevt_reset <>= subroutine subevt_reset (subevt, n_active) type(subevt_t), intent(inout) :: subevt integer, intent(in) :: n_active subevt%n_active = n_active if (subevt%n_active > size (subevt%prt)) then deallocate (subevt%prt) allocate (subevt%prt (subevt%n_active)) end if end subroutine subevt_reset @ %def subevt_reset @ Output. No prefix for the headline 'subevt', because this will usually be printed appending to a previous line. <>= public :: subevt_write <>= procedure :: write => subevt_write <>= subroutine subevt_write (object, unit, prefix, pacified) class(subevt_t), intent(in) :: object integer, intent(in), optional :: unit character(*), intent(in), optional :: prefix logical, intent(in), optional :: pacified integer :: u, i u = given_output_unit (unit); if (u < 0) return write (u, "(1x,A)") "subevent:" do i = 1, object%n_active if (present (prefix)) write (u, "(A)", advance="no") prefix write (u, "(1x,I0)", advance="no") i call prt_write (object%prt(i), unit = unit, testflag = pacified) end do end subroutine subevt_write @ %def subevt_write @ Defined assignment: transfer only meaningful entries. This is a deep copy (as would be default assignment). <>= interface assignment(=) module procedure subevt_assign end interface @ %def = <>= subroutine subevt_assign (subevt, subevt_in) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: subevt_in if (.not. allocated (subevt%prt)) then call subevt_init (subevt, subevt_in%n_active) else call subevt_reset (subevt, subevt_in%n_active) end if subevt%prt(:subevt%n_active) = subevt_in%prt(:subevt%n_active) end subroutine subevt_assign @ %def subevt_assign @ \subsubsection{Fill contents} Store incoming/outgoing particles which are completely defined. <>= public :: subevt_set_beam public :: subevt_set_incoming public :: subevt_set_outgoing public :: subevt_set_composite <>= subroutine subevt_set_beam (subevt, i, pdg, p, p2, src) type(subevt_t), intent(inout) :: subevt integer, intent(in) :: i integer, intent(in) :: pdg type(vector4_t), intent(in) :: p real(default), intent(in) :: p2 integer, dimension(:), intent(in), optional :: src if (present (src)) then call prt_init_beam (subevt%prt(i), pdg, p, p2, src) else call prt_init_beam (subevt%prt(i), pdg, p, p2, [i]) end if end subroutine subevt_set_beam subroutine subevt_set_incoming (subevt, i, pdg, p, p2, src) type(subevt_t), intent(inout) :: subevt integer, intent(in) :: i integer, intent(in) :: pdg type(vector4_t), intent(in) :: p real(default), intent(in) :: p2 integer, dimension(:), intent(in), optional :: src if (present (src)) then call prt_init_incoming (subevt%prt(i), pdg, p, p2, src) else call prt_init_incoming (subevt%prt(i), pdg, p, p2, [i]) end if end subroutine subevt_set_incoming subroutine subevt_set_outgoing (subevt, i, pdg, p, p2, src) type(subevt_t), intent(inout) :: subevt integer, intent(in) :: i integer, intent(in) :: pdg type(vector4_t), intent(in) :: p real(default), intent(in) :: p2 integer, dimension(:), intent(in), optional :: src if (present (src)) then call prt_init_outgoing (subevt%prt(i), pdg, p, p2, src) else call prt_init_outgoing (subevt%prt(i), pdg, p, p2, [i]) end if end subroutine subevt_set_outgoing subroutine subevt_set_composite (subevt, i, p, src) type(subevt_t), intent(inout) :: subevt integer, intent(in) :: i type(vector4_t), intent(in) :: p integer, dimension(:), intent(in) :: src call prt_init_composite (subevt%prt(i), p, src) end subroutine subevt_set_composite @ %def subevt_set_incoming subevt_set_outgoing subevt_set_composite @ Separately assign flavors, simultaneously for all incoming/outgoing particles. <>= public :: subevt_set_pdg_beam public :: subevt_set_pdg_incoming public :: subevt_set_pdg_outgoing <>= subroutine subevt_set_pdg_beam (subevt, pdg) type(subevt_t), intent(inout) :: subevt integer, dimension(:), intent(in) :: pdg integer :: i, j j = 1 do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_BEAM) then call prt_set_pdg (subevt%prt(i), pdg(j)) j = j + 1 if (j > size (pdg)) exit end if end do end subroutine subevt_set_pdg_beam subroutine subevt_set_pdg_incoming (subevt, pdg) type(subevt_t), intent(inout) :: subevt integer, dimension(:), intent(in) :: pdg integer :: i, j j = 1 do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_INCOMING) then call prt_set_pdg (subevt%prt(i), pdg(j)) j = j + 1 if (j > size (pdg)) exit end if end do end subroutine subevt_set_pdg_incoming subroutine subevt_set_pdg_outgoing (subevt, pdg) type(subevt_t), intent(inout) :: subevt integer, dimension(:), intent(in) :: pdg integer :: i, j j = 1 do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_OUTGOING) then call prt_set_pdg (subevt%prt(i), pdg(j)) j = j + 1 if (j > size (pdg)) exit end if end do end subroutine subevt_set_pdg_outgoing @ %def subevt_set_pdg_beam @ %def subevt_set_pdg_incoming @ %def subevt_set_pdg_outgoing @ Separately assign momenta, simultaneously for all incoming/outgoing particles. <>= public :: subevt_set_p_beam public :: subevt_set_p_incoming public :: subevt_set_p_outgoing <>= subroutine subevt_set_p_beam (subevt, p) type(subevt_t), intent(inout) :: subevt type(vector4_t), dimension(:), intent(in) :: p integer :: i, j j = 1 do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_BEAM) then call prt_set_p (subevt%prt(i), p(j)) j = j + 1 if (j > size (p)) exit end if end do end subroutine subevt_set_p_beam subroutine subevt_set_p_incoming (subevt, p) type(subevt_t), intent(inout) :: subevt type(vector4_t), dimension(:), intent(in) :: p integer :: i, j j = 1 do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_INCOMING) then call prt_set_p (subevt%prt(i), p(j)) j = j + 1 if (j > size (p)) exit end if end do end subroutine subevt_set_p_incoming subroutine subevt_set_p_outgoing (subevt, p) type(subevt_t), intent(inout) :: subevt type(vector4_t), dimension(:), intent(in) :: p integer :: i, j j = 1 do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_OUTGOING) then call prt_set_p (subevt%prt(i), p(j)) j = j + 1 if (j > size (p)) exit end if end do end subroutine subevt_set_p_outgoing @ %def subevt_set_p_beam @ %def subevt_set_p_incoming @ %def subevt_set_p_outgoing @ Separately assign the squared invariant mass, simultaneously for all incoming/outgoing particles. <>= public :: subevt_set_p2_beam public :: subevt_set_p2_incoming public :: subevt_set_p2_outgoing <>= subroutine subevt_set_p2_beam (subevt, p2) type(subevt_t), intent(inout) :: subevt real(default), dimension(:), intent(in) :: p2 integer :: i, j j = 1 do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_BEAM) then call prt_set_p2 (subevt%prt(i), p2(j)) j = j + 1 if (j > size (p2)) exit end if end do end subroutine subevt_set_p2_beam subroutine subevt_set_p2_incoming (subevt, p2) type(subevt_t), intent(inout) :: subevt real(default), dimension(:), intent(in) :: p2 integer :: i, j j = 1 do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_INCOMING) then call prt_set_p2 (subevt%prt(i), p2(j)) j = j + 1 if (j > size (p2)) exit end if end do end subroutine subevt_set_p2_incoming subroutine subevt_set_p2_outgoing (subevt, p2) type(subevt_t), intent(inout) :: subevt real(default), dimension(:), intent(in) :: p2 integer :: i, j j = 1 do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_OUTGOING) then call prt_set_p2 (subevt%prt(i), p2(j)) j = j + 1 if (j > size (p2)) exit end if end do end subroutine subevt_set_p2_outgoing @ %def subevt_set_p2_beam @ %def subevt_set_p2_incoming @ %def subevt_set_p2_outgoing @ Set polarization for an entry <>= public :: subevt_polarize <>= subroutine subevt_polarize (subevt, i, h) type(subevt_t), intent(inout) :: subevt integer, intent(in) :: i, h call prt_polarize (subevt%prt(i), h) end subroutine subevt_polarize @ %def subevt_polarize @ Set color-flow indices for an entry <>= public :: subevt_colorize <>= subroutine subevt_colorize (subevt, i, col, acl) type(subevt_t), intent(inout) :: subevt integer, intent(in) :: i, col, acl if (col > 0 .and. acl > 0) then call prt_colorize (subevt%prt(i), [col], [acl]) else if (col > 0) then call prt_colorize (subevt%prt(i), [col], [integer ::]) else if (acl > 0) then call prt_colorize (subevt%prt(i), [integer ::], [acl]) else call prt_colorize (subevt%prt(i), [integer ::], [integer ::]) end if end subroutine subevt_colorize @ %def subevt_colorize @ \subsubsection{Accessing contents} Return true if the subevent has entries. <>= public :: subevt_is_nonempty <>= function subevt_is_nonempty (subevt) result (flag) logical :: flag type(subevt_t), intent(in) :: subevt flag = subevt%n_active /= 0 end function subevt_is_nonempty @ %def subevt_is_nonempty @ Return the number of entries <>= public :: subevt_get_length <>= function subevt_get_length (subevt) result (length) integer :: length type(subevt_t), intent(in) :: subevt length = subevt%n_active end function subevt_get_length @ %def subevt_get_length @ Return a specific particle. The index is not checked for validity. <>= public :: subevt_get_prt <>= function subevt_get_prt (subevt, i) result (prt) type(prt_t) :: prt type(subevt_t), intent(in) :: subevt integer, intent(in) :: i prt = subevt%prt(i) end function subevt_get_prt @ %def subevt_get_prt @ Return the partonic energy squared. We take the particles with flag [[PRT_INCOMING]] and compute their total invariant mass. <>= public :: subevt_get_sqrts_hat <>= function subevt_get_sqrts_hat (subevt) result (sqrts_hat) type(subevt_t), intent(in) :: subevt real(default) :: sqrts_hat type(vector4_t) :: p integer :: i do i = 1, subevt%n_active if (subevt%prt(i)%type == PRT_INCOMING) then p = p + prt_get_momentum (subevt%prt(i)) end if end do sqrts_hat = p ** 1 end function subevt_get_sqrts_hat @ %def subevt_get_sqrts_hat @ Return the number of incoming (outgoing) particles, respectively. Beam particles or composites are not counted. <>= public :: subevt_get_n_in public :: subevt_get_n_out <>= function subevt_get_n_in (subevt) result (n_in) type(subevt_t), intent(in) :: subevt integer :: n_in n_in = count (subevt%prt(:subevt%n_active)%type == PRT_INCOMING) end function subevt_get_n_in function subevt_get_n_out (subevt) result (n_out) type(subevt_t), intent(in) :: subevt integer :: n_out n_out = count (subevt%prt(:subevt%n_active)%type == PRT_OUTGOING) end function subevt_get_n_out @ %def subevt_get_n_in @ %def subevt_get_n_out @ <>= interface c_prt module procedure c_prt_from_subevt module procedure c_prt_array_from_subevt end interface @ %def c_prt <>= function c_prt_from_subevt (subevt, i) result (c_prt) type(c_prt_t) :: c_prt type(subevt_t), intent(in) :: subevt integer, intent(in) :: i c_prt = c_prt_from_prt (subevt%prt(i)) end function c_prt_from_subevt function c_prt_array_from_subevt (subevt) result (c_prt_array) type(subevt_t), intent(in) :: subevt type(c_prt_t), dimension(subevt%n_active) :: c_prt_array c_prt_array = c_prt_from_prt (subevt%prt(1:subevt%n_active)) end function c_prt_array_from_subevt @ %def c_prt_from_subevt @ %def c_prt_array_from_subevt @ \subsubsection{Operations with subevents} The join operation joins two subevents. When appending the elements of the second list, we check for each particle whether it is already in the first list. If yes, it is discarded. The result list should be initialized already. If a mask is present, it refers to the second subevent. Particles where the mask is not set are discarded. <>= public :: subevt_join <>= subroutine subevt_join (subevt, pl1, pl2, mask2) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl1, pl2 logical, dimension(:), intent(in), optional :: mask2 integer :: n1, n2, i, n n1 = pl1%n_active n2 = pl2%n_active call subevt_reset (subevt, n1 + n2) subevt%prt(:n1) = pl1%prt(:n1) n = n1 if (present (mask2)) then do i = 1, pl2%n_active if (mask2(i)) then if (disjoint (i)) then n = n + 1 subevt%prt(n) = pl2%prt(i) end if end if end do else do i = 1, pl2%n_active if (disjoint (i)) then n = n + 1 subevt%prt(n) = pl2%prt(i) end if end do end if subevt%n_active = n contains function disjoint (i) result (flag) integer, intent(in) :: i logical :: flag integer :: j do j = 1, pl1%n_active if (.not. are_disjoint (pl1%prt(j), pl2%prt(i))) then flag = .false. return end if end do flag = .true. end function disjoint end subroutine subevt_join @ %def subevt_join @ The combine operation makes a subevent whose entries are the result of adding (the momenta of) each pair of particles in the input lists. We trace the particles from which a particles is built by storing a [[src]] array. Each particle entry in the [[src]] list contains a list of indices which indicates its building blocks. The indices refer to an original list of particles. Index lists are sorted, and they contain no element more than once. We thus require that in a given pseudoparticle, each original particle occurs at most once. <>= public :: subevt_combine <>= subroutine subevt_combine (subevt, pl1, pl2, mask12) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl1, pl2 logical, dimension(:,:), intent(in), optional :: mask12 integer :: n1, n2, i1, i2, n, j logical :: ok n1 = pl1%n_active n2 = pl2%n_active call subevt_reset (subevt, n1 * n2) n = 1 do i1 = 1, n1 do i2 = 1, n2 if (present (mask12)) then ok = mask12(i1,i2) else ok = .true. end if if (ok) call prt_combine & (subevt%prt(n), pl1%prt(i1), pl2%prt(i2), ok) if (ok) then CHECK_DOUBLES: do j = 1, n - 1 if (subevt%prt(n) .match. subevt%prt(j)) then ok = .false.; exit CHECK_DOUBLES end if end do CHECK_DOUBLES if (ok) n = n + 1 end if end do end do subevt%n_active = n - 1 end subroutine subevt_combine @ %def subevt_combine @ The collect operation makes a single-entry subevent which results from combining (the momenta of) all particles in the input list. As above, the result does not contain an original particle more than once; this is checked for each particle when it is collected. Furthermore, each entry has a mask; where the mask is false, the entry is dropped. (Thus, if the input particles are already composite, there is some chance that the result depends on the order of the input list and is not as expected. This situation should be avoided.) <>= public :: subevt_collect <>= subroutine subevt_collect (subevt, pl1, mask1) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl1 logical, dimension(:), intent(in) :: mask1 type(prt_t) :: prt integer :: i logical :: ok call subevt_reset (subevt, 1) subevt%n_active = 0 do i = 1, pl1%n_active if (mask1(i)) then if (subevt%n_active == 0) then subevt%n_active = 1 subevt%prt(1) = pl1%prt(i) else call prt_combine (prt, subevt%prt(1), pl1%prt(i), ok) if (ok) subevt%prt(1) = prt end if end if end do end subroutine subevt_collect @ %def subevt_collect @ The cluster operation is similar to [[collect]], but applies a jet algorithm. The result is a subevent consisting of jets and, possibly, unclustered extra particles. As above, the result does not contain an original particle more than once; this is checked for each particle when it is collected. Furthermore, each entry has a mask; where the mask is false, the entry is dropped. The algorithm: first determine the (pseudo)particles that participate in the clustering. They should not overlap, and the mask entry must be set. We then cluster the particles, using the given jet definition. The result particles are retrieved from the cluster sequence. We still have to determine the source indices for each jet: for each input particle, we get the jet index. Accumulating the source entries for all particles that are part of a given jet, we derive the jet source entries. Finally, we delete the C structures that have been constructed by FastJet and its interface. <>= public :: subevt_cluster <>= subroutine subevt_cluster (subevt, pl1, dcut, mask1, jet_def, & keep_jets, exclusive) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl1 real(default), intent(in) :: dcut logical, dimension(:), intent(in) :: mask1 type(jet_definition_t), intent(in) :: jet_def logical, intent(in) :: keep_jets, exclusive integer, dimension(:), allocatable :: map, jet_index type(pseudojet_t), dimension(:), allocatable :: jet_in, jet_out type(pseudojet_vector_t) :: jv_in, jv_out type(cluster_sequence_t) :: cs integer :: i, n_src, n_active call map_prt_index (pl1, mask1, n_src, map) n_active = count (map /= 0) allocate (jet_in (n_active)) allocate (jet_index (n_active)) do i = 1, n_active call jet_in(i)%init (prt_get_momentum (pl1%prt(map(i)))) end do call jv_in%init (jet_in) call cs%init (jv_in, jet_def) if (exclusive) then jv_out = cs%exclusive_jets (dcut) else jv_out = cs%inclusive_jets () end if call cs%assign_jet_indices (jv_out, jet_index) allocate (jet_out (jv_out%size ())) jet_out = jv_out call fill_pseudojet (subevt, pl1, jet_out, jet_index, n_src, map) do i = 1, size (jet_out) call jet_out(i)%final () end do call jv_out%final () call cs%final () call jv_in%final () do i = 1, size (jet_in) call jet_in(i)%final () end do contains ! Uniquely combine sources and add map those new indices to the old ones subroutine map_prt_index (pl1, mask1, n_src, map) type(subevt_t), intent(in) :: pl1 logical, dimension(:), intent(in) :: mask1 integer, intent(out) :: n_src integer, dimension(:), allocatable, intent(out) :: map integer, dimension(:), allocatable :: src, src_tmp integer :: i allocate (src(0)) allocate (map (pl1%n_active), source = 0) n_active = 0 do i = 1, pl1%n_active if (.not. mask1(i)) cycle call combine_index_lists (src_tmp, src, pl1%prt(i)%src) if (.not. allocated (src_tmp)) cycle call move_alloc (from=src_tmp, to=src) n_active = n_active + 1 map(n_active) = i end do n_src = size (src) end subroutine map_prt_index ! Retrieve source(s) of a jet and fill corresponding subevent subroutine fill_pseudojet (subevt, pl1, jet_out, jet_index, n_src, map) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl1 type(pseudojet_t), dimension(:), intent(in) :: jet_out integer, dimension(:), intent(in) :: jet_index integer, dimension(:), intent(in) :: map integer, intent(in) :: n_src integer, dimension(n_src) :: src_fill integer :: i, jet, k, combined_pdg, pdg, n_quarks, n_src_fill logical :: is_b, is_c call subevt_reset (subevt, size (jet_out)) do jet = 1, size (jet_out) pdg = 0; src_fill = 0; n_src_fill = 0; combined_pdg = 0; n_quarks = 0 is_b = .false.; is_c = .false. PARTICLE: do i = 1, size (jet_index) if (jet_index(i) /= jet) cycle PARTICLE associate (prt => pl1%prt(map(i)), n_src_prt => size(pl1%prt(map(i))%src)) do k = 1, n_src_prt src_fill(n_src_fill + k) = prt%src(k) end do n_src_fill = n_src_fill + n_src_prt if (is_quark (prt%pdg)) then n_quarks = n_quarks + 1 if (.not. is_b) then if (abs (prt%pdg) == 5) then is_b = .true. is_c = .false. else if (abs (prt%pdg) == 4) then is_c = .true. end if end if if (combined_pdg == 0) combined_pdg = prt%pdg end if end associate end do PARTICLE if (keep_jets .and. n_quarks == 1) pdg = combined_pdg call prt_init_pseudojet (subevt%prt(jet), jet_out(jet), & src_fill(:n_src_fill), pdg, is_b, is_c) end do end subroutine fill_pseudojet end subroutine subevt_cluster @ %def subevt_cluster @ 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. <>= public :: subevt_recombine <>= 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), 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_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) 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 = pl end if else subevt = pl end if end if end subroutine subevt_recombine @ %def subevt_recombine @ Return a list of all particles for which the mask is true. <>= public :: subevt_select <>= subroutine subevt_select (subevt, pl, mask1) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl logical, dimension(:), intent(in) :: mask1 integer :: i, n call subevt_reset (subevt, pl%n_active) n = 0 do i = 1, pl%n_active if (mask1(i)) then n = n + 1 subevt%prt(n) = pl%prt(i) end if end do subevt%n_active = n end subroutine subevt_select @ %def subevt_select @ Return a subevent which consists of the single particle with specified [[index]]. If [[index]] is negative, count from the end. If it is out of bounds, return an empty list. <>= public :: subevt_extract <>= subroutine subevt_extract (subevt, pl, index) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl integer, intent(in) :: index if (index > 0) then if (index <= pl%n_active) then call subevt_reset (subevt, 1) subevt%prt(1) = pl%prt(index) else call subevt_reset (subevt, 0) end if else if (index < 0) then if (abs (index) <= pl%n_active) then call subevt_reset (subevt, 1) subevt%prt(1) = pl%prt(pl%n_active + 1 + index) else call subevt_reset (subevt, 0) end if else call subevt_reset (subevt, 0) end if end subroutine subevt_extract @ %def subevt_extract @ Return the list of particles sorted according to increasing values of the provided integer or real array. If no array is given, sort by PDG value. <>= public :: subevt_sort <>= interface subevt_sort module procedure subevt_sort_pdg module procedure subevt_sort_int module procedure subevt_sort_real end interface <>= subroutine subevt_sort_pdg (subevt, pl) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl integer :: n n = subevt%n_active call subevt_sort_int (subevt, pl, abs (3 * subevt%prt(:n)%pdg - 1)) end subroutine subevt_sort_pdg subroutine subevt_sort_int (subevt, pl, ival) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl integer, dimension(:), intent(in) :: ival call subevt_reset (subevt, pl%n_active) subevt%n_active = pl%n_active subevt%prt = pl%prt( order (ival) ) end subroutine subevt_sort_int subroutine subevt_sort_real (subevt, pl, rval) type(subevt_t), intent(inout) :: subevt type(subevt_t), intent(in) :: pl real(default), dimension(:), intent(in) :: rval integer :: i integer, dimension(size(rval)) :: idx call subevt_reset (subevt, pl%n_active) subevt%n_active = pl%n_active if (allocated (subevt%prt)) deallocate (subevt%prt) allocate (subevt%prt (size(pl%prt))) idx = order (rval) do i = 1, size (idx) subevt%prt(i) = pl%prt (idx(i)) end do end subroutine subevt_sort_real @ %def subevt_sort @ Return the list of particles which have any of the specified PDG codes (and optionally particle type: beam, incoming, outgoing). <>= public :: subevt_select_pdg_code <>= subroutine subevt_select_pdg_code (subevt, aval, subevt_in, prt_type) type(subevt_t), intent(inout) :: subevt type(pdg_array_t), intent(in) :: aval type(subevt_t), intent(in) :: subevt_in integer, intent(in), optional :: prt_type integer :: n_active, n_match logical, dimension(:), allocatable :: mask integer :: i, j n_active = subevt_in%n_active allocate (mask (n_active)) forall (i = 1:n_active) & mask(i) = aval .match. subevt_in%prt(i)%pdg if (present (prt_type)) & mask = mask .and. subevt_in%prt(:n_active)%type == prt_type n_match = count (mask) call subevt_reset (subevt, n_match) j = 0 do i = 1, n_active if (mask(i)) then j = j + 1 subevt%prt(j) = subevt_in%prt(i) end if end do end subroutine subevt_select_pdg_code @ %def subevt_select_pdg_code @ \subsection{Eliminate numerical noise} This is useful for testing purposes: set entries to zero that are smaller in absolute values than a given tolerance parameter. Note: instead of setting the tolerance in terms of EPSILON (kind-dependent), we fix it to $10^{-16}$, which is the typical value for double precision. The reason is that there are situations where intermediate representations (external libraries, files) are limited to double precision, even if the main program uses higher precision. <>= public :: pacify <>= interface pacify module procedure pacify_prt module procedure pacify_subevt end interface pacify @ %def pacify <>= subroutine pacify_prt (prt) class(prt_t), intent(inout) :: prt real(default) :: e e = max (1E-10_default * energy (prt%p), 1E-13_default) call pacify (prt%p, e) call pacify (prt%p2, 1E3_default * e) end subroutine pacify_prt subroutine pacify_subevt (subevt) class(subevt_t), intent(inout) :: subevt integer :: i do i = 1, subevt%n_active call pacify (subevt%prt(i)) end do end subroutine pacify_subevt @ %def pacify_prt @ %def pacify_subevt @ \clearpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Analysis tools} This module defines structures useful for data analysis. These include observables, histograms, and plots. Observables are quantities that are calculated and summed up event by event. At the end, one can compute the average and error. Histograms have their bins in addition to the observable properties. Histograms are usually written out in tables and displayed graphically. In plots, each record creates its own entry in a table. This can be used for scatter plots if called event by event, or for plotting dependencies on parameters if called once per integration run. Graphs are container for histograms and plots, which carry their own graphics options. The type layout is still somewhat obfuscated. This would become much simpler if type extension could be used. <<[[analysis.f90]]>>= <> module analysis <> <> use io_units use format_utils, only: quote_underscore, tex_format use system_defs, only: TAB use diagnostics use os_interface use ifiles <> <> <> <> <> <> contains <> end module analysis @ %def analysis @ \subsection{Output formats} These formats share a common field width (alignment). <>= character(*), parameter, public :: HISTOGRAM_HEAD_FORMAT = "1x,A15,3x" character(*), parameter, public :: HISTOGRAM_INTG_FORMAT = "3x,I9,3x" character(*), parameter, public :: HISTOGRAM_DATA_FORMAT = "ES19.12" @ %def HISTOGRAM_HEAD_FORMAT HISTOGRAM_INTG_FORMAT HISTOGRAM_DATA_FORMAT @ \subsection{Graph options} These parameters are used for displaying data. They apply to a whole graph, which may contain more than one plot element. The GAMELAN code chunks are part of both [[graph_options]] and [[drawing_options]]. The [[drawing_options]] copy is used in histograms and plots, also as graph elements. The [[graph_options]] copy is used for [[graph]] objects as a whole. Both copies are usually identical. <>= public :: graph_options_t <>= type :: graph_options_t private type(string_t) :: id type(string_t) :: title type(string_t) :: description type(string_t) :: x_label type(string_t) :: y_label integer :: width_mm = 130 integer :: height_mm = 90 logical :: x_log = .false. logical :: y_log = .false. real(default) :: x_min = 0 real(default) :: x_max = 1 real(default) :: y_min = 0 real(default) :: y_max = 1 logical :: x_min_set = .false. logical :: x_max_set = .false. logical :: y_min_set = .false. logical :: y_max_set = .false. type(string_t) :: gmlcode_bg type(string_t) :: gmlcode_fg end type graph_options_t @ %def graph_options_t @ Initialize the record, all strings are empty. The limits are undefined. <>= public :: graph_options_init <>= subroutine graph_options_init (graph_options) type(graph_options_t), intent(out) :: graph_options graph_options%id = "" graph_options%title = "" graph_options%description = "" graph_options%x_label = "" graph_options%y_label = "" graph_options%gmlcode_bg = "" graph_options%gmlcode_fg = "" end subroutine graph_options_init @ %def graph_options_init @ Set individual options. <>= public :: graph_options_set <>= subroutine graph_options_set (graph_options, id, & title, description, x_label, y_label, width_mm, height_mm, & x_log, y_log, x_min, x_max, y_min, y_max, & gmlcode_bg, gmlcode_fg) type(graph_options_t), intent(inout) :: graph_options type(string_t), intent(in), optional :: id type(string_t), intent(in), optional :: title type(string_t), intent(in), optional :: description type(string_t), intent(in), optional :: x_label, y_label integer, intent(in), optional :: width_mm, height_mm logical, intent(in), optional :: x_log, y_log real(default), intent(in), optional :: x_min, x_max, y_min, y_max type(string_t), intent(in), optional :: gmlcode_bg, gmlcode_fg if (present (id)) graph_options%id = id if (present (title)) graph_options%title = title if (present (description)) graph_options%description = description if (present (x_label)) graph_options%x_label = x_label if (present (y_label)) graph_options%y_label = y_label if (present (width_mm)) graph_options%width_mm = width_mm if (present (height_mm)) graph_options%height_mm = height_mm if (present (x_log)) graph_options%x_log = x_log if (present (y_log)) graph_options%y_log = y_log if (present (x_min)) graph_options%x_min = x_min if (present (x_max)) graph_options%x_max = x_max if (present (y_min)) graph_options%y_min = y_min if (present (y_max)) graph_options%y_max = y_max if (present (x_min)) graph_options%x_min_set = .true. if (present (x_max)) graph_options%x_max_set = .true. if (present (y_min)) graph_options%y_min_set = .true. if (present (y_max)) graph_options%y_max_set = .true. if (present (gmlcode_bg)) graph_options%gmlcode_bg = gmlcode_bg if (present (gmlcode_fg)) graph_options%gmlcode_fg = gmlcode_fg end subroutine graph_options_set @ %def graph_options_set @ Write a simple account of all options. <>= public :: graph_options_write <>= subroutine graph_options_write (gro, unit) type(graph_options_t), intent(in) :: gro integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) 1 format (A,1x,'"',A,'"') 2 format (A,1x,L1) 3 format (A,1x,ES19.12) 4 format (A,1x,I0) 5 format (A,1x,'[undefined]') write (u, 1) "title =", char (gro%title) write (u, 1) "description =", char (gro%description) write (u, 1) "x_label =", char (gro%x_label) write (u, 1) "y_label =", char (gro%y_label) write (u, 2) "x_log =", gro%x_log write (u, 2) "y_log =", gro%y_log if (gro%x_min_set) then write (u, 3) "x_min =", gro%x_min else write (u, 5) "x_min =" end if if (gro%x_max_set) then write (u, 3) "x_max =", gro%x_max else write (u, 5) "x_max =" end if if (gro%y_min_set) then write (u, 3) "y_min =", gro%y_min else write (u, 5) "y_min =" end if if (gro%y_max_set) then write (u, 3) "y_max =", gro%y_max else write (u, 5) "y_max =" end if write (u, 4) "width_mm =", gro%width_mm write (u, 4) "height_mm =", gro%height_mm write (u, 1) "gmlcode_bg =", char (gro%gmlcode_bg) write (u, 1) "gmlcode_fg =", char (gro%gmlcode_fg) end subroutine graph_options_write @ %def graph_options_write @ Write a \LaTeX\ header/footer for the analysis file. <>= subroutine graph_options_write_tex_header (gro, unit) type(graph_options_t), intent(in) :: gro integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) if (gro%title /= "") then write (u, "(A)") write (u, "(A)") "\section{" // char (gro%title) // "}" else write (u, "(A)") "\section{" // char (quote_underscore (gro%id)) // "}" end if if (gro%description /= "") then write (u, "(A)") char (gro%description) write (u, *) write (u, "(A)") "\vspace*{\baselineskip}" end if write (u, "(A)") "\vspace*{\baselineskip}" write (u, "(A)") "\unitlength 1mm" write (u, "(A,I0,',',I0,A)") & "\begin{gmlgraph*}(", & gro%width_mm, gro%height_mm, & ")[dat]" end subroutine graph_options_write_tex_header subroutine graph_options_write_tex_footer (gro, unit) type(graph_options_t), intent(in) :: gro integer, intent(in), optional :: unit integer :: u, width, height width = gro%width_mm - 10 height = gro%height_mm - 10 u = given_output_unit (unit) write (u, "(A)") " begingmleps ""Whizard-Logo.eps"";" write (u, "(A,I0,A,I0,A)") & " base := (", width, "*unitlength,", height, "*unitlength);" write (u, "(A)") " height := 9.6*unitlength;" write (u, "(A)") " width := 11.2*unitlength;" write (u, "(A)") " endgmleps;" write (u, "(A)") "\end{gmlgraph*}" end subroutine graph_options_write_tex_footer @ %def graph_options_write_tex_header @ %def graph_options_write_tex_footer @ Return the analysis object ID. <>= function graph_options_get_id (gro) result (id) type(string_t) :: id type(graph_options_t), intent(in) :: gro id = gro%id end function graph_options_get_id @ %def graph_options_get_id @ Create an appropriate [[setup]] command (linear/log). <>= function graph_options_get_gml_setup (gro) result (cmd) type(string_t) :: cmd type(graph_options_t), intent(in) :: gro type(string_t) :: x_str, y_str if (gro%x_log) then x_str = "log" else x_str = "linear" end if if (gro%y_log) then y_str = "log" else y_str = "linear" end if cmd = "setup (" // x_str // ", " // y_str // ");" end function graph_options_get_gml_setup @ %def graph_options_get_gml_setup @ Return the labels in GAMELAN form. <>= function graph_options_get_gml_x_label (gro) result (cmd) type(string_t) :: cmd type(graph_options_t), intent(in) :: gro cmd = 'label.bot (<' // '<' // gro%x_label // '>' // '>, out);' end function graph_options_get_gml_x_label function graph_options_get_gml_y_label (gro) result (cmd) type(string_t) :: cmd type(graph_options_t), intent(in) :: gro cmd = 'label.ulft (<' // '<' // gro%y_label // '>' // '>, out);' end function graph_options_get_gml_y_label @ %def graph_options_get_gml_x_label @ %def graph_options_get_gml_y_label @ Create an appropriate [[graphrange]] statement for the given graph options. Where the graph options are not set, use the supplied arguments, if any, otherwise set the undefined value. <>= function graph_options_get_gml_graphrange & (gro, x_min, x_max, y_min, y_max) result (cmd) type(string_t) :: cmd type(graph_options_t), intent(in) :: gro real(default), intent(in), optional :: x_min, x_max, y_min, y_max type(string_t) :: x_min_str, x_max_str, y_min_str, y_max_str character(*), parameter :: fmt = "(ES15.8)" if (gro%x_min_set) then x_min_str = "#" // trim (adjustl (real2string (gro%x_min, fmt))) else if (present (x_min)) then x_min_str = "#" // trim (adjustl (real2string (x_min, fmt))) else x_min_str = "??" end if if (gro%x_max_set) then x_max_str = "#" // trim (adjustl (real2string (gro%x_max, fmt))) else if (present (x_max)) then x_max_str = "#" // trim (adjustl (real2string (x_max, fmt))) else x_max_str = "??" end if if (gro%y_min_set) then y_min_str = "#" // trim (adjustl (real2string (gro%y_min, fmt))) else if (present (y_min)) then y_min_str = "#" // trim (adjustl (real2string (y_min, fmt))) else y_min_str = "??" end if if (gro%y_max_set) then y_max_str = "#" // trim (adjustl (real2string (gro%y_max, fmt))) else if (present (y_max)) then y_max_str = "#" // trim (adjustl (real2string (y_max, fmt))) else y_max_str = "??" end if cmd = "graphrange (" // x_min_str // ", " // y_min_str // "), " & // "(" // x_max_str // ", " // y_max_str // ");" end function graph_options_get_gml_graphrange @ %def graph_options_get_gml_graphrange @ Get extra GAMELAN code to be executed before and after the usual drawing commands. <>= function graph_options_get_gml_bg_command (gro) result (cmd) type(string_t) :: cmd type(graph_options_t), intent(in) :: gro cmd = gro%gmlcode_bg end function graph_options_get_gml_bg_command function graph_options_get_gml_fg_command (gro) result (cmd) type(string_t) :: cmd type(graph_options_t), intent(in) :: gro cmd = gro%gmlcode_fg end function graph_options_get_gml_fg_command @ %def graph_options_get_gml_bg_command @ %def graph_options_get_gml_fg_command @ Append the header for generic data output in ifile format. We print only labels, not graphics parameters. <>= subroutine graph_options_get_header (pl, header, comment) type(graph_options_t), intent(in) :: pl type(ifile_t), intent(inout) :: header type(string_t), intent(in), optional :: comment type(string_t) :: c if (present (comment)) then c = comment else c = "" end if call ifile_append (header, & c // "ID: " // pl%id) call ifile_append (header, & c // "title: " // pl%title) call ifile_append (header, & c // "description: " // pl%description) call ifile_append (header, & c // "x axis label: " // pl%x_label) call ifile_append (header, & c // "y axis label: " // pl%y_label) end subroutine graph_options_get_header @ %def graph_options_get_header @ \subsection{Drawing options} These options apply to an individual graph element (histogram or plot). <>= public :: drawing_options_t <>= type :: drawing_options_t type(string_t) :: dataset logical :: with_hbars = .false. logical :: with_base = .false. logical :: piecewise = .false. logical :: fill = .false. logical :: draw = .false. logical :: err = .false. logical :: symbols = .false. type(string_t) :: fill_options type(string_t) :: draw_options type(string_t) :: err_options type(string_t) :: symbol type(string_t) :: gmlcode_bg type(string_t) :: gmlcode_fg end type drawing_options_t @ %def drawing_options_t @ Write a simple account of all options. <>= public :: drawing_options_write <>= subroutine drawing_options_write (dro, unit) type(drawing_options_t), intent(in) :: dro integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit) 1 format (A,1x,'"',A,'"') 2 format (A,1x,L1) write (u, 2) "with_hbars =", dro%with_hbars write (u, 2) "with_base =", dro%with_base write (u, 2) "piecewise =", dro%piecewise write (u, 2) "fill =", dro%fill write (u, 2) "draw =", dro%draw write (u, 2) "err =", dro%err write (u, 2) "symbols =", dro%symbols write (u, 1) "fill_options=", char (dro%fill_options) write (u, 1) "draw_options=", char (dro%draw_options) write (u, 1) "err_options =", char (dro%err_options) write (u, 1) "symbol =", char (dro%symbol) write (u, 1) "gmlcode_bg =", char (dro%gmlcode_bg) write (u, 1) "gmlcode_fg =", char (dro%gmlcode_fg) end subroutine drawing_options_write @ %def drawing_options_write @ Init with empty strings and default options, appropriate for either histogram or plot. <>= public :: drawing_options_init_histogram public :: drawing_options_init_plot <>= subroutine drawing_options_init_histogram (dro) type(drawing_options_t), intent(out) :: dro dro%dataset = "dat" dro%with_hbars = .true. dro%with_base = .true. dro%piecewise = .true. dro%fill = .true. dro%draw = .true. dro%fill_options = "withcolor col.default" dro%draw_options = "" dro%err_options = "" dro%symbol = "fshape(circle scaled 1mm)()" dro%gmlcode_bg = "" dro%gmlcode_fg = "" end subroutine drawing_options_init_histogram subroutine drawing_options_init_plot (dro) type(drawing_options_t), intent(out) :: dro dro%dataset = "dat" dro%draw = .true. dro%fill_options = "withcolor col.default" dro%draw_options = "" dro%err_options = "" dro%symbol = "fshape(circle scaled 1mm)()" dro%gmlcode_bg = "" dro%gmlcode_fg = "" end subroutine drawing_options_init_plot @ %def drawing_options_init_histogram @ %def drawing_options_init_plot @ Set individual options. <>= public :: drawing_options_set <>= subroutine drawing_options_set (dro, dataset, & with_hbars, with_base, piecewise, fill, draw, err, symbols, & fill_options, draw_options, err_options, symbol, & gmlcode_bg, gmlcode_fg) type(drawing_options_t), intent(inout) :: dro type(string_t), intent(in), optional :: dataset logical, intent(in), optional :: with_hbars, with_base, piecewise logical, intent(in), optional :: fill, draw, err, symbols type(string_t), intent(in), optional :: fill_options, draw_options type(string_t), intent(in), optional :: err_options, symbol type(string_t), intent(in), optional :: gmlcode_bg, gmlcode_fg if (present (dataset)) dro%dataset = dataset if (present (with_hbars)) dro%with_hbars = with_hbars if (present (with_base)) dro%with_base = with_base if (present (piecewise)) dro%piecewise = piecewise if (present (fill)) dro%fill = fill if (present (draw)) dro%draw = draw if (present (err)) dro%err = err if (present (symbols)) dro%symbols = symbols if (present (fill_options)) dro%fill_options = fill_options if (present (draw_options)) dro%draw_options = draw_options if (present (err_options)) dro%err_options = err_options if (present (symbol)) dro%symbol = symbol if (present (gmlcode_bg)) dro%gmlcode_bg = gmlcode_bg if (present (gmlcode_fg)) dro%gmlcode_fg = gmlcode_fg end subroutine drawing_options_set @ %def drawing_options_set @ There are sepate commands for drawing the curve and for drawing errors. The symbols are applied to the latter. First of all, we may have to compute a baseline: <>= function drawing_options_get_calc_command (dro) result (cmd) type(string_t) :: cmd type(drawing_options_t), intent(in) :: dro if (dro%with_base) then cmd = "calculate " // dro%dataset // ".base (" // dro%dataset // ") " & // "(x, #0);" else cmd = "" end if end function drawing_options_get_calc_command @ %def drawing_options_get_calc_command @ Return the drawing command. <>= function drawing_options_get_draw_command (dro) result (cmd) type(string_t) :: cmd type(drawing_options_t), intent(in) :: dro if (dro%fill) then cmd = "fill" else if (dro%draw) then cmd = "draw" else cmd = "" end if if (dro%fill .or. dro%draw) then if (dro%piecewise) cmd = cmd // " piecewise" if (dro%draw .and. dro%with_base) cmd = cmd // " cyclic" cmd = cmd // " from (" // dro%dataset if (dro%with_base) then if (dro%piecewise) then cmd = cmd // ", " // dro%dataset // ".base/\" ! " else cmd = cmd // " ~ " // dro%dataset // ".base\" ! " end if end if cmd = cmd // ")" if (dro%fill) then cmd = cmd // " " // dro%fill_options if (dro%draw) cmd = cmd // " outlined" end if if (dro%draw) cmd = cmd // " " // dro%draw_options cmd = cmd // ";" end if end function drawing_options_get_draw_command @ %def drawing_options_get_draw_command @ The error command draws error bars, if any. <>= function drawing_options_get_err_command (dro) result (cmd) type(string_t) :: cmd type(drawing_options_t), intent(in) :: dro if (dro%err) then cmd = "draw piecewise " & // "from (" // dro%dataset // ".err)" & // " " // dro%err_options // ";" else cmd = "" end if end function drawing_options_get_err_command @ %def drawing_options_get_err_command @ The symbol command draws symbols, if any. <>= function drawing_options_get_symb_command (dro) result (cmd) type(string_t) :: cmd type(drawing_options_t), intent(in) :: dro if (dro%symbols) then cmd = "phantom" & // " from (" // dro%dataset // ")" & // " withsymbol (" // dro%symbol // ");" else cmd = "" end if end function drawing_options_get_symb_command @ %def drawing_options_get_symb_command @ Get extra GAMELAN code to be executed before and after the usual drawing commands. <>= function drawing_options_get_gml_bg_command (dro) result (cmd) type(string_t) :: cmd type(drawing_options_t), intent(in) :: dro cmd = dro%gmlcode_bg end function drawing_options_get_gml_bg_command function drawing_options_get_gml_fg_command (dro) result (cmd) type(string_t) :: cmd type(drawing_options_t), intent(in) :: dro cmd = dro%gmlcode_fg end function drawing_options_get_gml_fg_command @ %def drawing_options_get_gml_bg_command @ %def drawing_options_get_gml_fg_command @ \subsection{Observables} The observable type holds the accumulated observable values and weight sums which are necessary for proper averaging. <>= type :: observable_t private real(default) :: sum_values = 0 real(default) :: sum_squared_values = 0 real(default) :: sum_weights = 0 real(default) :: sum_squared_weights = 0 integer :: count = 0 type(string_t) :: obs_label type(string_t) :: obs_unit type(graph_options_t) :: graph_options end type observable_t @ %def observable_t @ Initialize with defined properties <>= subroutine observable_init (obs, obs_label, obs_unit, graph_options) type(observable_t), intent(out) :: obs type(string_t), intent(in), optional :: obs_label, obs_unit type(graph_options_t), intent(in), optional :: graph_options if (present (obs_label)) then obs%obs_label = obs_label else obs%obs_label = "" end if if (present (obs_unit)) then obs%obs_unit = obs_unit else obs%obs_unit = "" end if if (present (graph_options)) then obs%graph_options = graph_options else call graph_options_init (obs%graph_options) end if end subroutine observable_init @ %def observable_init @ Reset all numeric entries. <>= subroutine observable_clear (obs) type(observable_t), intent(inout) :: obs obs%sum_values = 0 obs%sum_squared_values = 0 obs%sum_weights = 0 obs%sum_squared_weights = 0 obs%count = 0 end subroutine observable_clear @ %def observable_clear @ Record a value. Always successful for observables. <>= interface observable_record_value module procedure observable_record_value_unweighted module procedure observable_record_value_weighted end interface <>= subroutine observable_record_value_unweighted (obs, value, success) type(observable_t), intent(inout) :: obs real(default), intent(in) :: value logical, intent(out), optional :: success obs%sum_values = obs%sum_values + value obs%sum_squared_values = obs%sum_squared_values + value**2 obs%sum_weights = obs%sum_weights + 1 obs%sum_squared_weights = obs%sum_squared_weights + 1 obs%count = obs%count + 1 if (present (success)) success = .true. end subroutine observable_record_value_unweighted subroutine observable_record_value_weighted (obs, value, weight, success) type(observable_t), intent(inout) :: obs real(default), intent(in) :: value, weight logical, intent(out), optional :: success obs%sum_values = obs%sum_values + value * weight obs%sum_squared_values = obs%sum_squared_values + value**2 * weight obs%sum_weights = obs%sum_weights + weight obs%sum_squared_weights = obs%sum_squared_weights + weight**2 obs%count = obs%count + 1 if (present (success)) success = .true. end subroutine observable_record_value_weighted @ %def observable_record_value @ Here are the statistics formulas: \begin{enumerate} \item Unweighted case: Given a sample of $n$ values $x_i$, the average is \begin{equation} \langle x \rangle = \frac{\sum x_i}{n} \end{equation} and the error estimate \begin{align} \Delta x &= \sqrt{\frac{1}{n-1}\langle{\sum(x_i - \langle x\rangle)^2}} \\ &= \sqrt{\frac{1}{n-1} \left(\frac{\sum x_i^2}{n} - \frac{(\sum x_i)^2}{n^2}\right)} \end{align} \item Weighted case: Instead of weight 1, each event comes with weight $w_i$. \begin{equation} \langle x \rangle = \frac{\sum x_i w_i}{\sum w_i} \end{equation} and \begin{equation} \Delta x = \sqrt{\frac{1}{n-1} \left(\frac{\sum x_i^2 w_i}{\sum w_i} - \frac{(\sum x_i w_i)^2}{(\sum w_i)^2}\right)} \end{equation} For $w_i=1$, this specializes to the previous formula. \end{enumerate} <>= function observable_get_n_entries (obs) result (n) integer :: n type(observable_t), intent(in) :: obs n = obs%count end function observable_get_n_entries function observable_get_average (obs) result (avg) real(default) :: avg type(observable_t), intent(in) :: obs if (obs%sum_weights /= 0) then avg = obs%sum_values / obs%sum_weights else avg = 0 end if end function observable_get_average function observable_get_error (obs) result (err) real(default) :: err type(observable_t), intent(in) :: obs real(default) :: var, n if (obs%sum_weights /= 0) then select case (obs%count) case (0:1) err = 0 case default n = obs%count var = obs%sum_squared_values / obs%sum_weights & - (obs%sum_values / obs%sum_weights) ** 2 err = sqrt (max (var, 0._default) / (n - 1)) end select else err = 0 end if end function observable_get_error @ %def observable_get_n_entries @ %def observable_get_sum @ %def observable_get_average @ %def observable_get_error @ Write label and/or physical unit to a string. <>= function observable_get_label (obs, wl, wu) result (string) type(string_t) :: string type(observable_t), intent(in) :: obs logical, intent(in) :: wl, wu type(string_t) :: obs_label, obs_unit if (wl) then if (obs%obs_label /= "") then obs_label = obs%obs_label else obs_label = "\textrm{Observable}" end if else obs_label = "" end if if (wu) then if (obs%obs_unit /= "") then if (wl) then obs_unit = "\;[" // obs%obs_unit // "]" else obs_unit = obs%obs_unit end if else obs_unit = "" end if else obs_unit = "" end if string = obs_label // obs_unit end function observable_get_label @ %def observable_get_label @ \subsection{Output} <>= subroutine observable_write (obs, unit) type(observable_t), intent(in) :: obs integer, intent(in), optional :: unit real(default) :: avg, err, relerr integer :: n integer :: u u = given_output_unit (unit); if (u < 0) return avg = observable_get_average (obs) err = observable_get_error (obs) if (avg /= 0) then relerr = err / abs (avg) else relerr = 0 end if n = observable_get_n_entries (obs) if (obs%graph_options%title /= "") then write (u, "(A,1x,3A)") & "title =", '"', char (obs%graph_options%title), '"' end if if (obs%graph_options%title /= "") then write (u, "(A,1x,3A)") & "description =", '"', char (obs%graph_options%description), '"' end if write (u, "(A,1x," // HISTOGRAM_DATA_FORMAT // ")", advance = "no") & "average =", avg call write_unit () write (u, "(A,1x," // HISTOGRAM_DATA_FORMAT // ")", advance = "no") & "error[abs] =", err call write_unit () write (u, "(A,1x," // HISTOGRAM_DATA_FORMAT // ")") & "error[rel] =", relerr write (u, "(A,1x,I0)") & "n_entries =", n contains subroutine write_unit () if (obs%obs_unit /= "") then write (u, "(1x,A)") char (obs%obs_unit) else write (u, *) end if end subroutine write_unit end subroutine observable_write @ %def observable_write @ \LaTeX\ output. <>= subroutine observable_write_driver (obs, unit, write_heading) type(observable_t), intent(in) :: obs integer, intent(in), optional :: unit logical, intent(in), optional :: write_heading real(default) :: avg, err integer :: n_digits logical :: heading integer :: u u = given_output_unit (unit); if (u < 0) return heading = .true.; if (present (write_heading)) heading = write_heading avg = observable_get_average (obs) err = observable_get_error (obs) if (avg /= 0 .and. err /= 0) then n_digits = max (2, 2 - int (log10 (abs (err / real (avg, default))))) else if (avg /= 0) then n_digits = 100 else n_digits = 1 end if if (heading) then write (u, "(A)") if (obs%graph_options%title /= "") then write (u, "(A)") "\section{" // char (obs%graph_options%title) & // "}" else write (u, "(A)") "\section{Observable}" end if if (obs%graph_options%description /= "") then write (u, "(A)") char (obs%graph_options%description) write (u, *) end if write (u, "(A)") "\begin{flushleft}" end if write (u, "(A)", advance="no") " $\langle{" ! $ sign write (u, "(A)", advance="no") char (observable_get_label (obs, wl=.true., wu=.false.)) write (u, "(A)", advance="no") "}\rangle = " write (u, "(A)", advance="no") char (tex_format (avg, n_digits)) write (u, "(A)", advance="no") "\pm" write (u, "(A)", advance="no") char (tex_format (err, 2)) write (u, "(A)", advance="no") "\;{" write (u, "(A)", advance="no") char (observable_get_label (obs, wl=.false., wu=.true.)) write (u, "(A)") "}" write (u, "(A)", advance="no") " \quad[n_{\text{entries}} = " write (u, "(I0)",advance="no") observable_get_n_entries (obs) write (u, "(A)") "]$" ! $ fool Emacs' noweb mode if (heading) then write (u, "(A)") "\end{flushleft}" end if end subroutine observable_write_driver @ %def observable_write_driver @ \subsection{Histograms} \subsubsection{Bins} <>= type :: bin_t private real(default) :: midpoint = 0 real(default) :: width = 0 real(default) :: sum_weights = 0 real(default) :: sum_squared_weights = 0 real(default) :: sum_excess_weights = 0 integer :: count = 0 end type bin_t @ %def bin_t <>= subroutine bin_init (bin, midpoint, width) type(bin_t), intent(out) :: bin real(default), intent(in) :: midpoint, width bin%midpoint = midpoint bin%width = width end subroutine bin_init @ %def bin_init <>= elemental subroutine bin_clear (bin) type(bin_t), intent(inout) :: bin bin%sum_weights = 0 bin%sum_squared_weights = 0 bin%sum_excess_weights = 0 bin%count = 0 end subroutine bin_clear @ %def bin_clear <>= subroutine bin_record_value (bin, normalize, weight, excess) type(bin_t), intent(inout) :: bin logical, intent(in) :: normalize real(default), intent(in) :: weight real(default), intent(in), optional :: excess real(default) :: w, e if (normalize) then if (bin%width /= 0) then w = weight / bin%width if (present (excess)) e = excess / bin%width else w = 0 if (present (excess)) e = 0 end if else w = weight if (present (excess)) e = excess end if bin%sum_weights = bin%sum_weights + w bin%sum_squared_weights = bin%sum_squared_weights + w ** 2 if (present (excess)) & bin%sum_excess_weights = bin%sum_excess_weights + abs (e) bin%count = bin%count + 1 end subroutine bin_record_value @ %def bin_record_value <>= function bin_get_midpoint (bin) result (x) real(default) :: x type(bin_t), intent(in) :: bin x = bin%midpoint end function bin_get_midpoint function bin_get_width (bin) result (w) real(default) :: w type(bin_t), intent(in) :: bin w = bin%width end function bin_get_width function bin_get_n_entries (bin) result (n) integer :: n type(bin_t), intent(in) :: bin n = bin%count end function bin_get_n_entries function bin_get_sum (bin) result (s) real(default) :: s type(bin_t), intent(in) :: bin s = bin%sum_weights end function bin_get_sum function bin_get_error (bin) result (err) real(default) :: err type(bin_t), intent(in) :: bin err = sqrt (bin%sum_squared_weights) end function bin_get_error function bin_get_excess (bin) result (excess) real(default) :: excess type(bin_t), intent(in) :: bin excess = bin%sum_excess_weights end function bin_get_excess @ %def bin_get_midpoint @ %def bin_get_width @ %def bin_get_n_entries @ %def bin_get_sum @ %def bin_get_error @ %def bin_get_excess <>= subroutine bin_write_header (unit) integer, intent(in), optional :: unit character(120) :: buffer integer :: u u = given_output_unit (unit); if (u < 0) return write (buffer, "(A,4(1x," //HISTOGRAM_HEAD_FORMAT // "),2x,A)") & "#", "bin midpoint", "value ", "error ", & "excess ", "n" write (u, "(A)") trim (buffer) end subroutine bin_write_header subroutine bin_write (bin, unit) type(bin_t), intent(in) :: bin integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(1x,4(1x," // HISTOGRAM_DATA_FORMAT // "),2x,I0)") & bin_get_midpoint (bin), & bin_get_sum (bin), & bin_get_error (bin), & bin_get_excess (bin), & bin_get_n_entries (bin) end subroutine bin_write @ %def bin_write_header @ %def bin_write @ \subsubsection{Histograms} <>= type :: histogram_t private real(default) :: lower_bound = 0 real(default) :: upper_bound = 0 real(default) :: width = 0 integer :: n_bins = 0 logical :: normalize_bins = .false. type(observable_t) :: obs type(observable_t) :: obs_within_bounds type(bin_t) :: underflow type(bin_t), dimension(:), allocatable :: bin type(bin_t) :: overflow type(graph_options_t) :: graph_options type(drawing_options_t) :: drawing_options end type histogram_t @ %def histogram_t @ \subsubsection{Initializer/finalizer} Initialize a histogram. We may provide either the bin width or the number of bins. A finalizer is not needed, since the histogram contains no pointer (sub)components. <>= interface histogram_init module procedure histogram_init_n_bins module procedure histogram_init_bin_width end interface <>= subroutine histogram_init_n_bins (h, id, & lower_bound, upper_bound, n_bins, normalize_bins, & obs_label, obs_unit, graph_options, drawing_options) type(histogram_t), intent(out) :: h type(string_t), intent(in) :: id real(default), intent(in) :: lower_bound, upper_bound integer, intent(in) :: n_bins logical, intent(in) :: normalize_bins type(string_t), intent(in), optional :: obs_label, obs_unit type(graph_options_t), intent(in), optional :: graph_options type(drawing_options_t), intent(in), optional :: drawing_options real(default) :: bin_width integer :: i call observable_init (h%obs_within_bounds, obs_label, obs_unit) call observable_init (h%obs, obs_label, obs_unit) h%lower_bound = lower_bound h%upper_bound = upper_bound h%n_bins = max (n_bins, 1) h%width = h%upper_bound - h%lower_bound h%normalize_bins = normalize_bins bin_width = h%width / h%n_bins allocate (h%bin (h%n_bins)) call bin_init (h%underflow, h%lower_bound, 0._default) do i = 1, h%n_bins call bin_init (h%bin(i), & h%lower_bound - bin_width/2 + i * bin_width, bin_width) end do call bin_init (h%overflow, h%upper_bound, 0._default) if (present (graph_options)) then h%graph_options = graph_options else call graph_options_init (h%graph_options) end if call graph_options_set (h%graph_options, id = id) if (present (drawing_options)) then h%drawing_options = drawing_options else call drawing_options_init_histogram (h%drawing_options) end if end subroutine histogram_init_n_bins subroutine histogram_init_bin_width (h, id, & lower_bound, upper_bound, bin_width, normalize_bins, & obs_label, obs_unit, graph_options, drawing_options) type(histogram_t), intent(out) :: h type(string_t), intent(in) :: id real(default), intent(in) :: lower_bound, upper_bound, bin_width logical, intent(in) :: normalize_bins type(string_t), intent(in), optional :: obs_label, obs_unit type(graph_options_t), intent(in), optional :: graph_options type(drawing_options_t), intent(in), optional :: drawing_options integer :: n_bins if (bin_width /= 0) then n_bins = nint ((upper_bound - lower_bound) / bin_width) else n_bins = 1 end if call histogram_init_n_bins (h, id, & lower_bound, upper_bound, n_bins, normalize_bins, & obs_label, obs_unit, graph_options, drawing_options) end subroutine histogram_init_bin_width @ %def histogram_init @ Initialize a histogram by copying another one. Since [[h]] has no pointer (sub)components, intrinsic assignment is sufficient. Optionally, we replace the drawing options. <>= subroutine histogram_init_histogram (h, h_in, drawing_options) type(histogram_t), intent(out) :: h type(histogram_t), intent(in) :: h_in type(drawing_options_t), intent(in), optional :: drawing_options h = h_in if (present (drawing_options)) then h%drawing_options = drawing_options end if end subroutine histogram_init_histogram @ %def histogram_init_histogram @ \subsubsection{Fill histograms} Clear the histogram contents, but do not modify the structure. <>= subroutine histogram_clear (h) type(histogram_t), intent(inout) :: h call observable_clear (h%obs) call observable_clear (h%obs_within_bounds) call bin_clear (h%underflow) if (allocated (h%bin)) call bin_clear (h%bin) call bin_clear (h%overflow) end subroutine histogram_clear @ %def histogram_clear @ Record a value. Successful if the value is within bounds, otherwise it is recorded as under-/overflow. Optionally, we may provide an excess weight that could be returned by the unweighting procedure. <>= subroutine histogram_record_value_unweighted (h, value, excess, success) type(histogram_t), intent(inout) :: h real(default), intent(in) :: value real(default), intent(in), optional :: excess logical, intent(out), optional :: success integer :: i_bin call observable_record_value (h%obs, value) if (h%width /= 0) then i_bin = floor (((value - h%lower_bound) / h%width) * h%n_bins) + 1 else i_bin = 0 end if if (i_bin <= 0) then call bin_record_value (h%underflow, .false., 1._default, excess) if (present (success)) success = .false. else if (i_bin <= h%n_bins) then call observable_record_value (h%obs_within_bounds, value) call bin_record_value & (h%bin(i_bin), h%normalize_bins, 1._default, excess) if (present (success)) success = .true. else call bin_record_value (h%overflow, .false., 1._default, excess) if (present (success)) success = .false. end if end subroutine histogram_record_value_unweighted @ %def histogram_record_value_unweighted @ Weighted events: analogous, but no excess weight. <>= subroutine histogram_record_value_weighted (h, value, weight, success) type(histogram_t), intent(inout) :: h real(default), intent(in) :: value, weight logical, intent(out), optional :: success integer :: i_bin call observable_record_value (h%obs, value, weight) if (h%width /= 0) then i_bin = floor (((value - h%lower_bound) / h%width) * h%n_bins) + 1 else i_bin = 0 end if if (i_bin <= 0) then call bin_record_value (h%underflow, .false., weight) if (present (success)) success = .false. else if (i_bin <= h%n_bins) then call observable_record_value (h%obs_within_bounds, value, weight) call bin_record_value (h%bin(i_bin), h%normalize_bins, weight) if (present (success)) success = .true. else call bin_record_value (h%overflow, .false., weight) if (present (success)) success = .false. end if end subroutine histogram_record_value_weighted @ %def histogram_record_value_weighted @ \subsubsection{Access contents} Inherited from the observable component (all-over average etc.) <>= function histogram_get_n_entries (h) result (n) integer :: n type(histogram_t), intent(in) :: h n = observable_get_n_entries (h%obs) end function histogram_get_n_entries function histogram_get_average (h) result (avg) real(default) :: avg type(histogram_t), intent(in) :: h avg = observable_get_average (h%obs) end function histogram_get_average function histogram_get_error (h) result (err) real(default) :: err type(histogram_t), intent(in) :: h err = observable_get_error (h%obs) end function histogram_get_error @ %def histogram_get_n_entries @ %def histogram_get_average @ %def histogram_get_error @ Analogous, but applied only to events within bounds. <>= function histogram_get_n_entries_within_bounds (h) result (n) integer :: n type(histogram_t), intent(in) :: h n = observable_get_n_entries (h%obs_within_bounds) end function histogram_get_n_entries_within_bounds function histogram_get_average_within_bounds (h) result (avg) real(default) :: avg type(histogram_t), intent(in) :: h avg = observable_get_average (h%obs_within_bounds) end function histogram_get_average_within_bounds function histogram_get_error_within_bounds (h) result (err) real(default) :: err type(histogram_t), intent(in) :: h err = observable_get_error (h%obs_within_bounds) end function histogram_get_error_within_bounds @ %def histogram_get_n_entries_within_bounds @ %def histogram_get_average_within_bounds @ %def histogram_get_error_within_bounds Get the number of bins <>= function histogram_get_n_bins (h) result (n) type(histogram_t), intent(in) :: h integer :: n n = h%n_bins end function histogram_get_n_bins @ %def histogram_get_n_bins @ Check bins. If the index is zero or above the limit, return the results for underflow or overflow, respectively. <>= function histogram_get_n_entries_for_bin (h, i) result (n) integer :: n type(histogram_t), intent(in) :: h integer, intent(in) :: i if (i <= 0) then n = bin_get_n_entries (h%underflow) else if (i <= h%n_bins) then n = bin_get_n_entries (h%bin(i)) else n = bin_get_n_entries (h%overflow) end if end function histogram_get_n_entries_for_bin function histogram_get_sum_for_bin (h, i) result (avg) real(default) :: avg type(histogram_t), intent(in) :: h integer, intent(in) :: i if (i <= 0) then avg = bin_get_sum (h%underflow) else if (i <= h%n_bins) then avg = bin_get_sum (h%bin(i)) else avg = bin_get_sum (h%overflow) end if end function histogram_get_sum_for_bin function histogram_get_error_for_bin (h, i) result (err) real(default) :: err type(histogram_t), intent(in) :: h integer, intent(in) :: i if (i <= 0) then err = bin_get_error (h%underflow) else if (i <= h%n_bins) then err = bin_get_error (h%bin(i)) else err = bin_get_error (h%overflow) end if end function histogram_get_error_for_bin function histogram_get_excess_for_bin (h, i) result (err) real(default) :: err type(histogram_t), intent(in) :: h integer, intent(in) :: i if (i <= 0) then err = bin_get_excess (h%underflow) else if (i <= h%n_bins) then err = bin_get_excess (h%bin(i)) else err = bin_get_excess (h%overflow) end if end function histogram_get_excess_for_bin @ %def histogram_get_n_entries_for_bin @ %def histogram_get_sum_for_bin @ %def histogram_get_error_for_bin @ %def histogram_get_excess_for_bin @ Return a pointer to the graph options. <>= function histogram_get_graph_options_ptr (h) result (ptr) type(graph_options_t), pointer :: ptr type(histogram_t), intent(in), target :: h ptr => h%graph_options end function histogram_get_graph_options_ptr @ %def histogram_get_graph_options_ptr @ Return a pointer to the drawing options. <>= function histogram_get_drawing_options_ptr (h) result (ptr) type(drawing_options_t), pointer :: ptr type(histogram_t), intent(in), target :: h ptr => h%drawing_options end function histogram_get_drawing_options_ptr @ %def histogram_get_drawing_options_ptr @ \subsubsection{Output} <>= subroutine histogram_write (h, unit) type(histogram_t), intent(in) :: h integer, intent(in), optional :: unit integer :: u, i u = given_output_unit (unit); if (u < 0) return call bin_write_header (u) if (allocated (h%bin)) then do i = 1, h%n_bins call bin_write (h%bin(i), u) end do end if write (u, "(A)") write (u, "(A,1x,A)") "#", "Underflow:" call bin_write (h%underflow, u) write (u, "(A)") write (u, "(A,1x,A)") "#", "Overflow:" call bin_write (h%overflow, u) write (u, "(A)") write (u, "(A,1x,A)") "#", "Summary: data within bounds" call observable_write (h%obs_within_bounds, u) write (u, "(A)") write (u, "(A,1x,A)") "#", "Summary: all data" call observable_write (h%obs, u) write (u, "(A)") end subroutine histogram_write @ %def histogram_write @ Write the GAMELAN reader for histogram contents. <>= subroutine histogram_write_gml_reader (h, filename, unit) type(histogram_t), intent(in) :: h type(string_t), intent(in) :: filename integer, intent(in), optional :: unit character(*), parameter :: fmt = "(ES15.8)" integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(2x,A)") 'fromfile "' // char (filename) // '":' write (u, "(4x,A)") 'key "# Histogram:";' write (u, "(4x,A)") 'dx := #' & // real2char (h%width / h%n_bins / 2, fmt) // ';' write (u, "(4x,A)") 'for i withinblock:' write (u, "(6x,A)") 'get x, y, y.d, y.n, y.e;' if (h%drawing_options%with_hbars) then write (u, "(6x,A)") 'plot (' // char (h%drawing_options%dataset) & // ') (x,y) hbar dx;' else write (u, "(6x,A)") 'plot (' // char (h%drawing_options%dataset) & // ') (x,y);' end if if (h%drawing_options%err) then write (u, "(6x,A)") 'plot (' // char (h%drawing_options%dataset) & // '.err) ' & // '(x,y) vbar y.d;' end if !!! Future excess options for plots ! write (u, "(6x,A)") 'if show_excess: ' // & ! & 'plot(dat.e)(x, y plus y.e) hbar dx; fi' write (u, "(4x,A)") 'endfor' write (u, "(2x,A)") 'endfrom' end subroutine histogram_write_gml_reader @ %def histogram_write_gml_reader @ \LaTeX\ and GAMELAN output. <>= subroutine histogram_write_gml_driver (h, filename, unit) type(histogram_t), intent(in) :: h type(string_t), intent(in) :: filename integer, intent(in), optional :: unit type(string_t) :: calc_cmd, bg_cmd, draw_cmd, err_cmd, symb_cmd, fg_cmd integer :: u u = given_output_unit (unit); if (u < 0) return call graph_options_write_tex_header (h%graph_options, unit) write (u, "(2x,A)") char (graph_options_get_gml_setup (h%graph_options)) write (u, "(2x,A)") char (graph_options_get_gml_graphrange & (h%graph_options, x_min=h%lower_bound, x_max=h%upper_bound)) call histogram_write_gml_reader (h, filename, unit) calc_cmd = drawing_options_get_calc_command (h%drawing_options) if (calc_cmd /= "") write (u, "(2x,A)") char (calc_cmd) bg_cmd = drawing_options_get_gml_bg_command (h%drawing_options) if (bg_cmd /= "") write (u, "(2x,A)") char (bg_cmd) draw_cmd = drawing_options_get_draw_command (h%drawing_options) if (draw_cmd /= "") write (u, "(2x,A)") char (draw_cmd) err_cmd = drawing_options_get_err_command (h%drawing_options) if (err_cmd /= "") write (u, "(2x,A)") char (err_cmd) symb_cmd = drawing_options_get_symb_command (h%drawing_options) if (symb_cmd /= "") write (u, "(2x,A)") char (symb_cmd) fg_cmd = drawing_options_get_gml_fg_command (h%drawing_options) if (fg_cmd /= "") write (u, "(2x,A)") char (fg_cmd) write (u, "(2x,A)") char (graph_options_get_gml_x_label (h%graph_options)) write (u, "(2x,A)") char (graph_options_get_gml_y_label (h%graph_options)) call graph_options_write_tex_footer (h%graph_options, unit) write (u, "(A)") "\vspace*{2\baselineskip}" write (u, "(A)") "\begin{flushleft}" write (u, "(A)") "\textbf{Data within bounds:} \\" call observable_write_driver (h%obs_within_bounds, unit, & write_heading=.false.) write (u, "(A)") "\\[0.5\baselineskip]" write (u, "(A)") "\textbf{All data:} \\" call observable_write_driver (h%obs, unit, write_heading=.false.) write (u, "(A)") "\end{flushleft}" end subroutine histogram_write_gml_driver @ %def histogram_write_gml_driver @ Return the header for generic data output as an ifile. <>= subroutine histogram_get_header (h, header, comment) type(histogram_t), intent(in) :: h type(ifile_t), intent(inout) :: header type(string_t), intent(in), optional :: comment type(string_t) :: c if (present (comment)) then c = comment else c = "" end if call ifile_append (header, c // "WHIZARD histogram data") call graph_options_get_header (h%graph_options, header, comment) call ifile_append (header, & c // "range: " // real2string (h%lower_bound) & // " - " // real2string (h%upper_bound)) call ifile_append (header, & c // "counts total: " & // int2char (histogram_get_n_entries_within_bounds (h))) call ifile_append (header, & c // "total average: " & // real2string (histogram_get_average_within_bounds (h)) // " +- " & // real2string (histogram_get_error_within_bounds (h))) end subroutine histogram_get_header @ %def histogram_get_header @ \subsection{Plots} \subsubsection{Points} <>= type :: point_t private real(default) :: x = 0 real(default) :: y = 0 real(default) :: yerr = 0 real(default) :: xerr = 0 type(point_t), pointer :: next => null () end type point_t @ %def point_t <>= interface point_init module procedure point_init_contents module procedure point_init_point end interface <>= subroutine point_init_contents (point, x, y, yerr, xerr) type(point_t), intent(out) :: point real(default), intent(in) :: x, y real(default), intent(in), optional :: yerr, xerr point%x = x point%y = y if (present (yerr)) point%yerr = yerr if (present (xerr)) point%xerr = xerr end subroutine point_init_contents subroutine point_init_point (point, point_in) type(point_t), intent(out) :: point type(point_t), intent(in) :: point_in point%x = point_in%x point%y = point_in%y point%yerr = point_in%yerr point%xerr = point_in%xerr end subroutine point_init_point @ %def point_init <>= function point_get_x (point) result (x) real(default) :: x type(point_t), intent(in) :: point x = point%x end function point_get_x function point_get_y (point) result (y) real(default) :: y type(point_t), intent(in) :: point y = point%y end function point_get_y function point_get_xerr (point) result (xerr) real(default) :: xerr type(point_t), intent(in) :: point xerr = point%xerr end function point_get_xerr function point_get_yerr (point) result (yerr) real(default) :: yerr type(point_t), intent(in) :: point yerr = point%yerr end function point_get_yerr @ %def point_get_x @ %def point_get_y @ %def point_get_xerr @ %def point_get_yerr <>= subroutine point_write_header (unit) integer, intent(in) :: unit character(120) :: buffer integer :: u u = given_output_unit (unit); if (u < 0) return write (buffer, "(A,4(1x," // HISTOGRAM_HEAD_FORMAT // "))") & "#", "x ", "y ", "yerr ", "xerr " write (u, "(A)") trim (buffer) end subroutine point_write_header subroutine point_write (point, unit) type(point_t), intent(in) :: point integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(1x,4(1x," // HISTOGRAM_DATA_FORMAT // "))") & point_get_x (point), & point_get_y (point), & point_get_yerr (point), & point_get_xerr (point) end subroutine point_write @ %def point_write @ \subsubsection{Plots} <>= type :: plot_t private type(point_t), pointer :: first => null () type(point_t), pointer :: last => null () integer :: count = 0 type(graph_options_t) :: graph_options type(drawing_options_t) :: drawing_options end type plot_t @ %def plot_t @ \subsubsection{Initializer/finalizer} Initialize a plot. We provide the lower and upper bound in the $x$ direction. <>= interface plot_init module procedure plot_init_empty module procedure plot_init_plot end interface <>= subroutine plot_init_empty (p, id, graph_options, drawing_options) type(plot_t), intent(out) :: p type(string_t), intent(in) :: id type(graph_options_t), intent(in), optional :: graph_options type(drawing_options_t), intent(in), optional :: drawing_options if (present (graph_options)) then p%graph_options = graph_options else call graph_options_init (p%graph_options) end if call graph_options_set (p%graph_options, id = id) if (present (drawing_options)) then p%drawing_options = drawing_options else call drawing_options_init_plot (p%drawing_options) end if end subroutine plot_init_empty @ %def plot_init @ Initialize a plot by copying another one, optionally merging in a new set of drawing options. Since [[p]] has pointer (sub)components, we have to explicitly deep-copy the original. <>= subroutine plot_init_plot (p, p_in, drawing_options) type(plot_t), intent(out) :: p type(plot_t), intent(in) :: p_in type(drawing_options_t), intent(in), optional :: drawing_options type(point_t), pointer :: current, new current => p_in%first do while (associated (current)) allocate (new) call point_init (new, current) if (associated (p%last)) then p%last%next => new else p%first => new end if p%last => new current => current%next end do p%count = p_in%count p%graph_options = p_in%graph_options if (present (drawing_options)) then p%drawing_options = drawing_options else p%drawing_options = p_in%drawing_options end if end subroutine plot_init_plot @ %def plot_init_plot @ Finalize the plot by deallocating the list of points. <>= subroutine plot_final (plot) type(plot_t), intent(inout) :: plot type(point_t), pointer :: current do while (associated (plot%first)) current => plot%first plot%first => current%next deallocate (current) end do plot%last => null () end subroutine plot_final @ %def plot_final @ \subsubsection{Fill plots} Clear the plot contents, but do not modify the structure. <>= subroutine plot_clear (plot) type(plot_t), intent(inout) :: plot plot%count = 0 call plot_final (plot) end subroutine plot_clear @ %def plot_clear @ Record a value. Successful if the value is within bounds, otherwise it is recorded as under-/overflow. <>= subroutine plot_record_value (plot, x, y, yerr, xerr, success) type(plot_t), intent(inout) :: plot real(default), intent(in) :: x, y real(default), intent(in), optional :: yerr, xerr logical, intent(out), optional :: success type(point_t), pointer :: point plot%count = plot%count + 1 allocate (point) call point_init (point, x, y, yerr, xerr) if (associated (plot%first)) then plot%last%next => point else plot%first => point end if plot%last => point if (present (success)) success = .true. end subroutine plot_record_value @ %def plot_record_value @ \subsubsection{Access contents} The number of points. <>= function plot_get_n_entries (plot) result (n) integer :: n type(plot_t), intent(in) :: plot n = plot%count end function plot_get_n_entries @ %def plot_get_n_entries @ Return a pointer to the graph options. <>= function plot_get_graph_options_ptr (p) result (ptr) type(graph_options_t), pointer :: ptr type(plot_t), intent(in), target :: p ptr => p%graph_options end function plot_get_graph_options_ptr @ %def plot_get_graph_options_ptr @ Return a pointer to the drawing options. <>= function plot_get_drawing_options_ptr (p) result (ptr) type(drawing_options_t), pointer :: ptr type(plot_t), intent(in), target :: p ptr => p%drawing_options end function plot_get_drawing_options_ptr @ %def plot_get_drawing_options_ptr @ \subsubsection{Output} This output format is used by the GAMELAN driver below. <>= subroutine plot_write (plot, unit) type(plot_t), intent(in) :: plot integer, intent(in), optional :: unit type(point_t), pointer :: point integer :: u u = given_output_unit (unit); if (u < 0) return call point_write_header (u) point => plot%first do while (associated (point)) call point_write (point, unit) point => point%next end do write (u, *) write (u, "(A,1x,A)") "#", "Summary:" write (u, "(A,1x,I0)") & "n_entries =", plot_get_n_entries (plot) write (u, *) end subroutine plot_write @ %def plot_write @ Write the GAMELAN reader for plot contents. <>= subroutine plot_write_gml_reader (p, filename, unit) type(plot_t), intent(in) :: p type(string_t), intent(in) :: filename integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, "(2x,A)") 'fromfile "' // char (filename) // '":' write (u, "(4x,A)") 'key "# Plot:";' write (u, "(4x,A)") 'for i withinblock:' write (u, "(6x,A)") 'get x, y, y.err, x.err;' write (u, "(6x,A)") 'plot (' // char (p%drawing_options%dataset) & // ') (x,y);' if (p%drawing_options%err) then write (u, "(6x,A)") 'plot (' // char (p%drawing_options%dataset) & // '.err) (x,y) vbar y.err hbar x.err;' end if write (u, "(4x,A)") 'endfor' write (u, "(2x,A)") 'endfrom' end subroutine plot_write_gml_reader @ %def plot_write_gml_header @ \LaTeX\ and GAMELAN output. Analogous to histogram output. <>= subroutine plot_write_gml_driver (p, filename, unit) type(plot_t), intent(in) :: p type(string_t), intent(in) :: filename integer, intent(in), optional :: unit type(string_t) :: calc_cmd, bg_cmd, draw_cmd, err_cmd, symb_cmd, fg_cmd integer :: u u = given_output_unit (unit); if (u < 0) return call graph_options_write_tex_header (p%graph_options, unit) write (u, "(2x,A)") & char (graph_options_get_gml_setup (p%graph_options)) write (u, "(2x,A)") & char (graph_options_get_gml_graphrange (p%graph_options)) call plot_write_gml_reader (p, filename, unit) calc_cmd = drawing_options_get_calc_command (p%drawing_options) if (calc_cmd /= "") write (u, "(2x,A)") char (calc_cmd) bg_cmd = drawing_options_get_gml_bg_command (p%drawing_options) if (bg_cmd /= "") write (u, "(2x,A)") char (bg_cmd) draw_cmd = drawing_options_get_draw_command (p%drawing_options) if (draw_cmd /= "") write (u, "(2x,A)") char (draw_cmd) err_cmd = drawing_options_get_err_command (p%drawing_options) if (err_cmd /= "") write (u, "(2x,A)") char (err_cmd) symb_cmd = drawing_options_get_symb_command (p%drawing_options) if (symb_cmd /= "") write (u, "(2x,A)") char (symb_cmd) fg_cmd = drawing_options_get_gml_fg_command (p%drawing_options) if (fg_cmd /= "") write (u, "(2x,A)") char (fg_cmd) write (u, "(2x,A)") char (graph_options_get_gml_x_label (p%graph_options)) write (u, "(2x,A)") char (graph_options_get_gml_y_label (p%graph_options)) call graph_options_write_tex_footer (p%graph_options, unit) end subroutine plot_write_gml_driver @ %def plot_write_driver @ Append header for generic data output in ifile format. <>= subroutine plot_get_header (plot, header, comment) type(plot_t), intent(in) :: plot type(ifile_t), intent(inout) :: header type(string_t), intent(in), optional :: comment type(string_t) :: c if (present (comment)) then c = comment else c = "" end if call ifile_append (header, c // "WHIZARD plot data") call graph_options_get_header (plot%graph_options, header, comment) call ifile_append (header, & c // "number of points: " & // int2char (plot_get_n_entries (plot))) end subroutine plot_get_header @ %def plot_get_header @ \subsection{Graphs} A graph is a container for several graph elements. Each graph element is either a plot or a histogram. There is an appropriate base type below (the [[analysis_object_t]]), but to avoid recursion, we define a separate base type here. Note that there is no actual recursion: a graph is an analysis object, but a graph cannot contain graphs. (If we could use type extension, the implementation would be much more transparent.) \subsubsection{Graph elements} Graph elements cannot be filled by the [[record]] command directly. The contents are always copied from elementary histograms or plots. <>= type :: graph_element_t private integer :: type = AN_UNDEFINED type(histogram_t), pointer :: h => null () type(plot_t), pointer :: p => null () end type graph_element_t @ %def graph_element_t <>= subroutine graph_element_final (el) type(graph_element_t), intent(inout) :: el select case (el%type) case (AN_HISTOGRAM) deallocate (el%h) case (AN_PLOT) call plot_final (el%p) deallocate (el%p) end select el%type = AN_UNDEFINED end subroutine graph_element_final @ %def graph_element_final @ Return the number of entries in the graph element: <>= function graph_element_get_n_entries (el) result (n) integer :: n type(graph_element_t), intent(in) :: el select case (el%type) case (AN_HISTOGRAM); n = histogram_get_n_entries (el%h) case (AN_PLOT); n = plot_get_n_entries (el%p) case default; n = 0 end select end function graph_element_get_n_entries @ %def graph_element_get_n_entries @ Return a pointer to the graph / drawing options. <>= function graph_element_get_graph_options_ptr (el) result (ptr) type(graph_options_t), pointer :: ptr type(graph_element_t), intent(in) :: el select case (el%type) case (AN_HISTOGRAM); ptr => histogram_get_graph_options_ptr (el%h) case (AN_PLOT); ptr => plot_get_graph_options_ptr (el%p) case default; ptr => null () end select end function graph_element_get_graph_options_ptr function graph_element_get_drawing_options_ptr (el) result (ptr) type(drawing_options_t), pointer :: ptr type(graph_element_t), intent(in) :: el select case (el%type) case (AN_HISTOGRAM); ptr => histogram_get_drawing_options_ptr (el%h) case (AN_PLOT); ptr => plot_get_drawing_options_ptr (el%p) case default; ptr => null () end select end function graph_element_get_drawing_options_ptr @ %def graph_element_get_graph_options_ptr @ %def graph_element_get_drawing_options_ptr @ Output, simple wrapper for the plot/histogram writer. <>= subroutine graph_element_write (el, unit) type(graph_element_t), intent(in) :: el integer, intent(in), optional :: unit type(graph_options_t), pointer :: gro type(string_t) :: id integer :: u u = given_output_unit (unit); if (u < 0) return gro => graph_element_get_graph_options_ptr (el) id = graph_options_get_id (gro) write (u, "(A,A)") '#', repeat ("-", 78) select case (el%type) case (AN_HISTOGRAM) write (u, "(A)", advance="no") "# Histogram: " write (u, "(1x,A)") char (id) call histogram_write (el%h, unit) case (AN_PLOT) write (u, "(A)", advance="no") "# Plot: " write (u, "(1x,A)") char (id) call plot_write (el%p, unit) end select end subroutine graph_element_write @ %def graph_element_write <>= subroutine graph_element_write_gml_reader (el, filename, unit) type(graph_element_t), intent(in) :: el type(string_t), intent(in) :: filename integer, intent(in), optional :: unit select case (el%type) case (AN_HISTOGRAM); call histogram_write_gml_reader (el%h, filename, unit) case (AN_PLOT); call plot_write_gml_reader (el%p, filename, unit) end select end subroutine graph_element_write_gml_reader @ %def graph_element_write_gml_reader @ \subsubsection{The graph type} The actual graph type contains its own [[graph_options]], which override the individual settings. The [[drawing_options]] are set in the graph elements. This distinction motivates the separation of the two types. <>= type :: graph_t private type(graph_element_t), dimension(:), allocatable :: el type(graph_options_t) :: graph_options end type graph_t @ %def graph_t @ \subsubsection{Initializer/finalizer} The graph is created with a definite number of elements. The elements are filled one by one, optionally with modified drawing options. <>= subroutine graph_init (g, id, n_elements, graph_options) type(graph_t), intent(out) :: g type(string_t), intent(in) :: id integer, intent(in) :: n_elements type(graph_options_t), intent(in), optional :: graph_options allocate (g%el (n_elements)) if (present (graph_options)) then g%graph_options = graph_options else call graph_options_init (g%graph_options) end if call graph_options_set (g%graph_options, id = id) end subroutine graph_init @ %def graph_init <>= subroutine graph_insert_histogram (g, i, h, drawing_options) type(graph_t), intent(inout), target :: g integer, intent(in) :: i type(histogram_t), intent(in) :: h type(drawing_options_t), intent(in), optional :: drawing_options type(graph_options_t), pointer :: gro type(drawing_options_t), pointer :: dro type(string_t) :: id g%el(i)%type = AN_HISTOGRAM allocate (g%el(i)%h) call histogram_init_histogram (g%el(i)%h, h, drawing_options) gro => histogram_get_graph_options_ptr (g%el(i)%h) dro => histogram_get_drawing_options_ptr (g%el(i)%h) id = graph_options_get_id (gro) call drawing_options_set (dro, dataset = "dat." // id) end subroutine graph_insert_histogram @ %def graph_insert_histogram <>= subroutine graph_insert_plot (g, i, p, drawing_options) type(graph_t), intent(inout) :: g integer, intent(in) :: i type(plot_t), intent(in) :: p type(drawing_options_t), intent(in), optional :: drawing_options type(graph_options_t), pointer :: gro type(drawing_options_t), pointer :: dro type(string_t) :: id g%el(i)%type = AN_PLOT allocate (g%el(i)%p) call plot_init_plot (g%el(i)%p, p, drawing_options) gro => plot_get_graph_options_ptr (g%el(i)%p) dro => plot_get_drawing_options_ptr (g%el(i)%p) id = graph_options_get_id (gro) call drawing_options_set (dro, dataset = "dat." // id) end subroutine graph_insert_plot @ %def graph_insert_plot @ Finalizer. <>= subroutine graph_final (g) type(graph_t), intent(inout) :: g integer :: i do i = 1, size (g%el) call graph_element_final (g%el(i)) end do deallocate (g%el) end subroutine graph_final @ %def graph_final @ \subsubsection{Access contents} The number of elements. <>= function graph_get_n_elements (graph) result (n) integer :: n type(graph_t), intent(in) :: graph n = size (graph%el) end function graph_get_n_elements @ %def graph_get_n_elements @ Retrieve a pointer to the drawing options of an element, so they can be modified. (The [[target]] attribute is not actually needed because the components are pointers.) <>= function graph_get_drawing_options_ptr (g, i) result (ptr) type(drawing_options_t), pointer :: ptr type(graph_t), intent(in), target :: g integer, intent(in) :: i ptr => graph_element_get_drawing_options_ptr (g%el(i)) end function graph_get_drawing_options_ptr @ %def graph_get_drawing_options_ptr @ \subsubsection{Output} The default output format just writes histogram and plot data. <>= subroutine graph_write (graph, unit) type(graph_t), intent(in) :: graph integer, intent(in), optional :: unit integer :: i do i = 1, size (graph%el) call graph_element_write (graph%el(i), unit) end do end subroutine graph_write @ %def graph_write @ The GAMELAN driver is not a simple wrapper, but it writes the plot/histogram contents embedded the complete graph. First, data are read in, global background commands next, then individual elements, then global foreground commands. <>= subroutine graph_write_gml_driver (g, filename, unit) type(graph_t), intent(in) :: g type(string_t), intent(in) :: filename type(string_t) :: calc_cmd, bg_cmd, draw_cmd, err_cmd, symb_cmd, fg_cmd integer, intent(in), optional :: unit type(drawing_options_t), pointer :: dro integer :: u, i u = given_output_unit (unit); if (u < 0) return call graph_options_write_tex_header (g%graph_options, unit) write (u, "(2x,A)") & char (graph_options_get_gml_setup (g%graph_options)) write (u, "(2x,A)") & char (graph_options_get_gml_graphrange (g%graph_options)) do i = 1, size (g%el) call graph_element_write_gml_reader (g%el(i), filename, unit) calc_cmd = drawing_options_get_calc_command & (graph_element_get_drawing_options_ptr (g%el(i))) if (calc_cmd /= "") write (u, "(2x,A)") char (calc_cmd) end do bg_cmd = graph_options_get_gml_bg_command (g%graph_options) if (bg_cmd /= "") write (u, "(2x,A)") char (bg_cmd) do i = 1, size (g%el) dro => graph_element_get_drawing_options_ptr (g%el(i)) bg_cmd = drawing_options_get_gml_bg_command (dro) if (bg_cmd /= "") write (u, "(2x,A)") char (bg_cmd) draw_cmd = drawing_options_get_draw_command (dro) if (draw_cmd /= "") write (u, "(2x,A)") char (draw_cmd) err_cmd = drawing_options_get_err_command (dro) if (err_cmd /= "") write (u, "(2x,A)") char (err_cmd) symb_cmd = drawing_options_get_symb_command (dro) if (symb_cmd /= "") write (u, "(2x,A)") char (symb_cmd) fg_cmd = drawing_options_get_gml_fg_command (dro) if (fg_cmd /= "") write (u, "(2x,A)") char (fg_cmd) end do fg_cmd = graph_options_get_gml_fg_command (g%graph_options) if (fg_cmd /= "") write (u, "(2x,A)") char (fg_cmd) write (u, "(2x,A)") char (graph_options_get_gml_x_label (g%graph_options)) write (u, "(2x,A)") char (graph_options_get_gml_y_label (g%graph_options)) call graph_options_write_tex_footer (g%graph_options, unit) end subroutine graph_write_gml_driver @ %def graph_write_gml_driver @ Append header for generic data output in ifile format. <>= subroutine graph_get_header (graph, header, comment) type(graph_t), intent(in) :: graph type(ifile_t), intent(inout) :: header type(string_t), intent(in), optional :: comment type(string_t) :: c if (present (comment)) then c = comment else c = "" end if call ifile_append (header, c // "WHIZARD graph data") call graph_options_get_header (graph%graph_options, header, comment) call ifile_append (header, & c // "number of graph elements: " & // int2char (graph_get_n_elements (graph))) end subroutine graph_get_header @ %def graph_get_header @ \subsection{Analysis objects} This data structure holds all observables, histograms and such that are currently active. We have one global store; individual items are identified by their ID strings. (This should rather be coded by type extension.) <>= integer, parameter :: AN_UNDEFINED = 0 integer, parameter :: AN_OBSERVABLE = 1 integer, parameter :: AN_HISTOGRAM = 2 integer, parameter :: AN_PLOT = 3 integer, parameter :: AN_GRAPH = 4 <>= public :: AN_UNDEFINED, AN_HISTOGRAM, AN_OBSERVABLE, AN_PLOT, AN_GRAPH @ %def AN_UNDEFINED @ %def AN_OBSERVABLE AN_HISTOGRAM AN_PLOT AN_GRAPH <>= type :: analysis_object_t private type(string_t) :: id integer :: type = AN_UNDEFINED type(observable_t), pointer :: obs => null () type(histogram_t), pointer :: h => null () type(plot_t), pointer :: p => null () type(graph_t), pointer :: g => null () type(analysis_object_t), pointer :: next => null () end type analysis_object_t @ %def analysis_object_t @ \subsubsection{Initializer/finalizer} Allocate with the correct type but do not fill initial values. <>= subroutine analysis_object_init (obj, id, type) type(analysis_object_t), intent(out) :: obj type(string_t), intent(in) :: id integer, intent(in) :: type obj%id = id obj%type = type select case (obj%type) case (AN_OBSERVABLE); allocate (obj%obs) case (AN_HISTOGRAM); allocate (obj%h) case (AN_PLOT); allocate (obj%p) case (AN_GRAPH); allocate (obj%g) end select end subroutine analysis_object_init @ %def analysis_object_init <>= subroutine analysis_object_final (obj) type(analysis_object_t), intent(inout) :: obj select case (obj%type) case (AN_OBSERVABLE) deallocate (obj%obs) case (AN_HISTOGRAM) deallocate (obj%h) case (AN_PLOT) call plot_final (obj%p) deallocate (obj%p) case (AN_GRAPH) call graph_final (obj%g) deallocate (obj%g) end select obj%type = AN_UNDEFINED end subroutine analysis_object_final @ %def analysis_object_final @ Clear the analysis object, i.e., reset it to its initial state. Not applicable to graphs, which are always combinations of other existing objects. <>= subroutine analysis_object_clear (obj) type(analysis_object_t), intent(inout) :: obj select case (obj%type) case (AN_OBSERVABLE) call observable_clear (obj%obs) case (AN_HISTOGRAM) call histogram_clear (obj%h) case (AN_PLOT) call plot_clear (obj%p) end select end subroutine analysis_object_clear @ %def analysis_object_clear @ \subsubsection{Fill with data} Record data. The effect depends on the type of analysis object. <>= subroutine analysis_object_record_data (obj, & x, y, yerr, xerr, weight, excess, success) type(analysis_object_t), intent(inout) :: obj real(default), intent(in) :: x real(default), intent(in), optional :: y, yerr, xerr, weight, excess logical, intent(out), optional :: success select case (obj%type) case (AN_OBSERVABLE) if (present (weight)) then call observable_record_value_weighted (obj%obs, x, weight, success) else call observable_record_value_unweighted (obj%obs, x, success) end if case (AN_HISTOGRAM) if (present (weight)) then call histogram_record_value_weighted (obj%h, x, weight, success) else call histogram_record_value_unweighted (obj%h, x, excess, success) end if case (AN_PLOT) if (present (y)) then call plot_record_value (obj%p, x, y, yerr, xerr, success) else if (present (success)) success = .false. end if case default if (present (success)) success = .false. end select end subroutine analysis_object_record_data @ %def analysis_object_record_data @ Explicitly set the pointer to the next object in the list. <>= subroutine analysis_object_set_next_ptr (obj, next) type(analysis_object_t), intent(inout) :: obj type(analysis_object_t), pointer :: next obj%next => next end subroutine analysis_object_set_next_ptr @ %def analysis_object_set_next_ptr @ \subsubsection{Access contents} Return a pointer to the next object in the list. <>= function analysis_object_get_next_ptr (obj) result (next) type(analysis_object_t), pointer :: next type(analysis_object_t), intent(in) :: obj next => obj%next end function analysis_object_get_next_ptr @ %def analysis_object_get_next_ptr @ Return data as appropriate for the object type. <>= function analysis_object_get_n_elements (obj) result (n) integer :: n type(analysis_object_t), intent(in) :: obj select case (obj%type) case (AN_HISTOGRAM) n = 1 case (AN_PLOT) n = 1 case (AN_GRAPH) n = graph_get_n_elements (obj%g) case default n = 0 end select end function analysis_object_get_n_elements function analysis_object_get_n_entries (obj, within_bounds) result (n) integer :: n type(analysis_object_t), intent(in) :: obj logical, intent(in), optional :: within_bounds logical :: wb select case (obj%type) case (AN_OBSERVABLE) n = observable_get_n_entries (obj%obs) case (AN_HISTOGRAM) wb = .false.; if (present (within_bounds)) wb = within_bounds if (wb) then n = histogram_get_n_entries_within_bounds (obj%h) else n = histogram_get_n_entries (obj%h) end if case (AN_PLOT) n = plot_get_n_entries (obj%p) case default n = 0 end select end function analysis_object_get_n_entries function analysis_object_get_average (obj, within_bounds) result (avg) real(default) :: avg type(analysis_object_t), intent(in) :: obj logical, intent(in), optional :: within_bounds logical :: wb select case (obj%type) case (AN_OBSERVABLE) avg = observable_get_average (obj%obs) case (AN_HISTOGRAM) wb = .false.; if (present (within_bounds)) wb = within_bounds if (wb) then avg = histogram_get_average_within_bounds (obj%h) else avg = histogram_get_average (obj%h) end if case default avg = 0 end select end function analysis_object_get_average function analysis_object_get_error (obj, within_bounds) result (err) real(default) :: err type(analysis_object_t), intent(in) :: obj logical, intent(in), optional :: within_bounds logical :: wb select case (obj%type) case (AN_OBSERVABLE) err = observable_get_error (obj%obs) case (AN_HISTOGRAM) wb = .false.; if (present (within_bounds)) wb = within_bounds if (wb) then err = histogram_get_error_within_bounds (obj%h) else err = histogram_get_error (obj%h) end if case default err = 0 end select end function analysis_object_get_error @ %def analysis_object_get_n_elements @ %def analysis_object_get_n_entries @ %def analysis_object_get_average @ %def analysis_object_get_error @ Return pointers to the actual contents: <>= function analysis_object_get_observable_ptr (obj) result (obs) type(observable_t), pointer :: obs type(analysis_object_t), intent(in) :: obj select case (obj%type) case (AN_OBSERVABLE); obs => obj%obs case default; obs => null () end select end function analysis_object_get_observable_ptr function analysis_object_get_histogram_ptr (obj) result (h) type(histogram_t), pointer :: h type(analysis_object_t), intent(in) :: obj select case (obj%type) case (AN_HISTOGRAM); h => obj%h case default; h => null () end select end function analysis_object_get_histogram_ptr function analysis_object_get_plot_ptr (obj) result (plot) type(plot_t), pointer :: plot type(analysis_object_t), intent(in) :: obj select case (obj%type) case (AN_PLOT); plot => obj%p case default; plot => null () end select end function analysis_object_get_plot_ptr function analysis_object_get_graph_ptr (obj) result (g) type(graph_t), pointer :: g type(analysis_object_t), intent(in) :: obj select case (obj%type) case (AN_GRAPH); g => obj%g case default; g => null () end select end function analysis_object_get_graph_ptr @ %def analysis_object_get_observable_ptr @ %def analysis_object_get_histogram_ptr @ %def analysis_object_get_plot_ptr @ %def analysis_object_get_graph_ptr @ Return true if the object has a graphical representation: <>= function analysis_object_has_plot (obj) result (flag) logical :: flag type(analysis_object_t), intent(in) :: obj select case (obj%type) case (AN_HISTOGRAM); flag = .true. case (AN_PLOT); flag = .true. case (AN_GRAPH); flag = .true. case default; flag = .false. end select end function analysis_object_has_plot @ %def analysis_object_has_plot @ \subsubsection{Output} <>= subroutine analysis_object_write (obj, unit, verbose) type(analysis_object_t), intent(in) :: obj integer, intent(in), optional :: unit logical, intent(in), optional :: verbose logical :: verb integer :: u u = given_output_unit (unit); if (u < 0) return verb = .false.; if (present (verbose)) verb = verbose write (u, "(A)") repeat ("#", 79) select case (obj%type) case (AN_OBSERVABLE) write (u, "(A)", advance="no") "# Observable:" case (AN_HISTOGRAM) write (u, "(A)", advance="no") "# Histogram: " case (AN_PLOT) write (u, "(A)", advance="no") "# Plot: " case (AN_GRAPH) write (u, "(A)", advance="no") "# Graph: " case default write (u, "(A)") "# [undefined analysis object]" return end select write (u, "(1x,A)") char (obj%id) select case (obj%type) case (AN_OBSERVABLE) call observable_write (obj%obs, unit) case (AN_HISTOGRAM) if (verb) then call graph_options_write (obj%h%graph_options, unit) write (u, *) call drawing_options_write (obj%h%drawing_options, unit) write (u, *) end if call histogram_write (obj%h, unit) case (AN_PLOT) if (verb) then call graph_options_write (obj%p%graph_options, unit) write (u, *) call drawing_options_write (obj%p%drawing_options, unit) write (u, *) end if call plot_write (obj%p, unit) case (AN_GRAPH) call graph_write (obj%g, unit) end select end subroutine analysis_object_write @ %def analysis_object_write @ Write the object part of the \LaTeX\ driver file. <>= subroutine analysis_object_write_driver (obj, filename, unit) type(analysis_object_t), intent(in) :: obj type(string_t), intent(in) :: filename integer, intent(in), optional :: unit select case (obj%type) case (AN_OBSERVABLE) call observable_write_driver (obj%obs, unit) case (AN_HISTOGRAM) call histogram_write_gml_driver (obj%h, filename, unit) case (AN_PLOT) call plot_write_gml_driver (obj%p, filename, unit) case (AN_GRAPH) call graph_write_gml_driver (obj%g, filename, unit) end select end subroutine analysis_object_write_driver @ %def analysis_object_write_driver @ Return a data header for external formats, in ifile form. <>= subroutine analysis_object_get_header (obj, header, comment) type(analysis_object_t), intent(in) :: obj type(ifile_t), intent(inout) :: header type(string_t), intent(in), optional :: comment select case (obj%type) case (AN_HISTOGRAM) call histogram_get_header (obj%h, header, comment) case (AN_PLOT) call plot_get_header (obj%p, header, comment) end select end subroutine analysis_object_get_header @ %def analysis_object_get_header @ \subsection{Analysis object iterator} Analysis objects are containers which have iterable data structures: histograms/bins and plots/points. If they are to be treated on a common basis, it is useful to have an iterator which hides the implementation details. The iterator is used only for elementary analysis objects that contain plot data: histograms or plots. It is invalid for meta-objects (graphs) and non-graphical objects (observables). <>= public :: analysis_iterator_t <>= type :: analysis_iterator_t private integer :: type = AN_UNDEFINED type(analysis_object_t), pointer :: object => null () integer :: index = 1 type(point_t), pointer :: point => null () end type @ %def analysis_iterator_t @ The initializer places the iterator at the beginning of the analysis object. <>= subroutine analysis_iterator_init (iterator, object) type(analysis_iterator_t), intent(out) :: iterator type(analysis_object_t), intent(in), target :: object iterator%object => object if (associated (iterator%object)) then iterator%type = iterator%object%type select case (iterator%type) case (AN_PLOT) iterator%point => iterator%object%p%first end select end if end subroutine analysis_iterator_init @ %def analysis_iterator_init @ The iterator is valid as long as it points to an existing entry. An iterator for a data object without array data (observable) is always invalid. <>= public :: analysis_iterator_is_valid <>= function analysis_iterator_is_valid (iterator) result (valid) logical :: valid type(analysis_iterator_t), intent(in) :: iterator if (associated (iterator%object)) then select case (iterator%type) case (AN_HISTOGRAM) valid = iterator%index <= histogram_get_n_bins (iterator%object%h) case (AN_PLOT) valid = associated (iterator%point) case default valid = .false. end select else valid = .false. end if end function analysis_iterator_is_valid @ %def analysis_iterator_is_valid @ Advance the iterator. <>= public :: analysis_iterator_advance <>= subroutine analysis_iterator_advance (iterator) type(analysis_iterator_t), intent(inout) :: iterator if (associated (iterator%object)) then select case (iterator%type) case (AN_PLOT) iterator%point => iterator%point%next end select iterator%index = iterator%index + 1 end if end subroutine analysis_iterator_advance @ %def analysis_iterator_advance @ Retrieve the object type: <>= public :: analysis_iterator_get_type <>= function analysis_iterator_get_type (iterator) result (type) integer :: type type(analysis_iterator_t), intent(in) :: iterator type = iterator%type end function analysis_iterator_get_type @ %def analysis_iterator_get_type @ Use the iterator to retrieve data. We implement a common routine which takes the data descriptors as optional arguments. Data which do not occur in the selected type trigger to an error condition. The iterator must point to a valid entry. <>= public :: analysis_iterator_get_data <>= subroutine analysis_iterator_get_data (iterator, & x, y, yerr, xerr, width, excess, index, n_total) type(analysis_iterator_t), intent(in) :: iterator real(default), intent(out), optional :: x, y, yerr, xerr, width, excess integer, intent(out), optional :: index, n_total select case (iterator%type) case (AN_HISTOGRAM) if (present (x)) & x = bin_get_midpoint (iterator%object%h%bin(iterator%index)) if (present (y)) & y = bin_get_sum (iterator%object%h%bin(iterator%index)) if (present (yerr)) & yerr = bin_get_error (iterator%object%h%bin(iterator%index)) if (present (xerr)) & call invalid ("histogram", "xerr") if (present (width)) & width = bin_get_width (iterator%object%h%bin(iterator%index)) if (present (excess)) & excess = bin_get_excess (iterator%object%h%bin(iterator%index)) if (present (index)) & index = iterator%index if (present (n_total)) & n_total = histogram_get_n_bins (iterator%object%h) case (AN_PLOT) if (present (x)) & x = point_get_x (iterator%point) if (present (y)) & y = point_get_y (iterator%point) if (present (yerr)) & yerr = point_get_yerr (iterator%point) if (present (xerr)) & xerr = point_get_xerr (iterator%point) if (present (width)) & call invalid ("plot", "width") if (present (excess)) & call invalid ("plot", "excess") if (present (index)) & index = iterator%index if (present (n_total)) & n_total = plot_get_n_entries (iterator%object%p) case default call msg_bug ("analysis_iterator_get_data: called " & // "for unsupported analysis object type") end select contains subroutine invalid (typestr, objstr) character(*), intent(in) :: typestr, objstr call msg_bug ("analysis_iterator_get_data: attempt to get '" & // objstr // "' for type '" // typestr // "'") end subroutine invalid end subroutine analysis_iterator_get_data @ %def analysis_iterator_get_data @ \subsection{Analysis store} This data structure holds all observables, histograms and such that are currently active. We have one global store; individual items are identified by their ID strings and types. <>= type(analysis_store_t), save :: analysis_store @ %def analysis_store <>= type :: analysis_store_t private type(analysis_object_t), pointer :: first => null () type(analysis_object_t), pointer :: last => null () end type analysis_store_t @ %def analysis_store_t @ Delete the analysis store <>= public :: analysis_final <>= subroutine analysis_final () type(analysis_object_t), pointer :: current do while (associated (analysis_store%first)) current => analysis_store%first analysis_store%first => current%next call analysis_object_final (current) end do analysis_store%last => null () end subroutine analysis_final @ %def analysis_final @ Append a new analysis object <>= subroutine analysis_store_append_object (id, type) type(string_t), intent(in) :: id integer, intent(in) :: type type(analysis_object_t), pointer :: obj allocate (obj) call analysis_object_init (obj, id, type) if (associated (analysis_store%last)) then analysis_store%last%next => obj else analysis_store%first => obj end if analysis_store%last => obj end subroutine analysis_store_append_object @ %def analysis_store_append_object @ Return a pointer to the analysis object with given ID. <>= function analysis_store_get_object_ptr (id) result (obj) type(string_t), intent(in) :: id type(analysis_object_t), pointer :: obj obj => analysis_store%first do while (associated (obj)) if (obj%id == id) return obj => obj%next end do end function analysis_store_get_object_ptr @ %def analysis_store_get_object_ptr @ Initialize an analysis object: either reset it if present, or append a new entry. <>= subroutine analysis_store_init_object (id, type, obj) type(string_t), intent(in) :: id integer, intent(in) :: type type(analysis_object_t), pointer :: obj, next obj => analysis_store_get_object_ptr (id) if (associated (obj)) then next => analysis_object_get_next_ptr (obj) call analysis_object_final (obj) call analysis_object_init (obj, id, type) call analysis_object_set_next_ptr (obj, next) else call analysis_store_append_object (id, type) obj => analysis_store%last end if end subroutine analysis_store_init_object @ %def analysis_store_init_object @ Get the type of a analysis object <>= public :: analysis_store_get_object_type <>= function analysis_store_get_object_type (id) result (type) type(string_t), intent(in) :: id integer :: type type(analysis_object_t), pointer :: object object => analysis_store_get_object_ptr (id) if (associated (object)) then type = object%type else type = AN_UNDEFINED end if end function analysis_store_get_object_type @ %def analysis_store_get_object_type @ Return the number of objects in the store. <>= function analysis_store_get_n_objects () result (n) integer :: n type(analysis_object_t), pointer :: current n = 0 current => analysis_store%first do while (associated (current)) n = n + 1 current => current%next end do end function analysis_store_get_n_objects @ %def analysis_store_get_n_objects @ Allocate an array and fill it with all existing IDs. <>= public :: analysis_store_get_ids <>= subroutine analysis_store_get_ids (id) type(string_t), dimension(:), allocatable, intent(out) :: id type(analysis_object_t), pointer :: current integer :: i allocate (id (analysis_store_get_n_objects())) i = 0 current => analysis_store%first do while (associated (current)) i = i + 1 id(i) = current%id current => current%next end do end subroutine analysis_store_get_ids @ %def analysis_store_get_ids @ \subsection{\LaTeX\ driver file} Write a driver file for all objects in the store. <>= subroutine analysis_store_write_driver_all (filename_data, unit) type(string_t), intent(in) :: filename_data integer, intent(in), optional :: unit type(analysis_object_t), pointer :: obj call analysis_store_write_driver_header (unit) obj => analysis_store%first do while (associated (obj)) call analysis_object_write_driver (obj, filename_data, unit) obj => obj%next end do call analysis_store_write_driver_footer (unit) end subroutine analysis_store_write_driver_all @ %def analysis_store_write_driver_all @ Write a driver file for an array of objects. <>= subroutine analysis_store_write_driver_obj (filename_data, id, unit) type(string_t), intent(in) :: filename_data type(string_t), dimension(:), intent(in) :: id integer, intent(in), optional :: unit type(analysis_object_t), pointer :: obj integer :: i call analysis_store_write_driver_header (unit) do i = 1, size (id) obj => analysis_store_get_object_ptr (id(i)) if (associated (obj)) & call analysis_object_write_driver (obj, filename_data, unit) end do call analysis_store_write_driver_footer (unit) end subroutine analysis_store_write_driver_obj @ %def analysis_store_write_driver_obj @ The beginning of the driver file. <>= subroutine analysis_store_write_driver_header (unit) integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write (u, '(A)') "\documentclass[12pt]{article}" write (u, *) write (u, '(A)') "\usepackage{gamelan}" write (u, '(A)') "\usepackage{amsmath}" write (u, '(A)') "\usepackage{ifpdf}" write (u, '(A)') "\ifpdf" write (u, '(A)') " \DeclareGraphicsRule{*}{mps}{*}{}" write (u, '(A)') "\else" write (u, '(A)') " \DeclareGraphicsRule{*}{eps}{*}{}" write (u, '(A)') "\fi" write (u, *) write (u, '(A)') "\begin{document}" write (u, '(A)') "\begin{gmlfile}" write (u, *) write (u, '(A)') "\begin{gmlcode}" write (u, '(A)') " color col.default, col.excess;" write (u, '(A)') " col.default = 0.9white;" write (u, '(A)') " col.excess = red;" write (u, '(A)') " boolean show_excess;" !!! Future excess options for plots ! if (mcs(1)%plot_excess .and. mcs(1)%unweighted) then ! write (u, '(A)') " show_excess = true;" ! else write (u, '(A)') " show_excess = false;" ! end if write (u, '(A)') "\end{gmlcode}" write (u, *) end subroutine analysis_store_write_driver_header @ %def analysis_store_write_driver_header @ The end of the driver file. <>= subroutine analysis_store_write_driver_footer (unit) integer, intent(in), optional :: unit integer :: u u = given_output_unit (unit); if (u < 0) return write(u, *) write(u, '(A)') "\end{gmlfile}" write(u, '(A)') "\end{document}" end subroutine analysis_store_write_driver_footer @ %def analysis_store_write_driver_footer @ \subsection{API} \subsubsection{Creating new objects} The specific versions below: <>= public :: analysis_init_observable <>= subroutine analysis_init_observable (id, obs_label, obs_unit, graph_options) type(string_t), intent(in) :: id type(string_t), intent(in), optional :: obs_label, obs_unit type(graph_options_t), intent(in), optional :: graph_options type(analysis_object_t), pointer :: obj type(observable_t), pointer :: obs call analysis_store_init_object (id, AN_OBSERVABLE, obj) obs => analysis_object_get_observable_ptr (obj) call observable_init (obs, obs_label, obs_unit, graph_options) end subroutine analysis_init_observable @ %def analysis_init_observable <>= public :: analysis_init_histogram <>= interface analysis_init_histogram module procedure analysis_init_histogram_n_bins module procedure analysis_init_histogram_bin_width end interface <>= subroutine analysis_init_histogram_n_bins & (id, lower_bound, upper_bound, n_bins, normalize_bins, & obs_label, obs_unit, graph_options, drawing_options) type(string_t), intent(in) :: id real(default), intent(in) :: lower_bound, upper_bound integer, intent(in) :: n_bins logical, intent(in) :: normalize_bins type(string_t), intent(in), optional :: obs_label, obs_unit type(graph_options_t), intent(in), optional :: graph_options type(drawing_options_t), intent(in), optional :: drawing_options type(analysis_object_t), pointer :: obj type(histogram_t), pointer :: h call analysis_store_init_object (id, AN_HISTOGRAM, obj) h => analysis_object_get_histogram_ptr (obj) call histogram_init (h, id, & lower_bound, upper_bound, n_bins, normalize_bins, & obs_label, obs_unit, graph_options, drawing_options) end subroutine analysis_init_histogram_n_bins subroutine analysis_init_histogram_bin_width & (id, lower_bound, upper_bound, bin_width, normalize_bins, & obs_label, obs_unit, graph_options, drawing_options) type(string_t), intent(in) :: id real(default), intent(in) :: lower_bound, upper_bound, bin_width logical, intent(in) :: normalize_bins type(string_t), intent(in), optional :: obs_label, obs_unit type(graph_options_t), intent(in), optional :: graph_options type(drawing_options_t), intent(in), optional :: drawing_options type(analysis_object_t), pointer :: obj type(histogram_t), pointer :: h call analysis_store_init_object (id, AN_HISTOGRAM, obj) h => analysis_object_get_histogram_ptr (obj) call histogram_init (h, id, & lower_bound, upper_bound, bin_width, normalize_bins, & obs_label, obs_unit, graph_options, drawing_options) end subroutine analysis_init_histogram_bin_width @ %def analysis_init_histogram_n_bins @ %def analysis_init_histogram_bin_width <>= public :: analysis_init_plot <>= subroutine analysis_init_plot (id, graph_options, drawing_options) type(string_t), intent(in) :: id type(graph_options_t), intent(in), optional :: graph_options type(drawing_options_t), intent(in), optional :: drawing_options type(analysis_object_t), pointer :: obj type(plot_t), pointer :: plot call analysis_store_init_object (id, AN_PLOT, obj) plot => analysis_object_get_plot_ptr (obj) call plot_init (plot, id, graph_options, drawing_options) end subroutine analysis_init_plot @ %def analysis_init_plot <>= public :: analysis_init_graph <>= subroutine analysis_init_graph (id, n_elements, graph_options) type(string_t), intent(in) :: id integer, intent(in) :: n_elements type(graph_options_t), intent(in), optional :: graph_options type(analysis_object_t), pointer :: obj type(graph_t), pointer :: graph call analysis_store_init_object (id, AN_GRAPH, obj) graph => analysis_object_get_graph_ptr (obj) call graph_init (graph, id, n_elements, graph_options) end subroutine analysis_init_graph @ %def analysis_init_graph @ \subsubsection{Recording data} This procedure resets an object or the whole store to its initial state. <>= public :: analysis_clear <>= interface analysis_clear module procedure analysis_store_clear_obj module procedure analysis_store_clear_all end interface <>= subroutine analysis_store_clear_obj (id) type(string_t), intent(in) :: id type(analysis_object_t), pointer :: obj obj => analysis_store_get_object_ptr (id) if (associated (obj)) then call analysis_object_clear (obj) end if end subroutine analysis_store_clear_obj subroutine analysis_store_clear_all () type(analysis_object_t), pointer :: obj obj => analysis_store%first do while (associated (obj)) call analysis_object_clear (obj) obj => obj%next end do end subroutine analysis_store_clear_all @ %def analysis_clear @ There is one generic recording function whose behavior depends on the type of analysis object. <>= public :: analysis_record_data <>= subroutine analysis_record_data (id, x, y, yerr, xerr, & weight, excess, success, exist) type(string_t), intent(in) :: id real(default), intent(in) :: x real(default), intent(in), optional :: y, yerr, xerr, weight, excess logical, intent(out), optional :: success, exist type(analysis_object_t), pointer :: obj obj => analysis_store_get_object_ptr (id) if (associated (obj)) then call analysis_object_record_data (obj, x, y, yerr, xerr, & weight, excess, success) if (present (exist)) exist = .true. else if (present (success)) success = .false. if (present (exist)) exist = .false. end if end subroutine analysis_record_data @ %def analysis_record_data @ \subsubsection{Build a graph} This routine sets up the array of graph elements by copying the graph elements given as input. The object must exist and already be initialized as a graph. <>= public :: analysis_fill_graph <>= subroutine analysis_fill_graph (id, i, id_in, drawing_options) type(string_t), intent(in) :: id integer, intent(in) :: i type(string_t), intent(in) :: id_in type(drawing_options_t), intent(in), optional :: drawing_options type(analysis_object_t), pointer :: obj type(graph_t), pointer :: g type(histogram_t), pointer :: h type(plot_t), pointer :: p obj => analysis_store_get_object_ptr (id) g => analysis_object_get_graph_ptr (obj) obj => analysis_store_get_object_ptr (id_in) if (associated (obj)) then select case (obj%type) case (AN_HISTOGRAM) h => analysis_object_get_histogram_ptr (obj) call graph_insert_histogram (g, i, h, drawing_options) case (AN_PLOT) p => analysis_object_get_plot_ptr (obj) call graph_insert_plot (g, i, p, drawing_options) case default call msg_error ("Graph '" // char (id) // "': Element '" & // char (id_in) // "' is neither histogram nor plot.") end select else call msg_error ("Graph '" // char (id) // "': Element '" & // char (id_in) // "' is undefined.") end if end subroutine analysis_fill_graph @ %def analysis_fill_graph @ \subsubsection{Retrieve generic results} Check if a named object exists. <>= public :: analysis_exists <>= function analysis_exists (id) result (flag) type(string_t), intent(in) :: id logical :: flag type(analysis_object_t), pointer :: obj flag = .true. obj => analysis_store%first do while (associated (obj)) if (obj%id == id) return obj => obj%next end do flag = .false. end function analysis_exists @ %def analysis_exists @ The following functions should work for all kinds of analysis object: <>= public :: analysis_get_n_elements public :: analysis_get_n_entries public :: analysis_get_average public :: analysis_get_error <>= function analysis_get_n_elements (id) result (n) integer :: n type(string_t), intent(in) :: id type(analysis_object_t), pointer :: obj obj => analysis_store_get_object_ptr (id) if (associated (obj)) then n = analysis_object_get_n_elements (obj) else n = 0 end if end function analysis_get_n_elements function analysis_get_n_entries (id, within_bounds) result (n) integer :: n type(string_t), intent(in) :: id logical, intent(in), optional :: within_bounds type(analysis_object_t), pointer :: obj obj => analysis_store_get_object_ptr (id) if (associated (obj)) then n = analysis_object_get_n_entries (obj, within_bounds) else n = 0 end if end function analysis_get_n_entries function analysis_get_average (id, within_bounds) result (avg) real(default) :: avg type(string_t), intent(in) :: id type(analysis_object_t), pointer :: obj logical, intent(in), optional :: within_bounds obj => analysis_store_get_object_ptr (id) if (associated (obj)) then avg = analysis_object_get_average (obj, within_bounds) else avg = 0 end if end function analysis_get_average function analysis_get_error (id, within_bounds) result (err) real(default) :: err type(string_t), intent(in) :: id type(analysis_object_t), pointer :: obj logical, intent(in), optional :: within_bounds obj => analysis_store_get_object_ptr (id) if (associated (obj)) then err = analysis_object_get_error (obj, within_bounds) else err = 0 end if end function analysis_get_error @ %def analysis_get_n_elements @ %def analysis_get_n_entries @ %def analysis_get_average @ %def analysis_get_error @ Return true if any analysis object is graphical <>= public :: analysis_has_plots <>= interface analysis_has_plots module procedure analysis_has_plots_any module procedure analysis_has_plots_obj end interface <>= function analysis_has_plots_any () result (flag) logical :: flag type(analysis_object_t), pointer :: obj flag = .false. obj => analysis_store%first do while (associated (obj)) flag = analysis_object_has_plot (obj) if (flag) return end do end function analysis_has_plots_any function analysis_has_plots_obj (id) result (flag) logical :: flag type(string_t), dimension(:), intent(in) :: id type(analysis_object_t), pointer :: obj integer :: i flag = .false. do i = 1, size (id) obj => analysis_store_get_object_ptr (id(i)) if (associated (obj)) then flag = analysis_object_has_plot (obj) if (flag) return end if end do end function analysis_has_plots_obj @ %def analysis_has_plots @ \subsubsection{Iterators} Initialize an iterator for the given object. If the object does not exist or has wrong type, the iterator will be invalid. <>= public :: analysis_init_iterator <>= subroutine analysis_init_iterator (id, iterator) type(string_t), intent(in) :: id type(analysis_iterator_t), intent(out) :: iterator type(analysis_object_t), pointer :: obj obj => analysis_store_get_object_ptr (id) if (associated (obj)) call analysis_iterator_init (iterator, obj) end subroutine analysis_init_iterator @ %def analysis_init_iterator @ \subsubsection{Output} <>= public :: analysis_write <>= interface analysis_write module procedure analysis_write_object module procedure analysis_write_all end interface @ %def interface <>= subroutine analysis_write_object (id, unit, verbose) type(string_t), intent(in) :: id integer, intent(in), optional :: unit logical, intent(in), optional :: verbose type(analysis_object_t), pointer :: obj obj => analysis_store_get_object_ptr (id) if (associated (obj)) then call analysis_object_write (obj, unit, verbose) else call msg_error ("Analysis object '" // char (id) // "' not found") end if end subroutine analysis_write_object subroutine analysis_write_all (unit, verbose) integer, intent(in), optional :: unit logical, intent(in), optional :: verbose type(analysis_object_t), pointer :: obj integer :: u u = given_output_unit (unit); if (u < 0) return obj => analysis_store%first do while (associated (obj)) call analysis_object_write (obj, unit, verbose) obj => obj%next end do end subroutine analysis_write_all @ %def analysis_write_object @ %def analysis_write_all <>= public :: analysis_write_driver <>= subroutine analysis_write_driver (filename_data, id, unit) type(string_t), intent(in) :: filename_data type(string_t), dimension(:), intent(in), optional :: id integer, intent(in), optional :: unit if (present (id)) then call analysis_store_write_driver_obj (filename_data, id, unit) else call analysis_store_write_driver_all (filename_data, unit) end if end subroutine analysis_write_driver @ %def analysis_write_driver <>= public :: analysis_compile_tex <>= subroutine analysis_compile_tex (file, has_gmlcode, os_data) type(string_t), intent(in) :: file logical, intent(in) :: has_gmlcode type(os_data_t), intent(in) :: os_data integer :: status if (os_data%event_analysis_ps) then call os_system_call ("make compile " // os_data%makeflags // " -f " // & char (file) // "_ana.makefile", status) if (status /= 0) then call msg_error ("Unable to compile analysis output file") end if else call msg_warning ("Skipping results display because " & // "latex/mpost/dvips is not available") end if end subroutine analysis_compile_tex @ %def analysis_compile_tex @ Write header for generic data output to an ifile. <>= public :: analysis_get_header <>= subroutine analysis_get_header (id, header, comment) type(string_t), intent(in) :: id type(ifile_t), intent(inout) :: header type(string_t), intent(in), optional :: comment type(analysis_object_t), pointer :: object object => analysis_store_get_object_ptr (id) if (associated (object)) then call analysis_object_get_header (object, header, comment) end if end subroutine analysis_get_header @ %def analysis_get_header @ Write a makefile in order to do the compile steps. <>= public :: analysis_write_makefile <>= subroutine analysis_write_makefile (filename, unit, has_gmlcode, os_data) type(string_t), intent(in) :: filename integer, intent(in) :: unit logical, intent(in) :: has_gmlcode type(os_data_t), intent(in) :: os_data write (unit, "(3A)") "# WHIZARD: Makefile for analysis '", & char (filename), "'" write (unit, "(A)") "# Automatically generated file, do not edit" write (unit, "(A)") "" write (unit, "(A)") "# LaTeX setup" write (unit, "(A)") "LATEX = " // char (os_data%latex) write (unit, "(A)") "MPOST = " // char (os_data%mpost) write (unit, "(A)") "GML = " // char (os_data%gml) write (unit, "(A)") "DVIPS = " // char (os_data%dvips) write (unit, "(A)") "PS2PDF = " // char (os_data%ps2pdf) write (unit, "(A)") 'TEX_FLAGS = "$$TEXINPUTS:' // & char(os_data%whizard_texpath) // '"' write (unit, "(A)") 'MP_FLAGS = "$$MPINPUTS:' // & char(os_data%whizard_texpath) // '"' write (unit, "(A)") "" write (unit, "(5A)") "TEX_SOURCES = ", char (filename), ".tex" if (os_data%event_analysis_pdf) then write (unit, "(5A)") "TEX_OBJECTS = ", char (filename), ".pdf" else write (unit, "(5A)") "TEX_OBJECTS = ", char (filename), ".ps" end if if (os_data%event_analysis_ps) then if (os_data%event_analysis_pdf) then write (unit, "(5A)") char (filename), ".pdf: ", & char (filename), ".tex" else write (unit, "(5A)") char (filename), ".ps: ", & char (filename), ".tex" end if write (unit, "(5A)") TAB, "-TEXINPUTS=$(TEX_FLAGS) $(LATEX) " // & char (filename) // ".tex" if (has_gmlcode) then write (unit, "(5A)") TAB, "$(GML) " // char (filename) write (unit, "(5A)") TAB, "TEXINPUTS=$(TEX_FLAGS) $(LATEX) " // & char (filename) // ".tex" end if write (unit, "(5A)") TAB, "$(DVIPS) -o " // char (filename) // ".ps " // & char (filename) // ".dvi" if (os_data%event_analysis_pdf) then write (unit, "(5A)") TAB, "$(PS2PDF) " // char (filename) // ".ps" end if end if write (unit, "(A)") write (unit, "(A)") "compile: $(TEX_OBJECTS)" write (unit, "(A)") ".PHONY: compile" write (unit, "(A)") write (unit, "(5A)") "CLEAN_OBJECTS = ", char (filename), ".aux" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".log" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".dvi" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".out" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".[1-9]" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".[1-9][0-9]" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".[1-9][0-9][0-9]" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".t[1-9]" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".t[1-9][0-9]" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".t[1-9][0-9][0-9]" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".ltp" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".mp" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".mpx" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".dvi" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".ps" write (unit, "(5A)") "CLEAN_OBJECTS += ", char (filename), ".pdf" write (unit, "(A)") write (unit, "(A)") "# Generic cleanup targets" write (unit, "(A)") "clean-objects:" write (unit, "(A)") TAB // "rm -f $(CLEAN_OBJECTS)" write (unit, "(A)") "" write (unit, "(A)") "clean: clean-objects" write (unit, "(A)") ".PHONY: clean" end subroutine analysis_write_makefile @ %def analysis_write_makefile @ \subsection{Unit tests} Test module, followed by the corresponding implementation module. <<[[analysis_ut.f90]]>>= <> module analysis_ut use unit_tests use analysis_uti <> <> contains <> end module analysis_ut @ %def analysis_ut @ <<[[analysis_uti.f90]]>>= <> module analysis_uti <> <> use format_defs, only: FMT_19 use analysis <> <> contains <> end module analysis_uti @ %def analysis_ut @ API: driver for the unit tests below. <>= public :: analysis_test <>= subroutine analysis_test (u, results) integer, intent(in) :: u type(test_results_t), intent(inout) :: results <> end subroutine analysis_test @ %def analysis_test <>= call test (analysis_1, "analysis_1", & "check elementary analysis building blocks", & u, results) <>= public :: analysis_1 <>= subroutine analysis_1 (u) integer, intent(in) :: u type(string_t) :: id1, id2, id3, id4 integer :: i id1 = "foo" id2 = "bar" id3 = "hist" id4 = "plot" write (u, "(A)") "* Test output: Analysis" write (u, "(A)") "* Purpose: test the analysis routines" write (u, "(A)") call analysis_init_observable (id1) call analysis_init_observable (id2) call analysis_init_histogram & (id3, 0.5_default, 5.5_default, 1._default, normalize_bins=.false.) call analysis_init_plot (id4) do i = 1, 3 write (u, "(A,1x," // FMT_19 // ")") "data = ", real(i,default) call analysis_record_data (id1, real(i,default)) call analysis_record_data (id2, real(i,default), & weight=real(i,default)) call analysis_record_data (id3, real(i,default)) call analysis_record_data (id4, real(i,default), real(i,default)**2) end do write (u, "(A,10(1x,I5))") "n_entries = ", & analysis_get_n_entries (id1), & analysis_get_n_entries (id2), & analysis_get_n_entries (id3), & analysis_get_n_entries (id3, within_bounds = .true.), & analysis_get_n_entries (id4), & analysis_get_n_entries (id4, within_bounds = .true.) write (u, "(A,10(1x," // FMT_19 // "))") "average = ", & analysis_get_average (id1), & analysis_get_average (id2), & analysis_get_average (id3), & analysis_get_average (id3, within_bounds = .true.) write (u, "(A,10(1x," // FMT_19 // "))") "error = ", & analysis_get_error (id1), & analysis_get_error (id2), & analysis_get_error (id3), & analysis_get_error (id3, within_bounds = .true.) write (u, "(A)") write (u, "(A)") "* Clear analysis #2" write (u, "(A)") call analysis_clear (id2) do i = 4, 6 print *, "data = ", real(i,default) call analysis_record_data (id1, real(i,default)) call analysis_record_data (id2, real(i,default), & weight=real(i,default)) call analysis_record_data (id3, real(i,default)) call analysis_record_data (id4, real(i,default), real(i,default)**2) end do write (u, "(A,10(1x,I5))") "n_entries = ", & analysis_get_n_entries (id1), & analysis_get_n_entries (id2), & analysis_get_n_entries (id3), & analysis_get_n_entries (id3, within_bounds = .true.), & analysis_get_n_entries (id4), & analysis_get_n_entries (id4, within_bounds = .true.) write (u, "(A,10(1x," // FMT_19 // "))") "average = ", & analysis_get_average (id1), & analysis_get_average (id2), & analysis_get_average (id3), & analysis_get_average (id3, within_bounds = .true.) write (u, "(A,10(1x," // FMT_19 // "))") "error = ", & analysis_get_error (id1), & analysis_get_error (id2), & analysis_get_error (id3), & analysis_get_error (id3, within_bounds = .true.) write (u, "(A)") call analysis_write (u) write (u, "(A)") write (u, "(A)") "* Cleanup" call analysis_clear () call analysis_final () write (u, "(A)") write (u, "(A)") "* Test output end: analysis_1" end subroutine analysis_1 @ %def analysis_1