No OneTemporary
Actions

Size

254 KB

Subscribers

None

View Options

	Index: trunk/src/types/types.nw
	===================================================================
	--- trunk/src/types/types.nw (revision 7993)
	+++ trunk/src/types/types.nw (revision 7994)
	@@ -1,7669 +1,7673 @@
	%% -- 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]]>>=
	<<File header>>

	module particle_specifiers

	<<Use strings>>
	use io_units
	use diagnostics

	<<Standard module head>>

	<<Particle specifiers: public>>

	<<Particle specifiers: types>>

	<<Particle specifiers: interfaces>>

	contains

	<<Particle specifiers: procedures>>

	end module particle_specifiers
	@ %def particle_specifiers
	@
	\subsection{Base type}
	This is an abstract type which can hold a single particle or an expression.
	<<Particle specifiers: types>>=
	type, abstract :: prt_spec_expr_t
	contains
	<<Particle specifiers: prt spec expr: TBP>>
	end type prt_spec_expr_t

	@ %def prt_expr_t
	@ Output, as a string.
	<<Particle specifiers: prt spec expr: TBP>>=
	procedure (prt_spec_expr_to_string), deferred :: to_string
	<<Particle specifiers: interfaces>>=
	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).
	<<Particle specifiers: prt spec expr: TBP>>=
	procedure (prt_spec_expr_expand_sub), deferred :: expand_sub
	<<Particle specifiers: interfaces>>=
	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.
	<<Particle specifiers: public>>=
	public :: prt_expr_t
	<<Particle specifiers: types>>=
	type :: prt_expr_t
	class(prt_spec_expr_t), allocatable :: x
	contains
	<<Particle specifiers: prt expr: TBP>>
	end type prt_expr_t

	@ %def prt_expr_t
	@ Output as a string: delegate.
	<<Particle specifiers: prt expr: TBP>>=
	procedure :: to_string => prt_expr_to_string
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: prt expr: TBP>>=
	procedure :: init_spec => prt_expr_init_spec
	<<Particle specifiers: procedures>>=
	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
	<<Particle specifiers: prt expr: TBP>>=
	procedure :: init_list => prt_expr_init_list
	procedure :: init_sum => prt_expr_init_sum
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: prt expr: TBP>>=
	procedure :: get_n_terms => prt_expr_get_n_terms
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: prt expr: TBP>>=
	procedure :: term_to_array => prt_expr_term_to_array
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: public>>=
	public :: prt_spec_t
	<<Particle specifiers: types>>=
	type, extends (prt_spec_expr_t) :: prt_spec_t
	private
	type(string_t) :: name
	logical :: polarized = .false.
	type(string_t), dimension(:), allocatable :: decay
	contains
	<<Particle specifiers: prt spec: TBP>>
	end type prt_spec_t

	@ %def prt_spec_t
	@
	\subsubsection{I/O}
	Output. Old-style subroutines.
	<<Particle specifiers: public>>=
	public :: prt_spec_write
	<<Particle specifiers: interfaces>>=
	interface prt_spec_write
	module procedure prt_spec_write1
	module procedure prt_spec_write2
	end interface prt_spec_write
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: public>>=
	public :: prt_spec_read
	<<Particle specifiers: interfaces>>=
	interface prt_spec_read
	module procedure prt_spec_read1
	module procedure prt_spec_read2
	end interface prt_spec_read
	@ Read a single particle specifier
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: public>>=
	public :: new_prt_spec
	<<Particle specifiers: interfaces>>=
	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
	<<Particle specifiers: procedures>>=
	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
	<<Particle specifiers: prt spec: TBP>>=
	procedure :: get_name => prt_spec_get_name
	<<Particle specifiers: procedures>>=
	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
	<<Particle specifiers: prt spec: TBP>>=
	procedure :: to_string => prt_spec_to_string
	<<Particle specifiers: procedures>>=
	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
	<<Particle specifiers: prt spec: TBP>>=
	procedure :: is_polarized => prt_spec_is_polarized
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: prt spec: TBP>>=
	procedure :: is_unstable => prt_spec_is_unstable
	<<Particle specifiers: procedures>>=
	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
	<<Particle specifiers: prt spec: TBP>>=
	procedure :: get_n_decays => prt_spec_get_n_decays
	<<Particle specifiers: procedures>>=
	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
	<<Particle specifiers: prt spec: TBP>>=
	procedure :: get_decays => prt_spec_get_decays
	<<Particle specifiers: procedures>>=
	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:
	<<Particle specifiers: prt spec: TBP>>=
	procedure :: expand_sub => prt_spec_expand_sub
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: public>>=
	public :: prt_spec_list_t
	<<Particle specifiers: types>>=
	type, extends (prt_spec_expr_t) :: prt_spec_list_t
	type(prt_expr_t), dimension(:), allocatable :: expr
	contains
	<<Particle specifiers: prt spec list: TBP>>
	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.
	<<Particle specifiers: prt spec list: TBP>>=
	procedure :: to_string => prt_spec_list_to_string
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: prt spec list: TBP>>=
	procedure :: flatten => prt_spec_list_flatten
	<<Particle specifiers: procedures>>=
	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.)
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: prt spec list: TBP>>=
	procedure :: expand_sub => prt_spec_list_expand_sub
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: public>>=
	public :: prt_spec_sum_t
	<<Particle specifiers: types>>=
	type, extends (prt_spec_expr_t) :: prt_spec_sum_t
	type(prt_expr_t), dimension(:), allocatable :: expr
	contains
	<<Particle specifiers: prt spec sum: TBP>>
	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.
	<<Particle specifiers: prt spec sum: TBP>>=
	procedure :: to_string => prt_spec_sum_to_string
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: prt spec sum: TBP>>=
	procedure :: flatten => prt_spec_sum_flatten
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: prt spec sum: TBP>>=
	procedure :: expand_sub => prt_spec_sum_expand_sub
	<<Particle specifiers: procedures>>=
	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.
	<<Particle specifiers: prt expr: TBP>>=
	procedure :: expand => prt_expr_expand
	<<Particle specifiers: procedures>>=
	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]]>>=
	<<File header>>

	module particle_specifiers_ut
	use unit_tests
	use particle_specifiers_uti

	<<Standard module head>>

	<<Particle specifiers: public test>>

	contains

	<<Particle specifiers: test driver>>

	end module particle_specifiers_ut
	@ %def particle_specifiers_ut
	@
	<<[[particle_specifiers_uti.f90]]>>=
	<<File header>>

	module particle_specifiers_uti

	<<Use strings>>

	use particle_specifiers

	<<Standard module head>>

	<<Particle specifiers: test declarations>>

	contains

	<<Particle specifiers: tests>>

	end module particle_specifiers_uti
	@ %def particle_specifiers_ut
	@ API: driver for the unit tests below.
	<<Particle specifiers: public test>>=
	public :: particle_specifiers_test
	<<Particle specifiers: test driver>>=
	subroutine particle_specifiers_test (u, results)
	integer, intent(in) :: u
	type(test_results_t), intent(inout) :: results
	<<Particle specifiers: execute tests>>
	end subroutine particle_specifiers_test

	@ %def particle_specifiers_test
	@
	\subsubsection{Particle specifier array}
	Define, read and write an array of particle specifiers.
	<<Particle specifiers: execute tests>>=
	call test (particle_specifiers_1, "particle_specifiers_1", &
	"Handle particle specifiers", &
	u, results)
	<<Particle specifiers: test declarations>>=
	public :: particle_specifiers_1
	<<Particle specifiers: tests>>=
	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).
	<<Particle specifiers: execute tests>>=
	call test (particle_specifiers_2, "particle_specifiers_2", &
	"Particle specifier expressions", &
	u, results)
	<<Particle specifiers: test declarations>>=
	public :: particle_specifiers_2
	<<Particle specifiers: tests>>=
	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]]>>=
	<<File header>>

	module pdg_arrays

	use io_units
	use sorting
	use physics_defs, only: UNDEFINED

	<<Standard module head>>

	<<PDG arrays: public>>

	<<PDG arrays: types>>

	<<PDG arrays: interfaces>>

	contains

	<<PDG arrays: procedures>>

	end module pdg_arrays
	@ %def pdg_arrays
	@
	\subsection{Type definition}
	Using an allocatable array eliminates the need for initializer and/or
	finalizer.
	<<PDG arrays: public>>=
	public :: pdg_array_t
	<<PDG arrays: types>>=
	type :: pdg_array_t
	private
	integer, dimension(:), allocatable :: pdg
	contains
	<<PDG arrays: pdg array: TBP>>
	end type pdg_array_t

	@ %def pdg_array_t
	@ Output
	<<PDG arrays: public>>=
	public :: pdg_array_write
	<<PDG arrays: pdg array: TBP>>=
	procedure :: write => pdg_array_write
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: public>>=
	public :: pdg_array_write_set
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: public>>=
	public :: assignment(=)
	<<PDG arrays: interfaces>>=
	interface assignment(=)
	module procedure pdg_array_from_int_array
	module procedure pdg_array_from_int
	module procedure int_array_from_pdg_array
	end interface

	<<PDG arrays: procedures>>=
	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
	<<PDG arrays: public>>=
	public :: pdg_array_init
	<<PDG arrays: procedures>>=
	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
	<<PDG arrays: public>>=
	public :: pdg_array_delete
	<<PDG arrays: procedures>>=
	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
	<<PDG arrays: public>>=
	public :: pdg_array_merge
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: public>>=
	public :: pdg_array_get_length
	<<PDG arrays: pdg array: TBP>>=
	procedure :: get_length => pdg_array_get_length
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: public>>=
	public :: pdg_array_get
	<<PDG arrays: pdg array: TBP>>=
	procedure :: get => pdg_array_get
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg array: TBP>>=
	procedure :: set => pdg_array_set
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: pdg array: TBP>>=
	procedure :: add => pdg_array_add
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: public>>=
	public :: pdg_array_replace
	<<PDG arrays: pdg array: TBP>>=
	procedure :: replace => pdg_array_replace
	<<PDG arrays: procedures>>=
	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
	<<PDG arrays: public>>=
	public :: operator(//)
	<<PDG arrays: interfaces>>=
	interface operator(//)
	module procedure concat_pdg_arrays
	end interface

	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: public>>=
	public :: operator(.match.)
	<<PDG arrays: interfaces>>=
	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.
	<<PDG arrays: procedures>>=
	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
	<<PDG arrays: public>>=
	public :: is_quark
	<<PDG arrays: procedures>>=
	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
	<<PDG arrays: public>>=
	public :: is_gluon
	<<PDG arrays: procedures>>=
	elemental function is_gluon (pdg_nr)
	logical :: is_gluon
	integer, intent(in) :: pdg_nr
	if (pdg_nr == 21) then
	is_gluon = .true.
	else
	is_gluon = .false.
	end if
	end function is_gluon

	@ %def is_gluon
	@ Check if pdg-number corresponds to a colored particle
	<<PDG arrays: public>>=
	public :: is_colored
	<<PDG arrays: procedures>>=
	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
	<<PDG arrays: public>>=
	public :: is_lepton
	<<PDG arrays: procedures>>=
	elemental function is_lepton (pdg_nr)
	logical :: is_lepton
	integer, intent(in) :: pdg_nr
	if (abs (pdg_nr) >= 11 .and. abs (pdg_nr) <= 16) then
	is_lepton = .true.
	else
	is_lepton = .false.
	end if
	end function is_lepton

	@ %def is_lepton
	@
	<<PDG arrays: public>>=
	public :: is_fermion
	<<PDG arrays: procedures>>=
	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
	<<PDG arrays: public>>=
	public :: is_massless_vector
	<<PDG arrays: procedures>>=
	elemental function is_massless_vector (pdg_nr)
	integer, intent(in) :: pdg_nr
	logical :: is_massless_vector
	if (pdg_nr == 21 .or. pdg_nr == 22) 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
	<<PDG arrays: public>>=
	public :: is_massive_vector
	<<PDG arrays: procedures>>=
	elemental function is_massive_vector (pdg_nr)
	integer, intent(in) :: pdg_nr
	logical :: is_massive_vector
	if (abs (pdg_nr) == 23 .or. abs (pdg_nr) == 24) 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
	<<PDG arrays: public>>=
	public :: is_vector
	<<PDG arrays: procedures>>=
	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 strongly interacting
	<<PDG arrays: pdg array: TBP>>=
	procedure :: has_colored_particles => pdg_array_has_colored_particles
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: public>>=
	public :: operator(<)
	public :: operator(>)
	public :: operator(<=)
	public :: operator(>=)
	public :: operator(==)
	public :: operator(/=)
	<<PDG arrays: interfaces>>=
	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

	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: public>>=
	public :: operator(.eqv.)
	public :: operator(.neqv.)
	<<PDG arrays: interfaces>>=
	interface operator(.eqv.)
	module procedure pdg_array_equivalent
	end interface
	interface operator(.neqv.)
	module procedure pdg_array_inequivalent
	end interface

	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: public>>=
	public :: sort_abs
	<<PDG arrays: interfaces>>=
	interface sort_abs
	module procedure pdg_array_sort_abs
	end interface

	<<PDG arrays: pdg array: TBP>>=
	procedure :: sort_abs => pdg_array_sort_abs
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: pdg array: TBP>>=
	procedure :: intersect => pdg_array_intersect
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: pdg array: TBP>>=
	procedure :: search_for_particle => pdg_array_search_for_particle
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: pdg array: TBP>>=
	procedure :: invert => pdg_array_invert
	<<PDG arrays: procedures>>=
	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 (21, 22, 23, 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.
	<<PDG arrays: public>>=
	public :: pdg_list_t
	<<PDG arrays: types>>=
	type :: pdg_list_t
	type(pdg_array_t), dimension(:), allocatable :: a
	contains
	<<PDG arrays: pdg list: TBP>>
	end type pdg_list_t

	@ %def pdg_list_t
	@ Output, as a comma-separated list without advancing I/O.
	<<PDG arrays: pdg list: TBP>>=
	procedure :: write => pdg_list_write
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	generic :: init => pdg_list_init_size
	procedure, private :: pdg_list_init_size
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	generic :: init => pdg_list_init_int_array
	procedure, private :: pdg_list_init_int_array
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	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
	<<PDG arrays: procedures>>=
	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
	<<PDG arrays: pdg list: TBP>>=
	procedure :: get_size => pdg_list_get_size
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	procedure :: get => pdg_list_get
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	procedure :: is_regular => pdg_list_is_regular
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	procedure :: sort_abs => pdg_list_sort_abs
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	generic :: operator (==) => pdg_list_eq
	procedure, private :: pdg_list_eq
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	generic :: operator (<) => pdg_list_lt
	procedure, private :: pdg_list_lt
	<<PDG arrays: procedures>>=
	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]].
	<<PDG arrays: pdg list: TBP>>=
	procedure :: replace => pdg_list_replace
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: pdg list: TBP>>=
	procedure :: fusion => pdg_list_fusion
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: pdg list: TBP>>=
	procedure :: get_pdg_sizes => pdg_list_get_pdg_sizes
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	procedure :: match_replace => pdg_list_match_replace
	<<PDG arrays: procedures>>=
	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.
	<<PDG arrays: pdg list: TBP>>=
	generic :: operator (.match.) => pdg_list_match_pdg_array
	procedure, private :: pdg_list_match_pdg_array
	procedure :: find_match => pdg_list_find_match_pdg_array
	<<PDG arrays: procedures>>=
	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:
	<<PDG arrays: pdg list: TBP>>=
	procedure :: create_pdg_array => pdg_list_create_pdg_array
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: pdg list: TBP>>=
	procedure :: create_antiparticles => pdg_list_create_antiparticles
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: pdg list: TBP>>=
	procedure :: search_for_particle => pdg_list_search_for_particle
	<<PDG arrays: procedures>>=
	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
	@
	<<PDG arrays: pdg list: TBP>>=
	procedure :: contains_colored_particles => pdg_list_contains_colored_particles
	<<PDG arrays: procedures>>=
	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]]>>=
	<<File header>>

	module pdg_arrays_ut
	use unit_tests
	use pdg_arrays_uti

	<<Standard module head>>

	<<PDG arrays: public test>>

	contains

	<<PDG arrays: test driver>>

	end module pdg_arrays_ut
	@ %def pdg_arrays_ut
	@
	<<[[pdg_arrays_uti.f90]]>>=
	<<File header>>

	module pdg_arrays_uti

	use pdg_arrays

	<<Standard module head>>

	<<PDG arrays: test declarations>>

	contains

	<<PDG arrays: tests>>

	end module pdg_arrays_uti
	@ %def pdg_arrays_ut
	@ API: driver for the unit tests below.
	<<PDG arrays: public test>>=
	public :: pdg_arrays_test
	<<PDG arrays: test driver>>=
	subroutine pdg_arrays_test (u, results)
	integer, intent(in) :: u
	type (test_results_t), intent(inout) :: results
	<<PDG arrays: execute tests>>
	end subroutine pdg_arrays_test

	@ %def pdg_arrays_test
	@ Basic functionality.
	<<PDG arrays: execute tests>>=
	call test (pdg_arrays_1, "pdg_arrays_1", &
	"create and sort PDG array", &
	u, results)
	<<PDG arrays: test declarations>>=
	public :: pdg_arrays_1
	<<PDG arrays: tests>>=
	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.
	<<PDG arrays: execute tests>>=
	call test (pdg_arrays_2, "pdg_arrays_2", &
	"create and sort PDG lists", &
	u, results)
	<<PDG arrays: test declarations>>=
	public :: pdg_arrays_2
	<<PDG arrays: tests>>=
	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.
	<<PDG arrays: execute tests>>=
	call test (pdg_arrays_3, "pdg_arrays_3", &
	"check PDG lists", &
	u, results)
	<<PDG arrays: test declarations>>=
	public :: pdg_arrays_3
	<<PDG arrays: tests>>=
	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.
	<<PDG arrays: execute tests>>=
	call test (pdg_arrays_4, "pdg_arrays_4", &
	"compare PDG lists", &
	u, results)
	<<PDG arrays: test declarations>>=
	public :: pdg_arrays_4
	<<PDG arrays: tests>>=
	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.
	<<PDG arrays: execute tests>>=
	call test (pdg_arrays_5, "pdg_arrays_5", &
	"match PDG lists", &
	u, results)
	<<PDG arrays: test declarations>>=
	public :: pdg_arrays_5
	<<PDG arrays: tests>>=
	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]]>>=
	<<File header>>

	module jets

	use fastjet !NODEP!

	<<Standard module head>>

	<<Jets: public>>

	contains

	<<Jets: procedures>>

	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.
	<<Jets: public>>=
	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.
	<<Jets: procedures>>=
	subroutine fastjet_init ()
	call print_banner ()
	end subroutine fastjet_init

	@ %def fastjet_init
	@ The jet algorithm codes (actually, integers)
	<<Jets: public>>=
	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]]>>=
	<<File header>>

	module jets_ut
	use unit_tests
	use jets_uti

	<<Standard module head>>

	<<Jets: public test>>

	contains

	<<Jets: test driver>>

	end module jets_ut
	@ %def jets_ut
	@
	<<[[jets_uti.f90]]>>=
	<<File header>>

	module jets_uti

	<<Use kinds>>
	use fastjet !NODEP!

	use jets

	<<Standard module head>>

	<<Jets: test declarations>>

	contains

	<<Jets: tests>>

	end module jets_uti
	@ %def jets_ut
	@ API: driver for the unit tests below.
	<<Jets: public test>>=
	public :: jets_test
	<<Jets: test driver>>=
	subroutine jets_test (u, results)
	integer, intent(in) :: u
	type (test_results_t), intent(inout) :: results
	<<Jets: execute tests>>
	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.
	<<Jets: execute tests>>=
	call test (jets_1, "jets_1", &
	"basic FastJet functionality", &
	u, results)
	<<Jets: test declarations>>=
	public :: jets_1
	<<Jets: tests>>=
	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]]>>=
	<<File header>>

	module subevents

	use, intrinsic :: iso_c_binding !NODEP!

	<<Use kinds>>
	use io_units
	use format_defs, only: FMT_14, FMT_19
	use format_utils, only: pac_fmt
	use sorting
	use c_particles
	use lorentz
	use pdg_arrays
	use jets

	<<Standard module head>>

	<<Subevents: public>>

	<<Subevents: parameters>>

	<<Subevents: types>>

	<<Subevents: interfaces>>

	contains

	<<Subevents: procedures>>

	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.
	<<Subevents: parameters>>=
	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.
	<<Subevents: public>>=
	public :: prt_t
	<<Subevents: types>>=
	type :: prt_t
	private
	integer :: type = PRT_UNDEFINED
	integer :: pdg
	logical :: polarized = .false.
	integer :: h
	type(vector4_t) :: p
	real(default) :: p2
	integer, dimension(:), allocatable :: src
	end type prt_t

	@ %def prt_t
	@ Initializers. Polarization is set separately. Finalizers are not
	needed.
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: prt_init_combine
	<<Subevents: procedures>>=
	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.
	<<Subevents: procedures>>=
	subroutine prt_init_pseudojet (prt, jet, src, pdg)
	type(prt_t), intent(out) :: prt
	type(pseudojet_t), intent(in) :: jet
	integer, dimension(:), intent(in) :: src
	integer, intent(in) :: pdg
	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)
	end subroutine prt_init_pseudojet

	@ %def prt_init_pseudojet
	@
	\subsubsection{Accessing contents}
	<<Subevents: public>>=
	public :: prt_get_pdg
	<<Subevents: procedures>>=
	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
	<<Subevents: public>>=
	public :: prt_get_momentum
	<<Subevents: procedures>>=
	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
	<<Subevents: public>>=
	public :: prt_get_msq
	<<Subevents: procedures>>=
	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
	<<Subevents: public>>=
	public :: prt_is_polarized
	<<Subevents: procedures>>=
	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
	<<Subevents: public>>=
	public :: prt_get_helicity
	<<Subevents: procedures>>=
	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
	@
	\subsubsection{Setting data}
	Set the PDG, momentum and momentum squared, and ancestors. If
	allocate-on-assignment is available, this can be simplified.
	<<Subevents: procedures>>=
	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.
	<<Subevents: procedures>>=
	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.
	<<Subevents: procedures>>=
	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.
	<<Subevents: procedures>>=
	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).
	<<Subevents: procedures>>=
	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
	@
	\subsubsection{Conversion}
	Transform a [[prt_t]] object into a [[c_prt_t]] object.
	<<Subevents: public>>=
	public :: c_prt
	<<Subevents: interfaces>>=
	interface c_prt
	module procedure c_prt_from_prt
	end interface

	@ %def c_prt
	<<Subevents: procedures>>=
	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}
	<<Subevents: public>>=
	public :: prt_write
	<<Subevents: procedures>>=
	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)
	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
	write (u, "(A)") ")"
	end subroutine prt_write

	@ %def prt_write
	@
	\subsubsection{Tools}
	Two particles match if their [[src]] arrays are the same.
	<<Subevents: interfaces>>=
	interface operator(.match.)
	module procedure prt_match
	end interface

	@ %def .match.
	<<Subevents: procedures>>=
	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.
	<<Subevents: procedures>>=
	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
	call combine_index_lists (src, prt_in1%src, prt_in2%src)
	ok = allocated (src)
	if (ok) call prt_init_composite (prt, prt_in1%p + prt_in2%p, src)
	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).
	<<Subevents: public>>=
	public :: are_disjoint
	<<Subevents: procedures>>=
	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.
	<<Subevents: procedures>>=
	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).
	<<Subevents: procedures>>=
	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}
	<<Subevents: public>>=
	public :: subevt_t
	<<Subevents: types>>=
	type :: subevt_t
	private
	integer :: n_active = 0
	type(prt_t), dimension(:), allocatable :: prt
	contains
	<<Subevents: subevt: TBP>>
	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.
	<<Subevents: public>>=
	public :: subevt_init
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_reset
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_write
	<<Subevents: subevt: TBP>>=
	procedure :: write => subevt_write
	<<Subevents: procedures>>=
	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).
	<<Subevents: interfaces>>=
	interface assignment(=)
	module procedure subevt_assign
	end interface

	@ %def =
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_set_beam
	public :: subevt_set_incoming
	public :: subevt_set_outgoing
	public :: subevt_set_composite
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_set_pdg_beam
	public :: subevt_set_pdg_incoming
	public :: subevt_set_pdg_outgoing
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_set_p_beam
	public :: subevt_set_p_incoming
	public :: subevt_set_p_outgoing
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_set_p2_beam
	public :: subevt_set_p2_incoming
	public :: subevt_set_p2_outgoing
	<<Subevents: procedures>>=
	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
	<<Subevents: public>>=
	public :: subevt_polarize
	<<Subevents: procedures>>=
	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
	@
	\subsubsection{Accessing contents}
	Return true if the subevent has entries.
	<<Subevents: public>>=
	public :: subevt_is_nonempty
	<<Subevents: procedures>>=
	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
	<<Subevents: public>>=
	public :: subevt_get_length
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_get_prt
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_get_sqrts_hat
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_get_n_in
	public :: subevt_get_n_out
	<<Subevents: procedures>>=
	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
	@
	<<Subevents: interfaces>>=
	interface c_prt
	module procedure c_prt_from_subevt
	module procedure c_prt_array_from_subevt
	end interface

	@ %def c_prt
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_join
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_combine
	<<Subevents: procedures>>=
	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.)
	<<Subevents: public>>=
	public :: subevt_collect
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_cluster
	<<Subevents: procedures>>=
	subroutine subevt_cluster (subevt, pl1, mask1, jet_def, keep_jets)
	type(subevt_t), intent(inout) :: subevt
	type(subevt_t), intent(in) :: pl1
	logical, dimension(:), intent(in) :: mask1
	type(jet_definition_t), intent(in) :: jet_def
	logical, intent(in) :: keep_jets
	integer, dimension(:), allocatable :: src, src_tmp
	integer, dimension(:), allocatable :: map, jet_idx
	type(pseudojet_t), dimension(:), allocatable :: jet_in, jet_out
	type(pseudojet_vector_t) :: jv_in, jv_out
	type(cluster_sequence_t) :: cs
	integer :: i, prt_idx, k, n_src, n_active, this_jet, combined_pdg, pdg, n_quarks
	n_active = 0
	allocate (map (pl1%n_active), source = 0)
	allocate (src (0))
	do i = 1, pl1%n_active
	if (mask1(i)) then
	call combine_index_lists (src_tmp, src, pl1%prt(i)%src)
	if (allocated (src_tmp)) then
	call move_alloc (from=src_tmp, to=src)
	n_active = n_active + 1
	map(n_active) = i
	end if
	end if
	end do
	allocate (jet_in (count (map /= 0)))
	do i = 1, size (jet_in)
	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)
	jv_out = cs%inclusive_jets ()
	allocate (jet_idx (size (jet_in)))
	call cs%assign_jet_indices (jv_out, jet_idx)
	allocate (jet_out (jv_out%size ()))
	jet_out = jv_out
	call subevt_reset (subevt, size (jet_out))
	do this_jet = 1, size (jet_out)
	src = 0
	n_src = 0
	combined_pdg = 0
	n_quarks = 0
	do prt_idx = 1, size (jet_idx)
	if (jet_idx(prt_idx) == this_jet) then
	associate (prt => pl1%prt(map(prt_idx)))
	do k = 1, size (prt%src)
	src(n_src + k) = prt%src(k)
	end do
	n_src = n_src + size (prt%src)
	if (is_quark (prt%pdg)) then
	n_quarks = n_quarks + 1
	if (combined_pdg == 0) then
	combined_pdg = prt%pdg
	end if
	end if
	end associate
	end if
	end do
	if (keep_jets .and. n_quarks == 1) then
	pdg = combined_pdg
	else
	pdg = 0
	end if
	call prt_init_pseudojet (subevt%prt(this_jet), jet_out(this_jet), src(:n_src), pdg)
	end do
	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
	end subroutine subevt_cluster

	@ %def subevt_cluster
	@ Return a list of all particles for which the mask is true.
	<<Subevents: public>>=
	public :: subevt_select
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_extract
	<<Subevents: procedures>>=
	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.
	<<Subevents: public>>=
	public :: subevt_sort
	<<Subevents: interfaces>>=
	interface subevt_sort
	module procedure subevt_sort_pdg
	module procedure subevt_sort_int
	module procedure subevt_sort_real
	end interface

	<<Subevents: procedures>>=
	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)))
	- subevt%prt = pl%prt(order (rval))
	+ 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).

	The [[pack]] command was buggy in some gfortran versions, therefore it
	is unrolled. The unrolled version may be more efficient, actually.
	<<Subevents: public>>=
	public :: subevt_select_pdg_code
	<<Subevents: procedures>>=
	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)
	!!! !!! !!! Workaround for gfortran compiler bug
	! subevt%prt(:n_match) = pack (subevt_in%prt(:n_active), mask)
	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.
	<<Subevents: public>>=
	public :: pacify
	<<Subevents: interfaces>>=
	interface pacify
	module procedure pacify_prt
	module procedure pacify_subevt
	end interface pacify

	@ %def pacify
	<<Subevents: procedures>>=
	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]]>>=
	<<File header>>

	module analysis

	<<Use kinds>>
	<<Use strings>>
	use io_units
	use format_utils, only: quote_underscore, tex_format
	use system_defs, only: TAB
	use diagnostics
	use os_interface
	use ifiles

	<<Standard module head>>

	<<Analysis: public>>

	<<Analysis: parameters>>

	<<Analysis: types>>

	<<Analysis: interfaces>>

	<<Analysis: variables>>

	contains

	<<Analysis: procedures>>

	end module analysis
	@ %def analysis
	@
	\subsection{Output formats}
	These formats share a common field width (alignment).
	<<Analysis: parameters>>=
	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.
	<<Analysis: public>>=
	public :: graph_options_t
	<<Analysis: types>>=
	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.
	<<Analysis: public>>=
	public :: graph_options_init
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: graph_options_set
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: graph_options_write
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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).
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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).
	<<Analysis: public>>=
	public :: drawing_options_t
	<<Analysis: types>>=
	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.
	<<Analysis: public>>=
	public :: drawing_options_write
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: drawing_options_init_histogram
	public :: drawing_options_init_plot
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: drawing_options_set
	<<Analysis: procedures>>=
	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:
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: types>>=
	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
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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 a value. Always successful for observables.
	<<Analysis: interfaces>>=
	interface observable_record_value
	module procedure observable_record_value_unweighted
	module procedure observable_record_value_weighted
	end interface

	<<Analysis: procedures>>=
	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 + abs (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}
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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}
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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}
	<<Analysis: types>>=
	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
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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 + abs (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
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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}
	<<Analysis: types>>=
	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.
	<<Analysis: interfaces>>=
	interface histogram_init
	module procedure histogram_init_n_bins
	module procedure histogram_init_bin_width
	end interface

	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.)
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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}
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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}
	<<Analysis: types>>=
	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
	<<Analysis: interfaces>>=
	interface point_init
	module procedure point_init_contents
	module procedure point_init_point
	end interface
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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}
	<<Analysis: types>>=
	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.
	<<Analysis: interfaces>>=
	interface plot_init
	module procedure plot_init_empty
	module procedure plot_init_plot
	end interface
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: types>>=
	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
	<<Analysis: procedures>>=
	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:
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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.
	<<Analysis: 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.
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.)
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.)
	<<Analysis: parameters>>=
	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
	<<Analysis: public>>=
	public :: AN_UNDEFINED, AN_HISTOGRAM, AN_OBSERVABLE, AN_PLOT, AN_GRAPH

	@ %def AN_UNDEFINED
	@ %def AN_OBSERVABLE AN_HISTOGRAM AN_PLOT AN_GRAPH
	<<Analysis: types>>=
	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.
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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:
	<<Analysis: procedures>>=
	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:
	<<Analysis: procedures>>=
	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}
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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).
	<<Analysis: public>>=
	public :: analysis_iterator_t
	<<Analysis: types>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_iterator_is_valid
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_iterator_advance
	<<Analysis: procedures>>=
	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:
	<<Analysis: public>>=
	public :: analysis_iterator_get_type
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_iterator_get_data
	<<Analysis: procedures>>=
	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.
	<<Analysis: variables>>=
	type(analysis_store_t), save :: analysis_store

	@ %def analysis_store
	<<Analysis: types>>=
	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
	<<Analysis: public>>=
	public :: analysis_final
	<<Analysis: procedures>>=
	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
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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
	<<Analysis: public>>=
	public :: analysis_store_get_object_type
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_store_get_ids
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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.
	<<Analysis: procedures>>=
	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:
	<<Analysis: public>>=
	public :: analysis_init_observable
	<<Analysis: procedures>>=
	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
	<<Analysis: public>>=
	public :: analysis_init_histogram
	<<Analysis: interfaces>>=
	interface analysis_init_histogram
	module procedure analysis_init_histogram_n_bins
	module procedure analysis_init_histogram_bin_width
	end interface

	<<Analysis: procedures>>=
	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
	<<Analysis: public>>=
	public :: analysis_init_plot
	<<Analysis: procedures>>=
	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
	<<Analysis: public>>=
	public :: analysis_init_graph
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_clear
	<<Analysis: interfaces>>=
	interface analysis_clear
	module procedure analysis_store_clear_obj
	module procedure analysis_store_clear_all
	end interface

	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_record_data
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_fill_graph
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_exists
	<<Analysis: procedures>>=
	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:
	<<Analysis: public>>=
	public :: analysis_get_n_elements
	public :: analysis_get_n_entries
	public :: analysis_get_average
	public :: analysis_get_error
	<<Analysis: procedures>>=
	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
	<<Analysis: public>>=
	public :: analysis_has_plots
	<<Analysis: interfaces>>=
	interface analysis_has_plots
	module procedure analysis_has_plots_any
	module procedure analysis_has_plots_obj
	end interface

	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_init_iterator
	<<Analysis: procedures>>=
	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}
	<<Analysis: public>>=
	public :: analysis_write
	<<Analysis: interfaces>>=
	interface analysis_write
	module procedure analysis_write_object
	module procedure analysis_write_all
	end interface

	@ %def interface
	<<Analysis: procedures>>=
	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
	<<Analysis: public>>=
	public :: analysis_write_driver
	<<Analysis: procedures>>=
	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
	<<Analysis: public>>=
	public :: analysis_compile_tex
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_get_header
	<<Analysis: procedures>>=
	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.
	<<Analysis: public>>=
	public :: analysis_write_makefile
	<<Analysis: procedures>>=
	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]]>>=
	<<File header>>

	module analysis_ut
	use unit_tests
	use analysis_uti

	<<Standard module head>>

	<<Analysis: public test>>

	contains

	<<Analysis: test driver>>

	end module analysis_ut
	@ %def analysis_ut
	@
	<<[[analysis_uti.f90]]>>=
	<<File header>>

	module analysis_uti

	<<Use kinds>>
	<<Use strings>>
	use format_defs, only: FMT_19

	use analysis

	<<Standard module head>>

	<<Analysis: test declarations>>

	contains

	<<Analysis: tests>>

	end module analysis_uti
	@ %def analysis_ut
	@ API: driver for the unit tests below.
	<<Analysis: public test>>=
	public :: analysis_test
	<<Analysis: test driver>>=
	subroutine analysis_test (u, results)
	integer, intent(in) :: u
	type(test_results_t), intent(inout) :: results
	<<Analysis: execute tests>>
	end subroutine analysis_test

	@ %def analysis_test
	<<Analysis: execute tests>>=
	call test (analysis_1, "analysis_1", &
	"check elementary analysis building blocks", &
	u, results)
	<<Analysis: test declarations>>=
	public :: analysis_1
	<<Analysis: tests>>=
	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

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:46 PM (1 d, 21 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805016
Default Alt Text: (254 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions