Chapter 4  Running WHIZARD

• Setting up the process list
• Compilation, installation, and bundles
• Process selection
• Input data
• Cuts
• Integration
• Event generation
• Fragmentation
• Controlling WHIZARD runs

4.1  Setting up the process list

WHIZARD is able to cope with arbitrary scattering processes allowed by the selected physical model. The program is not a library of processes, but the user can compile a list of processes he is interested in, write a file whizard.prc and put it in the subdirectory conf. If no file is provided by the user, the file whizard.prc.default (Fig. 4.1) will be used. Note that as per default CompHEP and MadGraph are disabled all processes by these generators will be skipped. Only O’Mega and trivial test matrix elements will be produced. The user may take this file as a template, insert or delete processes as desired, and save the new version as whizard.prc.¹

---

                 # WHIZARD configuration file
                 
                 # The selected model
                 model   SM
                 
                 # Processes
                 #  Methods: chep=CompHEP, mad=MadGraph, omega=O'Mega, test=trivial)
                 #  Options: s      selected diagrams (CompHEP/MadGraph)
                 #           r      restricted intermediate state (O'Mega)
                 #           c      apply exact color algebra (O'Mega)
                 #           n:XXX  coupling order (MadGraph)
                 #           w:XXX  width scheme (O'Mega)
                 #           p      transfer polarization (test)
                 #           u      unit matrix element (test)
                 #
                 # Tag           In      Out             Method  Option
                 #=====================================================
                 ee_c            e1,E1   e1,E1           chep
                 ee_m            e1,E1   e1,E1           mad
                 ee_o            e1,E1   e1,E1           omega
                 ww_c            e1,E1   W+,W-           chep
                 ww_m            e1,E1   W+,W-           mad
                 ww_o            e1,E1   W+,W-           omega
                 zh_c            e1,E1   Z,H             chep
                 zh_m            e1,E1   Z,H             mad
                 zh_o            e1,E1   Z,H             omega
                 nnh_c           e1,E1   n1,N1,H         chep
                 nnh_m           e1,E1   n1,N1,H         mad
                 nnh_o           e1,E1   n1,N1,H         omega
                 nnbb_m          e1,E1   n1,N1,b,B       mad
                 nnbb_o          e1,E1   n1,N1,b,B       omega
                 

Figure 4.1: Default process configuration file. Only O’Mega processes will be generated per default, the others will be skipped.

---

4.1.1  Model selection

WHIZARD is not restricted to the Standard Model. The physics model is selected in the header of whizard.prc, respectively.

For each model, all particles, parameters, and the vertices which go into the phase space setup are defined in a single file with extension .mdl. The vertex definitions include the complete specification of Feynman rules for MadGraph and CompHEP, while the O’Mega vertices are currently hard-coded in the corresponding executable. The model files can be found in the conf/models subdirectory. Along with the master files SM.mdl, MSSM.mdl, etc., for each model XXX there is also a file parameters.XXX.omega.f90 which contains the translation to the parameter definitions used internally by O’Mega.

The distribution contains the following models (see Tab. 4.2):

• QED: QED with three lepton generations.
• QCD: QCD with three quark generations.
• SM: The Standard Model with e, sinθ_W and M_Z as independent parameters in the electroweak sector. Some remarks are in order²:
  • The u and d masses are zero, all other masses are non-zero by default and user-definable.
  • All fermions, except for the top quark, have zero width.
  • By default, the CKM matrix is the unit matrix, but for each model where it matters, a variant with nontrivial CKM matrix is available. Keep in mind that keeping the full CKM matrix may result in less efficient calculations.
  • The Higgs boson is coupled to all massive fermions, with one exception: O’Mega has no Higgs couplings to the strange quark.
• SM_ac: The Standard Model with anomalous couplings. This is usually considered in the context of dynamical electroweak symmetry breaking, hence, by default, the Higgs mass is set to a large value. O’Mega supports anomalous trilinear and quartic gauge couplings³.
• SM_km: This model has been included in version 1.90. It is intended for LHC or CLIC physics of Higgsless models or models without unitarizing UV completion. A K-matrix unitarization guarantees a physically sane description for LHC and CLIC energies. More details can be found in [13].
• MSSM: The Minimal Supersymmetric Standard Model (MSSM)⁴. The input data may be provided in the SUSY Les Houches Accord format (see Sec. 4.4.7). Typically, such a file is generated by automatic SUSY spectrum and decay packages.

  For convenience, there are two MSSM input files provided in the conf directory: sps1a.in implements the mass spectrum of the standard SUSY point SPS1a [14] with zero particle widths; this file has been used for generating the comparison results of Ref. [12]. The file sps1ap.in corresponds to the reference point SPS1a’ as it has been defined for the SPA convention [15]; this file includes width values for all particles that are computed, in some cases, with higher-order corrections and running couplings incorporated. Be aware that, for a consistent tree-level calculation, it may be necessary to compute widths in a different convention. This caveat is discussed further in Ref. [12].

• MSSM_Grav: This variant of the MSSM has been included in version 1.93. It contains the MSSM with an additional gravitino and is intended for the studies of GMSB models.
• NMSSM: This is the Next-to-Minimal Supersymmetric SM, which has been contributed to WHIZARD by Felix Braam [16].
• PSSSM: A non-minimal MSSM version motivated by an E₆ → Pati-Salam GUT group. This is a very extended model which has been implemented by Daniel Wiesler.

---

Model type	with CKM matrix	trivial CKM
QED with e,µ,τ,γ	–	QED
QCD with d,u,s,c,b,t,g	–	QCD
Standard Model	SM_CKM	SM
SM with anomalous couplings	–	SM_ac
SM with K matrix	–	SM_km
MSSM	MSSM_CKM	MSSM
MSSM with gravitino	–	MSSM_Grav
NMSSM	NMSSM_CKM	NMSSM
E6 SSM	–	PSSSM
Littlest Higgs model	–	Littlest
Littlest Higgs model (ungauged U(1))	–	Littlest_Eta
Simplest Little Higgs model	–	Simplest
Simplest Little Higgs (universal coupl.)	–	Simplest_univ
UED	–	UED
SM with Z′	–	Zprime
SM with graviton resonance	–	Xdim
SUSY toy model with gravitino	–	GravTest
SM as template	–	Template
[SM, all unitary gauge, CompHEP model	–	SM_ug ]

Figure 4.2: Summary of models currently supported by WHIZARD.

---

There are more BSM models, whose inclusion into WHIZARD started from version 1.51 on:

• Littlest: The Littlest Higgs model, which contains in addition to the SM a B′ (or γ′), a W′/Z′ triplet, a complex isospin-2 scalar multiplet ψ, and a t′.
• Littlest_Eta: The variant of the Littlest Higgs model which contains a pseudoscalar η instead of the B′ vector boson.
• Simplest: A different Little Higgs model which has recurrences of all SM fermions and new isodoublet gauge bosons.
• Simplest_univ: A variant of the previous model with universal, but not anomaly-free, fermion assignments.
• Xdim: The SM augmented by a spin-2 graviton as a toy model
• GravTest: The SM augmented by a photino and a gravitino as a toy model.
• UED: Universal Extra Dimensions.
• Zprime: The SM with a completely generic Z′ resonance.
• Template: This is just a copy of the SM. It might serve as a template if the user wants to implement his/her own model.

The user might wonder whether this list can be extended. Actually, it is trivial to set up, e.g., a SM version with a different input parameter scheme or different particle naming conventions, etc.: Copy SM.mdl to, e.g., SM_new.mdl, make the changes in that file, and provide copies of the corresponding files for use of the matrix element generators⁵. There is a model Template as described above, which is just the SM and might serve as a template to implement the user’s own model. However, for this task one needs to modify the module Template in omega-src/bundle/src/models5.ml (which requires programming in O’Caml) and then the corresponding files Template.mdl and parameters.Template.omega.f90 in conf/models. This is a rather cumbersome task. There is fortunately an easier way: using the interface of WHIZARD to the program FeynRules [17]. Here, one can write down a Mathematica file containing a Lagrangian in a very intuitive language, and then FeynRules derives the mass eigenstates, mixing and Feynman rules and imports the model into WHIZARD. Although this has been developped already for WHIZARD version 1 and is fully supported and functional for that version, we refer to the manual of version 2 for more information, or alternatively the WHIZARD Wiki page http://projects.hepforge.org/whizard/trac/wiki.

4.1.2  Process list

WHIZARD is capable of computing 2→ n scattering processes and 1→ n decay processes. In the former case, the overall result (the integral) displayed by the program will be the total cross section in fb. In the latter case, it will be the partical decay width in GeV.

For scattering processes, n=1 (single-particle production) is allowed⁶⁷. Obviously, this requires at least one beam to be structured, such that a energy spectrum is available for the initial state.

Each process is defined by a single line in whizard.prc, for instance:

h_production    e1,E1   n1,N1,H         omega
z_decay         Z       e1,E1           omega
z_production    e1,E1   Z               omega

(see Fig. 4.1). The meaning of the entries within a line is as follows:

The first entry is a unique alphanumeric tag which will be used to identify the process. It can be arbitrarily chosen, but should be a valid Fortran 95 token, lower case only. In particular, one may use lower-case characters and digits as well as the underscore, where the first character is not a digit.

The second entry defines the (partonic) initial state. It consists of two particle names, separated by a comma, with no space in between. The particle names are those given in the model file mentioned above, where typically for each particle several possible names are listed. As an example, in Fig. 4.3, the Standard Model particles are displayed as they appear in SM.mdl. This file also contains alternative names, quantum number assignments, and the physical parameters of the model.

---

Particle	Name	Antiparticle	Mass	Width
d-quark	d	D
u-quark	u	U
s-quark	s	S	ms
c-quark	c	C	mc
b-quark	b	B	mb
t-quark	t	T	mtop	wtop
electron	e1	E1	me
muon	e2	E2	mmu
tau-lepton	e3	E3	mtau
e-neutrino	n1	N1
µ-neutrino	n2	N2
τ-neutrino	n3	N3
gluon	G
photon	A
Z-boson	Z		mZ	wZ
W-boson	W+	W-	mW	wW
Higgs	H		mH	wH

Figure 4.3: Particle names (CompHEP naming scheme) and parameters for the Standard Model.

---

The third entry defines the final state. It consists of one, two or more particle names, separated by commas, with no space in between.

Flavor sums⁸ are defined by particle names separated by colons. For instance, the line

  eeqq  e1,E1   u:d:s,U:D:S     omega

will create a sum of all processes e⁻e⁺→ qq, with q being any of the light quarks. Similarly,

  qqmm  u:d:s,U:D:S   e2,E2     omega

creates a sum of all processes qq→ µ⁻µ⁺. This is especially suitable for hadronic collisions. One should note that the mirror process qq→ µ⁻µ⁺ is not automatically included. To get the correct results for a hadron collider, one can either double cross section values and symmetrize final-state distributions⁹ (for a symmetric initial state, e.g. pp) or define the process as¹⁰

  qqmm u:d:s:U:D:S,u:d:s:U:D:S e2,E2 omega

Note that in the diagnostics messages at runtime, only a single flavor assignment will be shown, which is used as a template for generating phase space parameterizations. For the matrix element evaluation and event generation, however, all valid flavor assignments are used.

To make this more readable, there is the possibility to specify aliases: The above definition is equivalent to

  alias  q  u:d:s:U:D:S
  qqmm   q,q   e2,E2    omega

Alias definitions can be arbitrarily mixed with process definitions. They are applied consecutively to all in- and out-state specifications below their definition.

The particles to be summed over must have identical spin, mass and color representation; otherwise the summation will not work. This feature is therefore useful mainly for summing over fermions with identical quantum numbers (quarks or leptons).

The order of particles in the final state can matter if PYTHIA/JETSET is used for fragmenting events. If no information on color connections is available in simulated events (for O’Mega matrix elements if the ’c’ option is not set ¹¹), the (dominant) color configuration has to be inferred from the particle ordering, if at all possible ¹². To transfer the color connections to JETSET, particles that make up a color-singlet string should be put together, beginning with the quark, possibly intermediate gluons, and ending with the antiquark. This suffices for well-defined flavor states. If flavor sums are used and more than one quark pair is present, all possible orders of particles will show up in the final state.

The fourth entry selects the method for the generation of the matrix element, while the fifth entry is either absent or consists of a string of code letters separated by commas (blanks in between are allowed). Some code letters are followed by a colon ’:’ that indicates a value. As a default there are two options for matrix elements, either the method omega using the intrinsic generator O’Mega or the method test which produces dummy matrix elements for certain consistency tests described below.

For cross checks and debugging there are also (when enabled during the configuration of WHIZARD) the methods chep and mad for using matrix elements generated by CompHEP or MadGraph, respectively. These features have become obsolete for the user at latest with version 1.50 or 1.90, when O’Mega become the main or default generator WHIZARD. For completeness reasons the details and options for these two generators will be described at the end of this subsection.

O’Mega can generate matrix elements with an arbitrary number of final-state particles, limited only by the host resources (WHIZARD has extensively been tested only for processes up to 2→ 6, some very special cases of 2→ 7 or 2→ 8 have been studied). For multiparticle processes, code generated by O’Mega is expected to be the fastest, since more redundancies are eliminated from the matrix element evaluation.

There are several options (entries in the fifth column) for the O’Mega matrix element generator:

• r: The letter r indicates restrictions on intermediate states for the amplitude. Effectively, this also selects classes of Feynman diagrams. Example:

      nnh        e1,E1   n1,N1,H   omega
      nnh_zh     e1,E1   n1,N1,H   omega  r: 3+4~Z
    

  For more details, see Sec. 5.2.
• w: The letter w modifies the width treatment for unstable particles. The default behavior is the expected one: the width is inserted only for s-channel propagators. This may lead to gauge-invariance issues. The option w:fudged (or abbreviated w:f) applies the ’fudge-factor’ prescription for unstable particles where all terms are collected over the respective propagator denominators (in contrast to the usual fixed-width prescription). This prescription is useful if the process includes the exchange of photons and W bosons in the t-channel, for instance:

      enw        e1,E1  n1,E1,W-    omega    w:fudged
    

  A fixed-width option can be enforced by specifying w:constant (or w:c). Finally, all widths can be forced to zero both in the s- and t-channel by the option w:zero (or w:z).
• O’Mega comes in two versions: By default, only the leading-1/N_c color factor is generated, and no color-flow information is available. If the option flag c is given, WHIZARD uses the tensor-product feature for process definitions of O’Mega originally intended for flavor summation to treat color exactly and to generate color-flow information. Example:

      uudd_leading   u,U  d,D      omega
      uudd_full      u,U  d,D      omega   c
    

  The modified version also eliminates some redundancies in summing over flavor.

  In the current implementation, the default (leading-color) version can be more efficient (unless many flavor states need to be summed over), so it is still useful if color is not an issue. This is the reason why we decided to keep the leading-color version as the default. Note that the O’Mega generator within WHIZARD 2 takes care of the color-flow information completely by itself and delivers matrix elements that contain the full color correlations as a default.

The test method generates a constant matrix element, independent of the particles supplied for the process. There are two choices for the normalization:

1. If the option ’u’ is given, the matrix element is unity. With zero masses and no cuts, the integral should be exactly

   σn(s) =

   (2π)4 2 s

   Φn(s) =

   1 16π s

   Φn(s) Φ2(s)

   ,     (1)

   where the volume of the massless n-particle phase space is given by

   Φn(s) =

   1 4(2π)5

   ⎛ ⎜ ⎜ ⎝

   s 16π2

   ⎞ ⎟ ⎟ ⎠

   n−2

   1 (n−1)!(n−2)!

   .     (2)

   For n≠2 the phase space volume is dimensionful, so the units of the integral are fb×GeV²⁽ⁿ⁻²⁾.
2. Without option (default), the result of the integration is not a cross section, but the phase space volume, normalized by the massless phase space volume (2). Thus, if all final-state particles are massless and no cuts are specified, the integral is unity regardless of the number of particles.

Cuts and nonvanishing masses will affect the result accordingly. If events are generated, the final-state particles are distributed uniformly in phase space.

Note that the phase-space integration for the test matrix element is organized in the same way as it would be for the real 2→ n process. Since such a phase space parameterization is not optimized for the constant matrix element that is supplied instead, good convergence is not guaranteed. (Setting stratified = F may be helpful here.)

The possibility to call a dummy matrix element in this way allows to histogram spectra or structure functions: Choose a trivial process such as uu→ dd, select the test method, switch on structure functions for one (or both) beams, and generate events. The distribution of the final-state mass squared reflects the x dependence of the selected structure function.

For the test method there are two more options:

• p:: The letter p indicates that beam polarization should be transferred to the final state unchanged. (Otherwise, it would be ignored.) This is only possible if initial and final state coincide, e.g, e⁺e⁻→ e⁺e⁻. The option is useful for testing polarized spectra and structure functions.

      ee_test_pol    e1,E1  e1,E1   test   p
    

• u:: The letter u indicates that the matrix element is not just constant, but unity. Such code may serve as the initial setup for inserting some nontrivial matrix element by hand.

For historical and completeness reasons we mention here the technical details and available options when using CompHEP or MadGraph matrix elements.

The CompHEP version included in WHIZARD can only generate matrix elements for up to four particles in the final state. No polarization or flavor summation is possible (in contrast to the official CompHEP version), and color-flow information is not explicitly generated in the version used by WHIZARD. However, arbitrary models can be defined and the predefined models can be modified to include anomalous couplings etc. CompHEP matrix elements, evaluated by the trace method, are most efficient for simple 2→ 2 or 2→ 3 processes, while for 2→ 4 processes the other two packages typically generate faster code.

The MadGraph version inside WHIZARD has the Standard Model implemented with the exception of the trilinear or quartic Higgs couplings. The user may choose the maximal order of the strong (QCD), electromagnetic (QED), or weak (QFD) coupling for the matrix element evaluation. The implementation included in WHIZARD can generate matrix elements with up to six final-state particles. No flavor summation is possible.

There are several options in the fifth entry of the input file for the two generators: The letter s indicates that only a selection of particular Feynman diagrams is taken. Example:

    nnh        e1,E1   n1,N1,H   chep
    nnh_zh     e1,E1   n1,N1,H   chep   s
  

The setup of the necessary diagram selection file is described in Sec. 5.1.

For MadGraphonly it is possible to define the maximal order in each of the coupling constants: ‘QCD’ counts powers of g_s (gluon couplings), ‘QED’ counts powers of e (photon couplings), and ‘QFD’ counts powers of the remaining electroweak couplings. As an example, the following specifications generate the QCD and the electroweak contribution to uū→ dd separately:

    uudd_qcd   u,U   d,D    mad    n:QED<=0, n:QFD<=0
    uudd_ew    u,U   d,D    mad    n:QCD<=0

There is only an upper bound on the number of couplings of a given species supported, so <= is the only available relation.

MadGraph always uses a fixed-width prescription, both in the s- and the t-channel, so there is no option. The “fudged” width option mentioned above can be activated for CompHEP by setting the flag gwidth. However, this is a runtime flag that enters the input file whizard.in in the block parameter_input.

4.2  Compilation, installation, and bundles

After the processes have been selected, the matrix elements and executable code are compiled and linked by the simple command

make prg

This will generate an executable named whizard in the bin subdirectory.

An alternative is

make bundle-src
make bundle

which will create a restricted bundle containing the Fortran 95 source code of the selected matrix elements and the libraries which have not been explicitly disabled.¹³ This bundle will be put as a gzipped tar file in the main directory under the name whizard-bundle-yymmdd-hhmm.tgz, where yymmdd-hhmm specifies the current date and time. The bundle can be unpacked on a different platform. There, it can be configured, compiled and run without the need for a working O’Mega (i.e. O’Caml) installation. Important: In the bundle, you have to use the commands make prg_bundle install_bundle instead of make prg install, otherwise it won’t work. (The necessary O’Mega library will be contained in the bundle, of course.) The configuration and input files (whizard.prc, whizard.in, etc.) will be copied from the conf directory of the distribution (not from the results directory). The site-specific configuration config.site will not be copied (precisely because it is site-specific). The restricted bundle has all capabilities of the WHIZARD system except for the possibility to re-generate the matrix elements.

If only the source code for the process list should be generated, say

make proc-src

The files can be found in the subdirectory processes-src. If desired, the source code can be modified before a further make call compiles and links them together.

A run of the Monte Carlo generator requires, apart from the executable, the following input and configuration files. In their initial form they are found in the conf subdirectory:

• whizard.in containing all input data and parameters.
• whizard.mdl containing vertex definitions used in phase space setup. If an appropriate phase space configuration file whizard.phs exists, the model file is not necessary.
• whizard.cut1 containing a-priori cut definitions. This file may be empty.
• whizard.cut5 containing cut and histogram definitions for the data analysis step (the 5th pass of the program). This file may be empty.

When the command

make install

is issued, all necessary configuration files, together with the whizard executable, are copied into the results subdirectory. If no user-defined files are found, default ones will be used. The contents of the results directory may be copied anywhere without affecting the functionality of the whizard executable within, leaving room for the simultaneous use of multiple WHIZARD copies with different process content. (’make install’ implies ’make prg’, so the latter is in fact redundant.)

Alternatively, one can omit the install step and copy the WHIZARD executable to any place in the system, or just leave it in the bin subdirectory. For a successful run of WHIZARD, one needs, of course, copies of the input files in the working directory.

Dealing with compilation problems

Although the authors have tried to make the program compile and run in various environments, it may happen that something goes wrong. The following problems are, at least, known:

• The O’Caml compiler is missing on the host computer. This is not a fatal problem since the compiler is part of most Linux distributions and available via MAC OS X port. Furthermore, it is easily available via http://pauillac.inria.fr/ocaml. Note that the O'Caml programming language can be installed on virtually any UNIX and MAC OS X system, and that it is perfectly possible to do this on a user account without superuser access.
• One of the PERL scripts which organize matrix element generation fails. This may be due to inconsistent shell return codes on some systems (IRIX, for instance). Since the OS world seems to converge towards Linux, this problem rarely occurs in reality. (Please notify the author if you really need a fix for a different OS. However, for WHIZARD 1 we basically consider Linux, MAC OS X and ALPHA as supported architectures.)
• A matrix element generator (O’Mega, CompHEP, MadGraph) fails. The reason is probably an attempt to generate a process that is not supported or which is physically impossible. ¹⁴ If this happens, there could be some trace in the temporary directory that has been created (in /tmp if TMPDIR has not been set). Note that you should remove this directory by hand if the program terminates abnormally during process creation.
• The Fortran 95 compiler fails. You may have found an actual WHIZARD bug, or you may have discovered a bug in the Fortran 95 compiler (this has happened more than once during the development of WHIZARD!). Make sure that you have the most recent compiler version and complain to the compiler vendor, if nothing else helps. Note however, that WHIZARD 1 explores no Fortran 03 features and most Fortran 95 are very mature nowadays).
• The linker fails due to incompatible object file formats. The reason could be that the WHIZARD configuration has chosen a Fortran 77 compiler different from the one used in generating the precompiled Fortran 77 libraries (PYTHIA, for instance). Set the correct compiler using the environment variable F77, reconfigure and recompile.
• The program crashes with a runtime error. This can be due to an incorrect input syntax, e.g., in a NAMELIST block. A frequent error is the specification of parameters that are misspelled or unknown to the current model. Another source of runtime errors can be incompatible external library formats (PYTHIA etc.). Otherwise, in particular if there is no reasonable error message, it may be a WHIZARD bug. Please report this to one of the email addresses below.
• Floating underflow messages after successful program execution are usually harmless.
• CompHEP cannot be compiled, so the executable x_comphep is missing when WHIZARD tries to generate matrix elements. This is likely due to missing C libraries needed by CompHEP. In particular, although WHIZARD will call CompHEP in command-line mode, the CompHEP executable needs the X libraries and include files installed on the host system. Note that CompHEP is just for self-testing purposes and is considered to be no longer maintained by us.

In any case, whether you could not solve a configuration or compilation problem or you found a workaround, please contact the author

kilian@physik.uni-siegen.de, ohl@physik.uni-wuerzburg.de, juergen.reuter@desy.de, christian.speckner@physik.uni-freiburg.de

so that the issue can be dealt with in a future release. (Please make sure that you are running the most recent version of WHIZARD, or at least look in the current CHANGES file whether a bug has been fixed in the meantime.)

4.3  Process selection

For running WHIZARD, the user can freely choose among the processes he has defined in the configuration file and compiled in the executable. The parameter process_id in the main input file (see below) is a string which contains the list of process tags, e.g.

  process_id = "zh, zz, zw"

If there is more than one process, integration of the cross sections will be performed in sequence. If a subsequent simulation step is requested, the events will be mixed according to their relative total cross sections. If fragmentation is enabled, one can even mix WHIZARD processes with processes generated by PYTHIA (see Sec. 4.8).

Since the integration grids are written to separate files, they can be reused later if events are to be simulated for the same or a different collection of processes, provided the input parameters and cuts are identical. This may save a lot of adaptation and integration time.

If the parameter process_id is left empty, WHIZARD will just print the list of processes, write a logfile containing the complete parameter list (including user settings), and exit.

4.4  Input data

When the whizard executable is started, it will use the parameters from the copy of the file whizard.in in the working directory. In this file, the physical input parameters and runtime switches are set. There are very many parameters, and all of them are optional. WHIZARD will try to insert sensible default values if the user has not specified a value. Of course, if no process tag is given, WHIZARD can do nothing, so it will just display the list of processes and exit.

The file whizard.in is written using the Fortran 95 NAMELIST conventions. It consists of several blocks marked by a keyword beginning with the & character and closed by a slash /. Each block contains a list of assignments of the form variable = value. Variables left out are assigned their default value. The assignments are separated by commas or newlines. Whitespace can be inserted anywhere. Comments can be inserted in this file; they begin with ! and extend up to the end of the line.

An example for the input file is shown in Fig. 4.4. In many cases, it suffices to fill the process_input block and leave all others empty. If polarization, beamstrahlung etc. are intended, the beam_input blocks must also be considered.

As an alternative to the NAMELIST format, WHIZARD supports the SUSY Les Houches Accord format for its input file. See Sec. 4.4.7 for details.

---

                      &process_input
                      process_id = "nnh"
                      sqrts = 500
                      luminosity = 100
                      /
                     
                      &integration_input
                      calls = 
                        1  10000
                        5  10000
                        2  20000
                      /
                     
                      &simulation_input
                      write_events = T
                      /
                     
                      &diagnostics_input /
                     
                      &parameter_input
                      Me = 0
                      Ms = 0
                      Mc = 0
                      MH = 115
                      wH = 0.3228E-02
                      /
                     
                      &beam_input 
                      /
                     
                      &beam_input 
                      /

Figure 4.4: Sample input file.

---

4.4.1  Files

The process_input block allows to specify a generic filename. This can be used to organize WHIZARD run data by giving specific filenames. For instance, to indicate the process energy in the filename for all output files, one may write

&process_input
 process_id = "nnh"
 sqrts = 800
 filename = "whizard_800"
/

Then, the output files will be whizard_800.out and whizard.800.nnh.out, etc. You do not need to write a separate file whizard_800.in, however: if a specific file is not found, the program will always look for the default one, in this case whizard.in.

Similarly, the working directory can be specified by the directory parameter. Of course, the first input file to be read must reside in the initial directory, otherwise it cannot be found (unless you define it on the command line, see Sec. 5.6).

The filename setting is overridden by specific filenames such as input_file, or e.g. in the block integration_input the variables read_phase_space_file or write_phase_space_file (cf. below). Such filenames will be interpreted relative to the chosen directory (if set), unless they begin with a slash, e.g. for read_model_file within the integration_input block ¹⁵,

  read_model_file = "/home/whizard/standard_model/whizard"

or with “./”, e.g.,

  read_model_file = "./whizard"

The latter form refers to the initial working directory. Note that in all cases the filename extension (here, .mdl) will be appended to the filename you specify.

It is possible to specify an additional input_file name in the process_input block. This input file will be read after the current one. New settings of parameters override old ones. In principle, any number of input files can be read consecutively and merged in this way. Reading input is finished when an input file does not exist, or neither directory nor input_file is changed in the last reading. Obviously, one should avoid loops in file reading.

All external file names used by WHIZARD can be similarly redefined. As a special rule, the filename for reading something (cuts, grids, events etc.) will automatically be set equal to the corresponding filename for writing, unless it is specified explicitly. The reason is that those files will be rewritten if the corresponding input file is not found, which could lead to accidentally overwriting files which is not intended.

4.4.2  The process_input block

These parameters should always be inspected.

4.4.3  The integration_input block

These parameters will not affect the cross section value (with the possible exception of the default cut parameters) but can improve or disprove the running time, the stability and accuracy of the result, and the efficiency of event generation.

The following parameters affect the details of the algorithm by which WHIZARD generates the phase space setup. Under normal circumstances, there should be no need to change them.

The following parameters affect the setup of kinematical cuts. Note that the default cuts are taken into account only if there is no appropriate entry in the user cut file whizard.cut1. However, even if such an entry exists, the cut parameters below affect the phase-space mappings, so adjusting them to order-of-magnitude may improve the convergence of the result.

4.4.4  The simulation_input block

These parameters control event generation, hadronization and output. If the luminosity (above) is left undefined or zero, one may still generate a fixed number of events by setting n_events nonzero. By default, the events are written to file in raw format and can be reused in that way. Writing events in another format (which uses up more disk space) must be requested explicitly. However, if this has not been done, one may use WHIZARD afterwards as a translator to any file format by reading pre-generated events and immediately writing them in another format:

&integration_input  read_grids = T /
&simulation_input   read_events = T  write_events = T /

WHIZARD also contains an experimental version of a k_T-ordered parton shower, written by Sebastian Schmidt. In the version it is implemented here, it is highly experimental and does work only for LHEF event formats (cf. below, this is the only one which allows for a varying number of particles in the event). A much more advanced version can be found in WHIZARD 2. For more details about the shower confer the WHIZARD 2 manual. The shower can be switched on with the shower flag.

4.4.5  The diagnostics_input block

These parameters do not affect the result, but the information displayed on screen and stored in files.

4.4.6  The parameter_input block

These are the physical constants used in evaluating the matrix elements. If this block is left empty, default values (see below) will be inserted.

Which constants are actually present, and which ones are dependent or derived, depends on the physical model. Consult the model files in conf/models for the corresponding parameters. For example, the constants relevant for the Standard Model SM_CKM are shown below together with their default values:

The dependent parameters (such as e, sinθ_W, V_cb etc.) will be shown in the output file whizard.out, but one should not try to reset them in the input file. The dependencies take into account the tree approximation only, such that gauge invariance is automatically respected. This may lead to unwanted side-effects in particular cases, e.g., the Higgs coupling to b quarks is proportional to the b mass, so setting m_b=0 removes all H→ bb subprocesses.

Setting fermion masses to zero will considerably speed up the matrix element evaluation since certain helicity combinations vanish identically (with O’Mega and MadGraph). In that case, cuts may be needed for a finite answer. However, if necessary, finite masses can be kept for all particles except the light quarks u and d.

Concerning quark masses in particular, they depend in fact on the chosen scale. The default values for the parameters include the strong coupling constant defined at the scale µ=M_Z, but the quark masses evaluated at the low scale µ=2 GeV. These default values are taken from the PDG 2000 compilation ¹⁶. For a consistent scheme, running quark masses should be inserted instead wherever this is of importance. For instance, the b mass is only about 2.9 GeV at µ=M_Z, thus affecting the Higgs coupling to b quarks.

Although they can be calculated from the other parameters in principle, the particle widths are retained as independent input parameters. Thus, the branching ratio of a particular resonance decay can be tuned by modifying the resonance width. When the Higgs mass is changed, the Higgs width has to be changed accordingly if Higgs decays are considered. The default values for the particle widths are given by the sum of the tree-level 1→ 2 decay channels, using the default masses and couplings as input. This choice is questionable as well. Since O’Mega(MadGraph and CompHEP, too) can only calculate cross sections at tree level, this is at least consistent; however, for each particular problem one should re-investigate the underlying assumption that 2-particle decays are dominant, and recalculate the widths accordingly. This cannot be done automatically by WHIZARD 1.

Finally, even if a parameter is defined and appears in the output, this does not necessarily mean that its value is used by all three programs CompHEP, MadGraph and O’Mega. For instance, the CKM matrix is not taken into account by the MadGraph version of the Standard Model.

For using CompHEP, there are two switches within the parameter_input block:

4.4.7  Input in SUSY Les Houches Accord format

The SUSY Les Houches Accord (SLHA) [10, 11] describes a format to pass MSSM/NMSSM physics data between computer programs. WHIZARD accepts its input file in SLHA format. This requires that either

• The first line of the input file begins with the string

  # SUSY Les Houches Accord
  

• or the parameter input_slha_format in the process_input block (see above) is .true. This flag can be set on the command line, or in an input file which is read before the SLHA data.

It is possible to

• read SLHA input after the usual input file: In the process_input block of the standard input file, set input_file to the name of the SLHA file (without the extension .in which each input file must have).
• embed standard WHIZARD input within the SLHA file: Enclose it in a block WHIZARD_INPUT. The text of this block is read by WHIZARD only and copied verbatim into an internal file which is read after the SLHA data. Note that, according to the SLHA standard, the first character of each line within this block has to be blank.

In this way, several input files in both formats may be read consecutively, which could be useful to set defaults and override some parameters computed by a SUSY spectrum calculator. However, the recommended procedure is to supply a standard whizard.in file where the process selection, collider information, and WHIZARD switches and runtime parameters are set, while containing all physical data in a separate file written in SLHA format which is read afterwards. The latter may be produced by an automatic spectrum and decay calculator.

In the output, if SLHA data have been used, the data will be repeated. WHIZARD reorders the blocks, where the blocks that actually have been used (which may include MODSEL, SMINPUTS, MINPAR, MASS, and mixing data) come first. The parameter slha_rewrite_input in the diagnostics_input block controls whether the input data are just copied into the output (preserving all comments, but not necessarily the ordering of blocks) or rewritten by WHIZARD. In the latter case, the relevant parameters in the output are made consistent between the standard WHIZARD format and the SLHA format, where applicable. Furthermore, for each process a block PROCESS will be appended which contains the integration results in a format which follows the SLHA conventions. The block CSINFO (cross section calculator information) is added, analogous to the SPINFO and DCINFO blocks described in the standard (see Fig. 4.5).

---

                 ! Begin SLHA data
                # SUSY Les Houches Accord Version 1.0
                Block MODSEL  # Select model
                    1    1    # mSugra
                Block SMINPUTS   # Standard Model inputs
                    1   1.279340000e+02  # alpha_em^(-1)(MZ)
                    2   1.166370000e-05  # G_Fermi
                    3   1.172000000e-01  # alpha_s(MZ)MSbar
                    4   9.118760000e+01  # MZ(pole)
                    5   4.250000000e+00  # Mb(mb)
                    6   1.743000000e+02  # Mtop(pole)
                    7   1.777000000e+00  # Mtau(pole)
                

⋮

                # WHIZARD: The following blocks have not been used:
                Block EXTPAR          # non-universal SUSY breaking parameters
                

⋮

                # WHIZARD: The following blocks have been added:
                BLOCK PROCESS   # Result of cross section calculation
                    1          2         # Number of incoming particles
                    2         11         #   Incoming particle: e
                    2        -11         #   Incoming particle: a-e
                   11          2         # Number of outgoing particles
                   12    1000024         #   Outgoing particle: ch1+
                   12   -1000024         #   Outgoing particle: a-ch1+
                   21    1.62815594E+02  # Cross section in fb
                   22    8.00542154E-02  # Cross section error in fb
                   23    1.56582525E+00  # Chi-squared value
                #
                BLOCK CSINFO   # Cross section calculator information
                    1    WHIZARD    # Calculator name
                    2    1.97       # Calculator version
                #
                 ! End SLHA data

Figure 4.5: Sample output in SLHA format

---

WHIZARD will use the MASS and DECAY information in the SLHA file to set the physical masses and widths of the particles it recognizes. Branching ratio information is not used. Note that the masses of some SM particles (Z, t, b, and τ) are defined in SMINPUTS. For completeness, they may also be set in the MASS block (where the W mass is located). If the values disagree, the MASS information takes precedence for WHIZARD’s needs. The widths of Z, W, and t may be set using the DECAY format. In all cases, default values (which may have been set in a previous input file) will be taken without warning where data are missing.

While the SLHA is intended for the MSSM and its extensions, it is possible to use this format to pass physical parameters for a SM calculation, too. If the SM is chosen when configuring WHIZARD, only the blocks MODSEL, SMINPUTS, MASS and DECAY information will be used. By convention, we define the model parameter (in the MODSEL block) to be negative for a non-SUSY model.

4.4.8  The beam_input blocks

The input file is finished by two blocks which describe the properties of the first and second incoming particle beam, respectively. It depends on the master flags in the process_input block above whether particular information within these blocks is actually used:

• If cm_frame is false, the beam_input blocks specify information on the beam energies and directions. If cm_frame is true (default), this information is discarded, the beam energies are given by the sqrts value, and the beam directions are assumed along the positive and negative z axis, respectively.
• If polarized_beams is true, the beam_input blocks specify information on the beam polarization. If polarized_beams is false (default), this information is discarded and the beams are assumed to be unpolarized.
• If structured_beams is true, the beam_input blocks specify information on the beam particle types, beam energy spectra and/or structure functions. If structured_beams is false (default), the beam particles are assumed to be the incoming particles of the hard process, and the partonic energies are set equal to the beam energies.

If all of the master flags are kept at their default values, the content of the beam_input blocks is completely discarded, so they may be left empty.

For decay processes, there are some obvious restrictions: There is only one “beam”, so the second block is ignored. While cm_frame and polarized_beams retain their meaning (i.e., one can simulate the decay of a moving and/or polarized particle), structured_beams is ignored, and no structure functions can be selected.

Energy and direction

If the master flag cm_frame is false, the beam energies and directions are set by

Polarization

If the master flag polarized_beams is true, the polarization is set for each beam by

For fermions, the default (helicity basis) polarization model describes longitudinally polarized beams with the specified fraction of left- and right-handed particles. In the case of massless vector bosons, this corresponds to left and right circular polarization along the beam axis. Massive vector bosons have three entries, where the middle entry specifies the fraction of longitudinal polarization. If the numbers sum up to less than 1, the remainder is filled with unpolarized particles. For instance, .5 0 and .75 .25 both stand for 50 % left-handed polarized electrons.

The quantization axis of a polarized particle is the momentum vector in the lab frame. For particles at rest, the quantization axis is the positive z axis.

The vector polarization model currently applies to fermions only. The three numbers denote the components of a three-vector P^→ with length less or equal to 1. The helicity density matrix for fermions is given by

where the polarization vector has components (P_x,P_y,P_z). For antifermions this is replaced by

so the y component changes sign.

While the z component of P^→ describes longitudinal polarization, the x and y components correspond to transversal polarization. Longitudinally polarized states can therefore be specified in both polarization models, while transversal polarization requires the vector model.

Spectra and structure functions

If the master flag structured_beams is true, WHIZARD will check if any of the structure functions defined below is applicable for this beam. For each beam, the beam particle is defined by either one of

If both entries are specified, the given PDG code takes precedence.

The further parameters specify which structure functions are used for each beam and fix structure function parameters where necessary. Structure functions are applied consecutively (Tab. 4.1), choosing a sensible chain of beam particle splittings: For instance, if the beam particle is an electron and the process is initiated by a quark, switching on both EPA and PDF will take the quark content of a photon which is emitted (in the Weizsäcker-Williams approximation) by the incoming electron. WHIZARD will not check whether the selected combination makes physically sense.

If it is desired to have the beam remnants in the event record, one should set the parameter keep_beam_remnants in the simulation_input block. The beam remnants will be photons (beamstrahlung, ISR), charged particles (EPA) or undefined hadron remnants (LHAPDF,PDF). A beam remnant is always massless, regardless which mass the particle has been assigned somewhere else, and it will have no transverse momentum. There will be exactly one beam remnant per structure function even though the number of photons radiated in beamstrahlung and ISR is not well-defined in reality. Beam remnants cannot be used for the built-in analysis since they are assigned no binary codes. Their main virtue is to make possible a reconstruction of individual splittings in the initial state in an external analysis.

Furthermore, one can set keep_initials. In that case, the beam particles and the partons which initiate the hard scattering are kept in the beam record. Together with the beam remnants, this allows for full mother-daughter relationships to be displayed. For fixed-energy beams, there are no beam remnants, and the initial partons are identical with the beam particles.

---

Structure function	in-particle	out-particle	user beam pol.
File events	any	= in	keep
User spectrum (Sec. 5.3)	any	any	keep/ignore
CIRCE (beamstrahlung)	e±	e±	keep
	e±	γ	ignore
CIRCE2 (beamstr.,γ spectra)	e±/γ	e±/γ	ignore
ISR (photon radiation)	charged particle	unchanged	keep
EPA (effective photon approximation)	charged particle	γ	ignore
EWA (effective W approximation)	weak doublet	weak doublet	ignore
LHAPDF (parton distribution functions)	p,p or γ	q or g	ignore
PDF (parton distribution functions)	p,p or γ	q or g	ignore
User structure function (Sec. 5.3)	any	any	keep/ignore

Table 4.1: Available structure functions. Each transforms an in-particle into a collinear out particle with a certain momentum fraction. They may be concatenated in the order given in the table. In the present implementation, the user-specified beam polarization is either kept unchanged for all events, or is overwritten.

---

Events from file

WHIZARD is able to read a sample of events from file and use it for integration and event generation. If both beams are read from file in this way and the filename is identical for both beams, the program expects a pair of real numbers in each line of the file. (Otherwise, it expects separate files with a single real number per line.) Each number corresponds to the energy carried by the beam particle. Empty lines and comment lines (marked by # or !) are skipped.

When the end of the beam event file is reached, a warning is issued and the file is reread from the beginning. Note that reusing the event sample in this way leads to systematic errors in the calculation. For instance, the fluctuations of the computed cross section between subsequent iterations can exceed the error estimate by orders of magnitude. In the final integration step, this could be visible as a χ² value which is anomalously large. To avoid this problem, it is recommended to use an event sample which is statistically distributed and contains a sufficient number of events so that no rewinding is necessary.

Parton distribution functions

There are two possibilities two access parton distribution functions inside hadrons: either via the LHAPDF library or via the PDFLIB inside the CERNLIB. Meanwhile we consider LHAPDF to be the default option, while we mention the PDFLIB setup just for backwards compatiblity and completeness. Distribution functions exist for all quarks and for the gluon. For LHAPDF, the parameters LHAPDF_on, PDF_scale (if PDF_running_scale is false) and LHAPDF_file for the PDF set are mandatory. For an overview over available sets and parameters confer the LHAPDF manual.

For PDFLIB inside the CERNLIB setup, the parameters are listed below (for a more detailed description, refer to the PDFLIB manual). For the numerical parameters, a zero value will cause the PDFLIB default value to be taken, so only PDF_on and PDF_scale (if PDF_running_scale is false) are mandatory.

Beamstrahlung (CIRCE1)

For beamstrahlung, the accelerator type should always be defined, and it is not recommended to use an energy different from the values for which valid CIRCE [6] parameterizations exist. There is the option to either use the generator mode of CIRCE, or to accept the distributions calculated by CIRCE as analytical spectra. For efficiency reasons, the generator mode is typically preferred. An example input file whizard.in.circe_isr for the default process ee_o can be found the conf directory.

In general, more information can be found in the CIRCE manual on the HepForge page:

http://projects.hepforge.org/whizard/circe1.pdf

Photon spectra (CIRCE2)

CIRCE2 is intended as a replacement for CIRCE1. While the previous implementation provided smooth spectra, the new version, which has to be used for photon collider spectra where no CIRCE1 version is available, provides histogrammed spectra which are read from a data file. These data files are partially contained in the WHIZARD distribution, in the subdirectory circe2-src/data. For more information on the CIRCE2 generator and possible updates on the spectra confer the manual/documented source code available from the HepForge page:

http://projects.hepforge.org/whizard/circe2.pdf

There is the option to either use the generator mode of CIRCE2, or to accept the distributions calculated by CIRCE2 as (histogrammed) spectra. For efficiency reasons, the generator mode is typically preferred.

The following parameters can be defined (identical for both beams):

If CIRCE2 is switched on, the CIRCE2_file parameter must always be given. The filename will be interpreted relative to the current working directory. The accelerator design should be specified if the data file contains more than one dataset. (The default matches the last dataset within the data file.)

It is important to note that polarization treatment is different from the standard WHIZARD conventions when CIRCE2 is used. The CIRCE2 data files introduce fixed polarization distributions of the beams which originate from realistic simulations. Therefore, the polarized_beams switch and the polarization settings should not be used together with CIRCE2 under normal circumstances. Instead, the appropriate data file with polarization should be taken and CIRCE2_polarized switched on. This will generate the appropriate polarized spectra, properly normalized.

On the other hand, if a polarization-averaged dataset is specified, CIRCE2_polarized should also be switched off.

Despite the fact that the polarization distributions are a fixed property of any given dataset, one can reweight the polarization contributions by hand. To this effect one can set the parameters polarized_beams and polarization. For instance, by switching on only one helicity combination, the properties of the corresponding cross section component can be studied. This feature should be used with care, since the results may be unphysical.

In the spectra mode (CIRCE2_generate=false), if no mapping is specified by the user, the program inserts a mapping at x=0 (power = 3) for photons, and a mapping at x=1 (power = 12) for electrons. This assumes that all particles originate from an incident e⁻ beam, where for electrons there is a strong singularity (in fact, a smeared-out delta function) at x=1. For a different setup, e.g., a positron beam, the mapping parameters have to be adjusted by hand. Note that without the x=1 mapping, the delta-function singularity for electrons will almost certainly be missed by the adaptive integration, and the result may stabilize for a wrong cross section.

A sample input file, whizard.in.photoncollider, for the assumed hard process γγ → e⁻ e⁺ can be found in the conf subdirectory.

ISR spectrum

The ISR spectrum for charged particles accounts for the leading-logarithmic effect of multiple photon emission. In the leading-logarithmic approximation, universal terms of order єⁿ with є=α/πlns/m² can be absorbed in a structure function. The leading singularity at x=1 can be calculated exactly, while the other terms of the structure function are known only to finite order.

The implementation in WHIZARD takes into account the LLA structure function ([7]), including terms up to order єⁿ, where n=0,1,2,3, depending on the value of ISR_LLA_order. The default value is n=3. The scale which enters the ISR spectrum is given by ISR_Q_max; the default is to take the collider c.m. energy.

Since α_QED is usually calculated from G_F as far as the hard scattering is concerned, it makes sense to reset it to 1/137 here. The incoming particle mass (e.g, m_e) must be defined here only if it has been set to zero in the parameter input block (e.g., to speed up the calculation).

You can find an example, whizard.in.circe_isr for the sample process ee_o in the conf subdirectory.

EPA spectrum

The EPA spectrum is the Weizsäcker-Williams approximation for a real photon radiated off the incoming beam. Here, Q_max and m_X must be supplied, the other parameters are optional. (If the produced system is massive, m_X should be set equal or less to the mass sum of all produced particles.)

Again, an example, whizard.in.epa for the sample process γγ → e⁻ e⁺ named aaee can be found in the conf subdirectory.

EWA spectrum

The EWA spectrum is equivalent of the Weizsäcker-Williams approximation for massive W and Z bosons radiated off the incoming beam. For details about the EWA and its implementation into WHIZARD confer [13].

Again, an example, whizard.in.ewa for the sample process W⁺W⁻ → ZZ named wwzz in the leptonic setup of a 3 TeV CLIC collider can be found in the conf subdirectory.

4.5  Cuts

WHIZARD allows for cuts in some standard kinematical variables. Cuts may be necessary in case the matrix element is infinite without them (which happens if massless particles are either exchanged or produced). In this case, WHIZARD will apply default cuts on the invariant mass of colored or charged particle pairs, on the energy of emitted photons or gluons, and on the momentum transfer to exchanged photons or gluons. The values of those cuts are controlled by parameters in the input file (see above).

---

! e-  e+  ->  e-  e+  gamma
! 16   8       1   2    4
process eeg
  cut Q of 10 within -99999 -1
  cut Q of 17 within -99999 -1
  cut M of  3 within 10 99999
  cut E of  4 within  5 99999
  cut PT of 4 within 19 99999
  cut THETA(DEG) of  4 1 within 5 180
  cut THETA(DEG) of  4 2 within 5 180

Figure 4.6: Cut configuration file

---

Alternatively, the user may specify his own set of cuts by supplying an entry in the file whizard.cut1. If such an entry exists, the program will ignore the default cuts, and the user is responsible for setting up his cuts in such a way that the cross section is still finite. The format of such an entry for the process e⁻e⁺→ e⁻ e⁺ γ is shown in Fig. 4.6. The contents of whizard.cut0, where the default cuts are logged by WHIZARD (if any), may be used as a template.

The entry is begun by the keyword process followed by the process tag (eeg in this case). If the same set of cuts should be applied to several processes, one may specify more than one tag per entry, e.g.

process uudd uuss uucc
    cut ...

The cut configuration file may contain any number of process entries. Comments are marked by the letter ! or #.

---

Code	Alternative code(s)	# Args	Description
-		0−2	No cut
M	Q	1	(Signed) invariant mass M=sgn(p2)√|p2|
LM	LQ	1	log10 |M|
MSQ	QSQ S T U	1	Squared invariant mass M2 = p2
E		1	Energy in the lab frame
LE		1	log10E
PT		1	Transverse momentum p⊥
LPT		1	log10 p⊥
PL		1	Longitudinal momentum pL
P		1	Absolute value of momentum |p→|
Y	RAP RAPIDITY	1	Rapidity y
ETA		1	Pseudorapidity η
DETA	DELTA-ETA	2	Pseudorapidity distance Δη
PH	PHI	1	Azimuthal angle φ (lab frame) in radians
PHD	PHID PHI(DEG)	1	Azimuthal angle φ (lab frame) in degrees
DPH	DPHI DELTA-PHI	2	Azimuthal distance Δφ (lab frame) in radians
DPHD	DPHID DELTA-PHI(DEG)	2	Azimuthal distance Δφ (lab frame) in degrees
AA	ANGLE-ABS TH-ABS THETA-ABS	1	Absolute polar angle θabs (lab frame) in radians. Reference axis is the z-axis.
AAD	ANGLE-ABS(DEG)
	TH-ABS(DEG) THETA-ABS(DEG)	1	Absolute polar angle θabs (lab frame) in degrees
CTA	COS(TH-ABS) COS(THETA-ABS)	1	cosθabs
A	ANGLE TH THETA	2	Relative polar angle θ (lab frame) in radians
AD	ANGLE(DEG)
	TH(DEG) THETA(DEG)	2	Relative polar angle θ (lab frame) in degrees
CT	COS(TH) COS(THETA)	2	cosθ
A*	ANGLE* TH* THETA*	2	Relative polar angle θ* (rest frame of part.#2) in radians
AD*	ANGLE*(DEG)
	TH*(DEG) THETA*(DEG)	2	Relative polar angle θ* (rest frame of part.#2) in degrees
CT*	COS(TH*) COS(THETA*)	2	cosθ*
DR	DELTA-R CONE	2	Distance in η-φ space, i.e. √Δη2 + Δφ2
LDR	LOG-DELTA-R LOG-CONE	2	log10√Δη2 + Δφ2
VBF	VBF_ETA HEMI HEMISPHERE	2	η1 × η2

Table 4.2: Cut and histogram code letters. Masses, energies and momenta are measured in GeV.

---

For a given process, each cut is indicated by the keyword cut, followed by a code letter, the keyword of, one or two binary codes, and cut window specifications. The code letters are listed in Table 4.2. The binary codes indicate the combination of particle momenta to which the cut should be applied. As the comment lines in Fig. 4.6 suggest, the outgoing particles are given codes 1,2,4,…. Two additional codes are assigned to the incoming particles, where the first incoming particle has the highest number. Momentum sums (particle combinations) are defined by adding binary codes, such that the code 5, for instance, stands for the sum of momenta of particles 1 and 4.

The initial momenta are always counted negative, such that (in Fig. 4.6) the number 17 stands for the difference of the momenta of particle 1 (outgoing) and 16 (incoming). The same applies to angles: The polar angle between 1 and 16 is not A of 1 16 but A of 1 8 since the initial particle 16 goes into the direction of −p(8).

By default, the inserted momenta are those of the elementary partonic process. Thus, one can assume momentum conservation. For instance, in a 2→ 2 scattering process, the invariant mass of 1+2=3 (the outgoing particles) is the same as the invariant mass of 4+8=12 (the incoming particles). In some circumstances, one may need the incoming beam momenta instead. To insert those in place of the incoming parton momenta, use a negative number. In the previous example, the invariant mass of −12 would be the invariant mass of the incoming beams which is equal to the nominal c.m. energy. Of course, there is no difference unless the process is a 2→ n scattering process with structure functions enabled.

Cut window specifications consist of the keyword within and two real numbers, optionally followed by additional specifications like

  cut PT of  4 within 0 150
  cut M  of  3 within 80 100 or 180 200 or 500 99999

The absence of an upper bound should be marked by a number outside the kinematical range.

4.6  Integration

4.6.1  Running WHIZARD

After the input parameters and cuts have been specified, the integration can be started by

make run

or, if we are in the working directory (results),

./whizard

After each iteration, the result will be displayed on screen and logged in the file whizard.XXX.out, where XXX is the process tag. An example is shown in Fig. 4.7. As it is expected, the integration converges towards a stable result, and at the same time the estimated reweighting efficiency increases (allowing for some fluctuations).

---

         ! WHIZARD 1.97 (May 31 2011)
         ! Reading process data from file whizard.in
         ! Wrote whizard.out
         !
         ! Process nnbb:
         !    e a-e -> nu_e a-nu_e   b a-b
         !   32  16 ->    1      2   4   8
        ....

         ! WHIZARD run for process nnbb:
         !=============================================================================
         ! It      Calls  Integral[fb]  Error[fb]   Err[%]    Acc  Eff[%]   Chi2 N[It]
         !-----------------------------------------------------------------------------
         ! Reading cut configuration data from file whizard.cut1
         ! No cut data found for process nnbb
         ! Preparing (fixed weights):   2 samples of      10000 calls ...
            1      20000  5.6761132E+01  1.07E+00    1.88    2.65*  1.45    1.52   2
         !-----------------------------------------------------------------------------
         ! Adapting (variable wgts.):  10 samples of      20000 calls ...
            2      20000  5.9620359E+01  8.15E-01    1.37    1.93*  6.58
            3      20000  5.8344071E+01  4.62E-01    0.79    1.12*  8.99
            4      20000  5.8023577E+01  4.17E-01    0.72    1.02*  9.04
            5      20000  5.7938879E+01  3.83E-01    0.66    0.94* 11.95
            6      20000  5.8152958E+01  3.73E-01    0.64    0.91* 11.85
            7      20000  5.8132996E+01  3.66E-01    0.63    0.89* 12.02
            8      20000  5.8252891E+01  3.50E-01    0.60    0.85* 13.61
            9      20000  5.8296069E+01  3.42E-01    0.59    0.83* 12.74
           10      20000  5.8391405E+01  3.32E-01    0.57    0.81* 15.93
           11      20000  5.8168510E+01  3.19E-01    0.55    0.78* 16.66
         !-----------------------------------------------------------------------------
         ! Integrating (fixed wgts.):   3 samples of      50000 calls ...
           12     150000  5.8232347E+01  1.17E-01    0.20    0.78  10.11    0.49   3
         !-----------------------------------------------------------------------------

Figure 4.7: Sample output of running WHIZARD.

---

The number of iterations and the number of calls per iterations are as given in the input file whizard.in (see Fig. 4.4). Adaptation proceeds in three passes which are separated by a horizontal line in the output:

In the first integration pass, the grids are adapted after each iteration, but the relative weight of the integration channels is kept fixed. In the example of Fig. 4.7, this pass consists of two iterations with 10,000 calls each. This pass is useful in particular for many-particle processes, where it takes some time before VAMP has found the physical region in each channel.

In the second pass, the grids within each phase space channel and the relative weights of the phase space channels are adapted after each iteration. The example of Fig. 4.7 shows how the integral converges, the error estimate improves, and the estimated reweighting efficiency increases, until some saturation is reached. Here, the sixth column is of main interest, which is given by the estimator for the relative error multiplied by the square root of the number of calls. One expects a number of the order of unity (or smaller) here for a well-adapted grid, independent of the number of calls. Each time this number improves, it is marked by an asterisk. The best grid will be used later for the final integration and event generation steps.

In the third pass, the relative weights are fixed again, and the remaining iterations are averaged for the final cross-section result.

4.6.2  The phase space configuration

The critical point for the integration is the phase space configuration. This is set up automatically by WHIZARD before starting the integration. When integrating a given process for the first time, a phase space configuration is generated and written to the file whizard.phx. In subsequent runs for the same process, this file is reused. If you would like to edit and modify this file or if you want to set up your own configuration, the new configuration should be stored in the file whizard.phs which takes precedence over whizard.phx and will not get overwritten. The automatically generated file whizard.phx will be deleted when the process list changes or, in any case, by a make clean command.

For each process, WHIZARD phase space consists of a set of channels which formally correspond to Feynman diagrams. In the automatically generated phase space configuration file whizard.phx, these channels are listed in order, using binary codes for the internal lines. The channels that have been generated can be drawn and viewed on screen or printed:

make channels

This generates a L^AT_EX file whizard-channels.tex which uses the FeynMF package for drawing Feynman graphs, and transforms it into a PostScript file whizard-channels.ps. These files are located in the results directory. Note that this is currently available only for the automatically generated phase space, not for a phase space setup supplied by the user. (Note further that the Feynman diagrams drawn are only sample diagrams for the corresponding phase space channels. Especially if the user edits phase space by removing channels from the file whizard.phx by hand, still the full matrix element is taken into account!)

Similarly, the bin distributions of the grids as they have been adapted during the WHIZARD run can be converted into histograms by

make grids

and viewed as the file whizard-grids.ps (note that this only works if the user has set the variable plot_grids_channels in the diagnostics_input block).

Again, it should be kept in mind that the choice of phase space channels does not influence the integration results (all Feynman graphs are included in the cross section calculation and event generation, but for the phase space setup only some of them are used — and some are even used twice). However, the accuracy and the speed of convergence depend on this choice. In many cases it will not be necessary to modify anything, but if you are not satisfied with the results of the default configuration, the following remarks may help.

• The channels are grouped in groves: Sets of channels that are equivalent in some sense and therefore should get equal weight. All channels within a grove are sampled with the same number of points per iteration which may increase or decrease during the adaptation procedure. To prevent channels from being dropped altogether, one may set min_calls_per_channel to a positive value, e.g. 100.
• When a Feynman graph is interpreted as a phase-space parameterization channel, the internal lines are considered either as resonant, log-enhanced, or off-shell. A resonant line either corresponds to a resonance which is kinematically allowed to go on-shell. A log-enhanced line is an (approximately) massless propagator, either in the s-channel (e.g., a photon which splits into massless quarks), or in the t channel (e.g., a Coulomb photon which is exchanged between the initial particles).
• The first criterion for selecting a channel is the number of off-shell lines. The default value is off_shell_lines=1. Since in dominant graphs most internal lines are either resonant or log-enhanced, this already covers many cases.
• By this criterion, some important graphs are missed, e.g., the Higgs radiation off a top quark in e⁺e⁻ annihiliation, which has two off-shell lines. The rule is that additional channels are selected with extra_off_shell_lines (default 1). For those, the log-enhanced lines count as off-shell. It turns out that this rule is sufficient for essentially all physically relevant cases. However, it may fail in more obscure cases (e.g., if non-resonant regions of phase space are selected by cuts), so increasing these parameters manually, which adds more Feynman diagram classes to the phase-space setup, is a possible option for the user.
• At high energies, multiperipheral graphs with particle exchange become important even though they may be subleading according to the former rules. This fact is accounted for by setting a threshold_mass_t value equal to 200 GeV, so W and Z bosons in the t-channel are considered light, and W/Z exchange graphs counted log-enhanced. For the s-channel, the threshold_mass value is 50 GeV, so W and Z are considered heavy. Of course, these values may be adjusted if the kinematics is completely different.
• If no phase space channels are found which satisfy all criteria, WHIZARD will give it another try with off_shell_lines increased.
• Five additional flags give finer control: single_off_shell_branchings, double_off _shell_branchings, single_off_shell_decays, double_off_shell_decays, massive _fsr. With the exception of the second flag, all are .true. by default, so switching them off will remove phase space channels.
• In particularly complicated situations you might wish to have full control and write your own configuration in whizard.phs, taking the pre-generated configuration whizard.phx as the starting point. To make this easier, for nonzero extra_off_shell_lines you will find extra channels in the output which are commented out. Instead of finding the binary coding of Feynman graphs yourself, you may just remove the comment characters for a inactive set of phase space channels.

The number of channels in the default configuration tends to be somewhat too large. Many channels turn out to be irrelevant and will be dropped during the self-adaptation procedure. If you can guess which ones those are, you may use the available parameters and options for manual intervention to optimize the channel selection. As an example, the effective Zee coupling is almost zero because of sin²θ_w∼ 1/4, so graphs with Z emission from electron lines are suppressed. You could remove such graphs, either manually, or by removing the Zee vertex from the whizard.mdl file before starting the program. (You have to be careful not to drop resonant Z→ ee decays, then.) However, in most cases, while this may save some running time for the adaptation, it is not worth the effort.

On the other hand, it may happen that important channels are missed by the default configuration. In this case, the restrictions should be weakened.

4.6.3  What if it does not converge?

Although WHIZARD is able to produce stable results with the default settings for a large variety of scattering processes, it is not possible to guarantee fast convergence of the cross section estimate in all circumstances. There is no simple rule for the achievable accuracy and event generation efficiency for a given process; the examples in Sec. 6 should provide some guideline. In any case the screen output of the adaptation procedure should be inspected, and questionable results should not be used for simulation.

Experience shows that in well-behaved 2→ 4 processes an accuracy (Err/Exp in the output) below 1 and an estimated reweighting efficiency above 10 % is reasonable. Due to the increased dimensionality of phase space, for 2→ 6 processes the numbers are typically worse by a factor 2 to 4.

If convergence is apparently not reached, i.e., there are sizable fluctuations in the result, this may be due to one of the following causes:

• The integral is actually infinite. One should check the cuts if they are sufficient for a finite result. Note that a singularity with a small coefficient is likely to produce an instability rather than an obvious divergence, if the number of calls per iteration is finite.

  A notorious example are processes where the final state contains µ⁺e⁻ν_eν_µ or µ⁺ūdν_µ: For zero muon mass, there is no problem, but if the muon mass is given the physical value, this final state can result from µ⁺µ⁻ production and muon decay. The muon width is zero, which is reasonable since the muon decay occurs outside the interaction point, but with cuts on visible particles only, real muon decay makes the integral divergent.

• The number of iterations is insufficient. For simple 2→ 3 processes, one reaches a stable result after 2−5 iterations, but in complicated cases (many particles with photon emission with weak cuts, etc.) 20−30 iterations may be needed. WHIZARD will make a guess, but this may be not sufficient in some cases.
• The number of calls per iteration is insufficient. With increasing dimension and complexity of the phase space, this number should be increased as well.
• The phase space parameterization is inadequate (see above). One may try to include or remove channels by modifying parameters such as the number of off-shell lines allowed in generating parameterizations.
• The parameterization is adequate in principle, but some subdominant (but still important) channels got dropped during the adaptation procedure. This may be prevented by setting min_calls_per_channel. The channel weight history can be inspected in the process-specific output file.
• The matrix element is numerically unstable. This is probably not an issue for O’Mega or MadGraph matrix elements with realistic kinematics, but cannot be excluded. CompHEP matrix elements, which result from squared amplitudes, are susceptible for numerical instabilities which occur for extreme kinematics. Actually, if the compiler supports it, one may switch to quadruple precision (configure with --enable-quadruple) in order to ameliorate this problem.
• One should be careful to avoid the effects of QED gauge invariance violation, which can be large if unstable particles and photons are around (e.g., single W production in e⁺e⁻ collisions). With CompHEP matrix elements, the gwidth switch should be set in such cases, while O’Mega matrix elements should be generated with the ‘f’ option. This enables the overall factor scheme, where the resonant contributions to the matrix element are evaluated for zero width of charged particles, and the Breit-Wigner correction factor is applied to the complete matrix element at the end. Another gauge-invariant prescription (the complex-mass scheme) is foreseen to be implemented in a future version of O’Mega within the context of WHIZARD 2.

4.6.4  Output files

All results are logged in a generic file whizard.out and in process-specific files whizard.XXX.out, where XXX is the process tag. The output files are rewritten after each iteration. Depending on the settings of the diagnostics parameters, they will contain additional information: The contents of the input file are repeated, including also the default values which have been inserted by the program. Furthermore, it shows the detailed history of the individual integration channels and of their relative weights.

After each iteration, the current VAMP grids are written to file as well: whizard.XXX.grc for the current grid and whizard.XXX.grb for the best grid. These may be reused in another run.

4.7  Event generation

4.7.1  Unweighted events

Event generation is activated in the input file, either by specifying a nonzero luminosity

&process_input
...
luminosity = 100
/

or by a nonzero number of events

&simulation_input
n_events = 10000
/

In the former case, the actual number of generated events depends on the calculated cross section. If both numbers are given, the one leading to the larger event sample will be used. It is also possible to specify n_calls, which is the number of matrix element calls. However, this is more useful for weighted event samples (see below).

The simulation makes sense only if WHIZARD could find a stable result for the cross section. This should be judged from inspecting the results displayed on screen and in the log file whizard.out.

With event generation activated, WHIZARD will not stop after the integration step, but start event generation. By default, WHIZARD will generate unweighted events, which result from reweighting the events produced by the VAMP sampling algorithm. In the event generation step, importance sampling is used, such that the events are random and have the proper distribution as required for a realistic simulation.

WHIZARD is able to write events to file in various formats (note that the manual of WHIZARD 2, http://projects.hepforge.org/whizard/manual.pdf), does contain examples for all the different event formats presented here, though the debug format does look different there):

• The file whizard.evx contains the events in a binary machine-dependent format. By setting read_events=T they can be re-read, bypassing the time-consuming event generation. Writing this file is suppressed if write_events_raw=F.
• The file whizard.evt (or, alternatively, whizard.pythia or whizard.stdhep, depending on the format) is written if write_events=T. The event file format is controlled by the parameter write_events_format. The supported formats are:
  • 0 Verbose format (ASCII), useful mainly for debugging.
  • 1 HEPEVT format (ASCII). The format corresponds to the HEPEVT standard, such as if the contents of a HEPEVT common block were written in their natural order (Fig. 4.8, 4.9). As an extension, the helicity (−1, 0, or 1) is written for each particle in the same line after JDAHEP. For unfragmented events, the first line will contain further information: the number of beam remnants and outgoing partons listed separately, and the event weight.
  • 2 SHORT format (ASCII). This is similar to HEPEVT, but unnecessary entries are suppressed (Fig. 4.10).
  • 3 LONG format (ASCII). The same as SHORT, but there is an additional line per record which contains the value of the integrand and the integrand ratio (Fig. 4.10). In a normal run (recalculate=F), the ratio is unity (see Sec.5.5).
  • 4 ATHENA format (ASCII). This is the input format used by the ATHENA software (ATLAS). A slight variation of the HEPEVT format above. (cf. Fig. 4.12)
  • 5 Les Houches Accord format LHA (ASCII). This format implements the Les Houches Accord as in the original MadGraph implementation for transferring parton-level events to hadronization packages. To use it, you must set fragment=T, and you may set fragmentation_method=0 since any built-in fragmentation is irrelevant for the output. (cf. Fig. 4.13)
  • 6 Les Houches Accord format LHEF (ASCII). The new version containing XML headers [18], cf. Fig. 4.14. For more details cf. the WHIZARD 2 manual,
    http://projects.hepforge.org/whizard/manual.pdf.
  • 11/12/13 PYLIST formats 1,2,3 (ASCII). This is useful in particular for fully fragmented events, but can be used as well if fragmentation is turned off (see the PYTHIA manual [1] for details).¹⁷
  • 20 STDHEP format for writing the HEPEVT common block (binary, machine-independent)
  • 21 STDHEP_UP format for writing the HEPRUP/HEPEUP common blocks (binary, machine-independent)

Note that the formats that require external libraries (PYTHIA, STDHEP) do not work if the corresponding libraries have not been linked to the main program.

Within each event, the particles will appear in the order determined by the process definition /footnotefor O’Mega and MadGraph; CompHEP may reorder them.. Hadronization packages such as JETSET require a particular order of colored particles to properly assign showers, be may necessary for the user to rearrange them before feeding them into the hadronization step. In general, if color-flow information is absent ¹⁸, one may set guess_color_flow to tell the program that the particles are to be combined as singlets in the order in which they appear in the event record.

Reweighting Monte-Carlo events requires the knowledge of the highest weight. This can never be known with certainty in advance, so there may be some fraction of events with weight greater than one. By default, WHIZARD will sum up the excess weight of those events and show the induced error in the total cross section, which should be small compared to the statistical error of event generation. If desired, the excess events can be further analyzed by WHIZARD’s own analysis capabilities (see below), to make sure that they have no significant impact on final results. By setting the safety_factor in the simulation_input block to a number greater than one, the amount and weight of excess events (if any) will be reduced. However, the time needed for event generation will be increased by the same factor.

WHIZARD will use the adaptation step to determine the highest weight beforehand. This requires sufficient statistics in the adaptation procedure. If many events have to be generated, the default choice for the number of calls per iteration may be too low, and too many excess events are produced. As a rule of thumb, the number of calls per adaptation iteration should be of the order of the number of unweighted events divided by the estimated efficiency.

Alternatively, a “test” run can be made by specifying a nonzero value of n_events_warmup in the input file. For instance, one may generate 200,000 (weighted) events to determine the highest weight before generating the actual sample to be used for the reweighting procedure:

&simulation_input
n_events = 10000
n_events_warmup = 200000
/

Fortunately, this is not necessary in most cases since the pre-determined highest weight is usually quite accurate.

---

  integer, parameter :: nmxhep = 4000
  integer :: nevhep, nhep
  integer, dimension(nmxhep) :: isthep, idhep
  integer, dimension(2, nmxhep) :: jmohep, jdahep
  real(kind=double), dimension(5, nmxhep) :: phep
  real(kind=double), dimension(4, nmxhep) :: vhep
  common /HEPEVT/ nevhep, nhep, isthep, idhep, &
       & jmohep, jdahep, phep, vhep

Figure 4.8: Definition of the HEPEVT common block in Fortran 95 format

---

---

 write(u, 11) nhep
 do i=1, nhep
    write (u,13) isthep(I), idhep(I), &
         & (jmohep(j,i), j=1,2), (jdahep(j,i), j=1,2)
    write (u,12) (phep(j,i), j=1,5)
    write (u,12) (vhep(j,i), j=1,min(5,ubound(vhep,1)))
 end do

 11 format(10I5)
 12 format(10F17.10)
 13 format(I9,I9,5I5)

Figure 4.9: Output routine which would produce the event file whizard.evt from the contents of a HEPEVT common block. u is an output unit.

---

---

write(u,15) evt%n_out + evt%n_beam_remnants, &
     &      evt%n_out, evt%n_beam_remnants, weight
do i=1, evt%n_out
   p = array(particle_momentum(evt%prt(i)))
   write (u,13) particle_code(evt%prt(i)), evt%hel_out(i)
   write (u,12) p(1:3), p(0), particle_mass(evt%prt(i))
end do
12 format(10F17.10)
13 format(I9,I9,5I5)
15 format(3I5,1X,F17.10)

Figure 4.10: Output routine for the SHORT event file format. Note that the number of beam remnants and the event weight are not part of the HEPEVT standard. (We here and in the sequel do not include the code for the beam remnants and/or initials.)

---

---

write(u,15) evt%n_out + evt%n_beam_remnants, &
     &      evt%n_out, evt%n_beam_remnants, weight
do i=1, evt%n_out
   p = array(particle_momentum(evt%prt(i)))
   write(u,13) particle_code(evt%prt(i)), evt%hel_out(i)
   write(u,12) p(1:3), p(0), particle_mass(evt%prt(i))
end do
write(u, 14) evt%mc_function_value, evt%mc_function_ratio
12 format(10F17.10)
13 format(I9,I9,5I5)
14 format(ES17.10,F17.10)
15 format(3I5,1X,F17.10)

Figure 4.11: Output routine for the LONG event file format. Note that the number of beam remnants and the event weight are not part of the HEPEVT standard.

---

---

write(u,11) evt%count, event_n_entries (evt)
k = 0
do i=1, evt%n_out
   k = k + 1
   p = array(particle_momentum(evt%prt(i)))
   write (u,13) k, 1, particle_code(evt%prt(i)), &
        & evt%n_in+1, 2*evt%n_in, (0, j=1,2)
   write (u,12) p(1:3), p(0), particle_mass(evt%prt(i))
   write (u,12) (0._default, j=1,4)
end do
11 format(10I5)
12 format(10F17.10)
13 format(I9,I9,5I5)

Figure 4.12: Output routine for the ATHENA event file format. (If the initial particles and/or the beam remnants are to be kept, these are written before the outgoing particles.)

---

---

write (u,16) nup, idprup, xwgtup, scalup, aqedup, aqcdup
write (u,17) idup(:nup)
write (u,17) mothup(1,:nup)
write (u,17) mothup(2,:nup)
write (u,17) icolup(1,:nup)
write (u,17) icolup(2,:nup)
write (u,17) istup(:nup)
write (u,17) ispinup(:nup)
do i=1, nup
   write (u,18) i, pup((/4,1,2,3/),i)
end do
16 format(2(1x,I5),1x,F17.10,3(1x,F13.6))
17 format(500(1x,I5))
18 format(1x,I5,4(1x,F17.10))

Figure 4.13: Output routine for the (old) Les Houches Accord (LHA) event file format.

---

---

  subroutine lhef_write_init (unit, process_id)
    integer, intent(in) :: unit
    integer :: i
    character(*), dimension(:), intent(in) :: process_id
    write (unit, "(A)")  '<LesHouchesEvents version="1.0">'
    write (unit, "(A)")  '<header>'
    write (unit, "(2x,A)")  '<generator>'
    write (unit, "(4x,A)")  &
         '<name version="<<Version>>" date="<<Date>>">' // &
         'WHIZARD</name>'
    do i = 1, size (process_id)
       write (unit, "(4x,A)") &
            '<process_id>' // trim (adjustl (process_id(i))) &
            // '</process_id>'
    end do
    write (unit, "(2x,A)")  '</generator>'
    write (unit, "(A)")  '</header>'
  end subroutine lhef_write_init

  subroutine lhef_write_heprup (unit, process_id)
    character(*), dimension(:), intent(in) :: process_id
    integer, intent(in) :: unit
    integer :: i
    write (unit, "(A)") '<init>'
    write (unit, "(2(2x,I10),2(2x,G22.15E3),8x,6(2x,I10))") &
         idbmup(1:2), ebmup(1:2), &
         pdfgup(1:2), pdfsup(1:2), &
         idwtup, nprup
    do i = 1, size (process_id)
       write (unit, "(3(2x,G22.15E3),8x,1(2x,I10))") &
            xsecup(i), xerrup(i), xmaxup(i), lprup(i)
    end do
    write (unit, "(A)") '</init>'
  end subroutine lhef_write_heprup

  subroutine lhef_write_hepeup (unit, pro_id)
    integer, intent(in) :: unit, pro_id
    integer :: i
    idprup = pro_id
    write (unit, "(A)") '<event>'
    write (unit, "(2(2x,I10),2(2x,G22.15E3),8x,2(2x,G22.15E3))") &
         nup, idprup, xwgtup, scalup, aqedup, aqcdup
    do i = 1, nup
       write (unit, "(6(2x,I10),8x,3(2x,G22.15E3),8x,2(2x,G22.15E3),32x,2(2x,G22.15E3))") &
            idup(i), istup(i), &
            mothup(1:2,i), icolup(1:2,i), &
            pup(1:5,i), vtimup(i), spinup(i)
    end do
    write (unit, "(A)") '</event>'
  end subroutine lhef_write_hepeup

  subroutine lhef_write_final (unit)
    integer, intent(in) :: unit
    write (unit, "(A)")  '</LesHouchesEvents>'
  end subroutine lhef_write_final

Figure 4.14: Output routine for the (new) Les Houches Accord (LHEF) event file format.

---

4.7.2  Weighted events

The default setup will generate unweighted events, which is the proper way of simulating actual physics runs. The cross section error of the simulation run is equal to the square root of the number of events, as it should be.

This requires an unweighting step where, depending on the ratio of the current event’s weight and the maximum weight, a large fraction of the generated event sample is thrown away. Clearly, information is lost in this step. To keep all available information, one can use weighted events instead:

&simulation_input
unweighted = F
...
/

The number of weighted events is calculated such that (approximately) the same unweighted event sample would result from unweighting them afterwards, as it would from directly generating unweighting events. If luminosity or n_events do not require a higher number, the total number of weighted events generated is approximately equal to n_calls.

If there is more than one process, the luminosity is distributed among them proportional to the individual cross sections. In contrast to unweighted generation, the events are ordered, such that all events belonging to the first process are generated first, then the next one, and so on. The exact number of weighted events for a given process depends on both the cross section and the efficiency. For a process with low efficiency, a larger number of weighted events is needed for the same cross section than for a high-efficiency process. (In unweighted event samples, the individual event numbers depend only on the cross sections.)

These events can also be written to file and analyzed by an external program. For each event, the event weight is a number between zero and one (if no excess events are present), and it is printed as the last number in the first line of the output. Hence, in this definition the weight is the probability with which the event would have been kept in the sample if an unweighting step had taken place. As a result, differential (binned) distributions should be identical for weighted and unweighted event samples, apart from fluctuations.

Unfortunately, the HEPEVT standard does not allow for event weights. Therefore, the PYLIST and STDHEP output formats which read out the HEPEVT common block cannot be used for weighted events. By contrast, event weights can be transferred using the HEPEUP common block, so the STDHEP_UP format does work for weighted events. Keep in mind that weighted events are ordered, so files containing those must not be truncated or extended.

4.7.3  Built-in analysis

WHIZARD has capabilities to analyze the events without reference to external packages. To this end, a configuration file whizard.cut5 should be set up which has a similar format as whizard.cut1. For instance, in Fig. 4.15 a configuration file is shown which tells WHIZARD to analyze the generated event sample in three different ways: to calculate the cross section within a window 80<M(e⁺e⁻)<100 GeV, then within 180<M(e⁺e⁻)<200 GeV with an additional cut on p_⊥(γ), and finally to generate a histogram of p_⊥(e⁻) with a cut on the photon energy.

---

! e-  e+  ->  e-  e+  gamma
! 16   8       1   2    4
process eeg
  cut M of 3 within  80 100
and
  cut M  of 3 within 180 200
  cut PT of 4 within 100 99999
and
  cut E of 4 within 0 100
  histogram PT of 1 within 0 500

Figure 4.15: Analysis configuration file.

---

The histogram specification may take either one of the form

histogram PT of 1 within 0 500
histogram PT of 1 within 0 500 nbin 50
histogram PT of 1 within 0 500 step 10

Instead of the number of bins (20 by default), the bin width can be specified by the keyword step.

The cross section results are displayed on screen. In addition, they are written together with the histogram data to the file whizard-plots.dat. If show_excess=T has been set in the input file, the summary and the histograms will contain also information about any events with weight greater than one, such that their properties can be examined.

Clearly, for a quick analysis it is not necessary to simulate proper unweighted events. Therefore, one may set unweighted=F in order to analyze a sample of weighted events, which saves a lot of execution time.

The contents of whizard-plots.dat can be visualized directly by means of the gamelan graphical analysis package which is contained in the distribution. This requires the existence of the MetaPost program, which is part of standard L^AT_EX distributions. Plots can then be made by typing

make plots

This will result in a Postscript file whizard-plots.ps which can be viewed and printed. Note that this works only if you have generated events, otherwise WHIZARD might fail with the message

No rule to make target plots.

The event sample can further be analyzed in terms of color flow, helicity and flavor content. Example

! e-  e+  ->  q  qbar  b  bbar
! 16   8      1   2    4   8
process qqbb
  cut M of 3 within 85 95
  color 0 0 4 0 6 0
  helicity  1 -1 -1 1
  helicity -1  1 -1 1
  flavor    2 -2 5 -5
  histogram M of 12 within 0 500

This would plot a histogram of the invariant mass of the bb pair with the additional requirement that the other quark pair must be uū, that the b (anti)quarks must be left- resp. right-handed, and that the uū and bb pairs are each in a color-singlet state.

Color flow is defined in a natural way: The particles are numbered consecutively beginning from the initial ones (do not confuse this with the binary codes used for specifying cuts and histograms). For an incoming quark or gluon, the index of the particle its color flows to is given, while for an outgoing quark or gluon it is the index of the particle its color comes from. (This differs from WHIZARD’s internal convention where incoming and outgoing particles are equally treated, but it is easier to use.) For colorless particles and for antiquarks, the entry is zero. So, in the above example, the q color comes from the q, and the b color from the b.

Helicities and flavors are indicated as usual: Fermions and massless vector bosons have helicity ± 1, massive vector bosons can come with helicity 0 in addition, and for scalars the helicity is always zero. Flavor is denoted by standard PDG codes.

The order should be as in the above example: first cut specifications (the event must pass all cuts), then color, helicity and flavor specifications (the event must match any of the given helicity combinations and any of the given flavor combinations), finally histogram specifications. Any of these may be absent; an empty entry will just display the total cross section again. More analysis configurations can be added for the same process, they must be separated by the keyword and.

4.8  Fragmentation

The Les Houches standard [9] for transferring events together with weight and color information allows for a consistent interface of partonic event generators and programs that account for showering, fragmentation, and hadronization. In order to use such an external program, WHIZARD has to write events to file in a format that supports this information. Currently, there is the LHA format (code: 5) that writes, for each event, the contents of the HEPEUP block in ASCII format, and the STDHEP_UP format (code: 21) that writes this information into a portable binary STDHEP file. The first record of this file will contain the HEPRUP common block; then, for each event, the HEPEUP common block is written.

Alternatively, WHIZARD provides an internal interface to routines for fragmentation, decays and hadronization. For each event, WHIZARD calls the corresponding external program that fills the HEPEVT common block and writes this to file. To activate this, the parameter fragment in the simulation_input block must be set to true. Of course, the corresponding library has to be properly linked. Only unweighted events can be fragmented. (This is because the HEPEVT standard does not allow for the event weight to be transferred.)

Once fragment is set, various options are available, depending on the value of the parameter fragmentation_method:

1. No fragmentation within WHIZARD (for completeness).
2. JETSET fragmentation. The final state is fragmented and hadronized by a call to the PYEXEC fragmentation routine inside PYTHIA, i.e., by the former JETSET program. If available, color flow information is used. Otherwise, i.e. for CompHEP matrix elements and O’Mega without the ’c’ option, the color connection properties of the final state are guessed by JETSET from the particle ordering, which must therefore carefully be chosen by the user when defining the process. Note that even for processes like e⁺e⁻ → jj where O’Mega without the option ’c’ gives the correct color factors, it only gives the correct color flows for (JETSET) fragmentation when using the ’c’ option.
3. PYTHIA fragmentation. The final state is fragmented and hadronized by making the process generated by WHIZARD an external process within PYTHIA. If possible and depending on the PYTHIA parameter settings, beam remnants are constructed, multiple interactions and pile-up events are added, and more.
4. This parameter value is intended to transfer control to fragmentation routines written by the user and included in the user module. Of course, such code can make use of the HEPRUP and HEPEUP common block contents. For additional flexibility, there is also a parameter user_fragmentation_mode (integer) that allows for selecting different versions of user fragmentation.

Full color information is available when using the ’c’ option for O’Mega matrix elements (or for MadGraph matrix elements). However, in the other cases, for PYTHIA fragmentation WHIZARD will provide missing color information if can be uniquely determined. Otherwise, it tries to guess color connections from the particle ordering in the event record (just as JETSET does). If color connection information cannot be determined, colored states cannot be fragmented in this way. If fragmentation is nevertheless desired (even though the results may not be accurate), one has to revert to JETSET fragmentation.

No effort is done to streamline WHIZARD and PYTHIA parameters, so changing a value in the WHIZARD input file will have no effect on the fragmentation step. However, PYTHIA parameters can be explicitly set by an interface to PYTHIA’s pygive routine. For instance, the line

  pythia_parameters = "MSTP(41)=0; MSTP(81)=0"

in the simulation_input block disables resonance decays and multiple interactions within PYTHIA. Similarly, the Higgs mass may be set by a call like

  pythia_parameters = "PMAS(25,1)=115"

The detailed list of parameters available is described in the PYTHIA manual [1].

The internal analysis routines described in Sec. 4.7.3 always refer to the event before fragmentation. To analyze the fragmented event, the user must use an external package. Note that WHIZARD 2 meanwhile allows for analyzing fragmented events with the internal analysis tools, however.

4.8.1  Parton showers and matching

In WHIZARD there is a version of a k_T-ordered parton shower included written by Sebastian Schmidt. This is, however, a highly experimental version mainly intended for testing purposes. It is quite similar to the PYTHIA parton shower [1], but does contain only very basic tunes to LEP Run 2 jet observables. Furthermore, it is only a shower for final state radiation (FSR); there is no implementation of an initial state shower in WHIZARD 1. The internal shower has been improved a lot and augmented by a different (so-called analytical) parton shower within WHIZARD 2. For more details about the shower algorithm confer the WHIZARD 2 manual http://projects.hepforge.org/whizard/manual.pdf. The internal parton shower can be switched on for event generation in the simulation_input block. Note that the only event format which supports showered events from the WHIZARD-internal parton shower is the (new) Les Houches event format (LHEF, event format 6) as this is the only format that allows for a variable number of particles.

There is, of course, also the possibility to use showering (and hadronization) using the PYTHIA interface. In the default settings, PYTHIA events are constructed with parton showers developing along colored and charged partons. In WHIZARD, fixed-order matrix elements are evaluated. The matching of both is nontrivial. For this reason, PYTHIA parton showers are disabled in the current WHIZARD version. In general, the availability of exact higher-order matrix elements (albeit tree-level) within WHIZARD makes matching without double-counting a delicate task. This problem of potential double counting has been (partially) solved in the framework of so called matching algorithms by either vetoing shower events entering hard matrix elements regimes or reweighting with the corresponding Sudakov factors according to the shower history. None of these algorithms have been implemented within WHIZARD 1, but they are (at least partially) available inside WHIZARD 2.

There is also an experimental Herwig interface, but we do not discuss the details here.

4.8.2  Mixing WHIZARD and PYTHIA events

Since WHIZARD provides an external process to PYTHIA when PYTHIA fragmentation is enabled, it is perfectly possible to mix PYTHIA and WHIZARD events. To this end, one should specify PYTHIA processes in the simulation_input block. For instance,

  pythia_processes = "22, 24"

will simulate e⁺e⁻→ ZH and ZZ in addition to the WHIZARD-generated processes. There are some caveats, however:

• PYTHIA processes are generated only if PYTHIA fragmentation is enabled. Consequently, the PYTHIA events are mixed with the fully fragmented WHIZARD events and optionally written to file. There is no way to use the built-in analysis for those events.
• Since parton showers (ISR and FSR) are disabled globally, this applies to the PYTHIA processes as well. However, one could switch the showering on by modifying the inputs.
• When the number of events to generate is derived from the luminosity input parameter, it is assumed that no PYTHIA events are present. To obtain the actual luminosity corresponding to an event sample in the presence of PYTHIA events, this number has to be multiplied by the fraction of WHIZARD events in this sample. The parameter n_events always refers to the total number of events. When show_pythia_statistics is true, a summary of the generated channels and cross sections is printed.
• One should be careful not to double-count processes, except if this is actually desired (e.g., for comparing WHIZARD and PYTHIA event samples for related processes).

4.8.3  PYTHIA event display

When fragmentation is enabled, the conventions for writing events to file are changed. The verbose format is then not available since some information (used internally by WHIZARD) is lost after fragmentation. The format codes 1 (HEPEVT), 2 (SHORT), 3 (LONG), and 20/21 (STDHEP) are available. In addition, there are three new format codes (11, 12, 13) which allow for displaying events using PYTHIA’s PYLIST call. The respective argument given to PYLIST is 1, 2, or 3, corresponding to three different verbosity levels. (For automatic processing, the more condensed HEPEVT or STDHEP formats are better suited.) Note again that these output formats might not work properly if you use a CERN library compiled with a different compiler than the WHIZARD core program, because in this case the output units of the two compiler libraries are independent of each other. The PYTHIA events then just appear in a file fort.<unit>, where <unit> is the current output unit.

4.9  Controlling WHIZARD runs

For a new process, it is not easy to predict how long it will take to reach a sufficiently precise result or to generate a sufficiently large event sample. The achievable efficiency and accuracy and the overall running time strongly depend on the complexity of the process, the singularity structure of the matrix element and how well WHIZARD is able to map it, the presence or absence of structure functions, and, clearly, the evaluation time of the matrix element for a single event.

Ideally, computer time and space is abundant, one calcaulates cross sections to the highest possible accuracy, and infinitely large event samples are generated. In practice, one does not want to wait too long for useable results, and one is limited by CPU time and space restrictions. Priorities have to be defined. To meet these needs, various switches and levers are available that control WHIZARD runs. All of them are set in the main input file and are listed in the tables of Sec. 4.4.

4.9.1  Grid adaptation

Parameters that control the grid adaptation and integration passes of WHIZARD are collected in the integration_input block in the input file, listed in Sec. 4.4.3

Obviously, the array calls is the principal place for improving the accuracy and efficiency of the adaptation and integration runs. In no case more iterations or calls than specified in calls will be executed, where WHIZARD uses default values for entries left empty by the user.

Increasing the number of adaptation iterations normally improves (decreases) the constant α in the error estimate for the total cross section,

until a limit is reached which depends on the matrix element singularity structure and on the choice of phase space integration channels. The estimate for α is printed as the accuracy in the screen output table. Increasing the number of points per iteration tends to decrease the number of iterations needed to reach this limit. In particular, this improves the precision in estimating α and reduces inevitable fluctuations around the ideal grid configuration. Simultaneously the reweighting efficiency є increases, which is the ratio of the average to the maximum event weight. Note that fluctuations in the є estimate are significantly stronger than for α.

If the program has been killed for any reason, it can be restarted with read_grids set, so adaptation is resumed after the last complete iteration. This mechanism can also be used to increase the number of iterations afterwards for a job that has already been completed.

However, one wants to avoid a unnecessary large number of iterations and calls. Adaptation iterations are skipped if any one of the following parameters is set:

accuracy_goal:

If this number is greater than zero, adaptation will stop when the accuracy estimate of the last iteration is smaller.

efficiency_goal:

If this number is greater than zero, adaptation will stop when the efficiency estimate of the last iteration is higher. Note that the efficiency is measured in percent.

time_limit_adaptation:

If this number is greater than zero, adaptation will stop when the real time after the last iteration is later than the specified time limit, measured in minutes from the program start.

In any case, the program will not stop but the third pass (integration) will be executed after the adaptation pass, where the number of integration steps and calls per iteration are taken from the corresponding calls entry.

4.9.2  Simulation

The reweighting efficiency in the final event generation pass is determined by the grid quality, measured by the efficiency є. Since the precision in estimating є is not very high, it may be appropriate to insert a safety_factor f>1 to avoid excess events with weight >1. The time needed for generating a fixed number of unweighted events is then proportional to 1/(fє). As a rule, investing more time in grid adaptation results in a better є value, and thus less time for event generation, up to a certain limit.

The number of unweighted events is specified by either setting the parameter n_events in the simulation_input block, or by a luminosity value in the process_input block. The number that leads to a larger event sample will take priority.

If the program is killed during event generation, the event sample is incomplete. The switch catch_signals can be used to catch terminate signals if the OS and compiler supports it, so it is guaranteed the the event files are properly closed before the program terminates. Then, WHIZARD can be restarted with read_events set (which implies read_grids). The events already generated will be reread from the raw event file whizard.evx until the end is reached, and further events will be appended to this file. Event output files in other formats are re-generated on the fly. Thisy method can also be used to increase an event sample afterwards, spending CPU time essentially on newly generated events only.

4.9.3  Hard limits

There are situations where one has to control the amount of memory, CPU time, and disk space which the program uses. In particular, batch installations put limits on these resources and kill processes which exceed them.

Some of these issues can be handled by WHIZARD:

time_limit:

This parameter (in minutes) sets an internal limit on the wallclock time the run may use. When it is exceeded, the program terminates. Depending whether the termination occured during integration or simulation, the program can be restarted such that it will resume after the last grid or event written, respectively.

events_per_file:

This parameter takes care of event files becoming too large: once it is exceeded, a new file is begun. An external analysis would then process the files consecutively. Note that this does not apply to the internal event file whizard.evx which cannot be split; if this file becomes too large, one has to switch it off by write_events_raw=F.

bytes_per_file:

The analogous parameter which counts the number of bytes allowed for an event file. This relies on certain assumptions on the number of bytes an event takes in a particular format.

1

A better method is to give the user-created file a different name and use a symbolic link to whizard.prc instead. The same applies to whizard.in.

2

The MadGraph version implemented in WHIZARD does not treat multi-scalar vertices correctly. Don’t use it if you are interested in double (or triple) Higgs production. CompHEP uses Feynman gauge (i.e., includes ghosts as external states) for a more efficient evaluation of electroweak processes, while the other generators use unitarity gauge. This should not be visible to the user. Note also that CompHEP provides polarization-averaged matrix elements only. There is a CompHEP model using unitarity gauge, SM_ug, which turns out to be less efficient in many cases).

3

The CompHEP version currently has only the anomalous quartic couplings. The MadGraph version is just the Standard Model.

4

The MSSM is supported only by the O’Mega matrix element generator, while the included MadGraph and CompHEP versions will fall back to the Standard Model.

5

Furthermore, by writing the Feynman rules into a Lagrangian file, one may easily modify or set up a new model for CompHEP, which then can be used by WHIZARD. Modifying the model for O’Mega requires changing program sources, however. The same holds for MadGraph, if you are interested in Lorentz or Dirac structures that are not present in the SM.

6

It does not work for CompHEP matrix elements, however.

7

The current implementation is not well suited for single-particle processes if the initial-state energy distributions are strongly peaked. A similar caveat holds for processes with an s-channel resonance, for instance, e⁺e⁻→µ⁺µ⁻ with a continuous beamstrahlung/ISR spectrum for the incoming electrons, integrated over the whole energy range. Nevertheless, with large statistics reasonable results can be obtained. Some improvement has been achieved by using a specific mapping for s-channel like topologies.

8

Not available for CompHEP or MadGraph matrix elements)

9

Naively summing over identical particles in the final state may result in double-counting, and becomes very inefficient if many combinations are allowed. The ‘c’ option takes care of both issues (see below).

10

If the list of flavor sums or the number of particles in the process becomes very long, it may happen that the compiler complains for the number of continuation lines becoming too large. This problem is easily fixed by increasing the number of continuation lines accepted by the compiler. You should consult your Fortran compiler documentation for the appropriate flag which may then be set in the config.site file. Don’t forget to clean up and configure before restarting the compilation.

11

And also for CompHEP.

12

CompHEP, unfortunately, reorders particles according its own conventions.

13

Those that are part of the WHIZARD system. Precompiled libraries such as the LHAPDF have to be present on the target system and can be enabled there if they are needed.

14

For MadGraph in particular, if there is an ambiguity you may have to set the QCD order as an option.

15

Note that if WHIZARD does not find the corresponding file it takes its default file.

16

Note that we do not keep this updated. The parameter file of version 2, however, contain the most recent PDG values

17

This works only if PYTHIA has been linked. Also it does not work (properly) if WHIZARD and the CERN library containing PYTHIA have been compiled by different compilers, because PYTHIA or JETSET cannot connect the proper output unit in that case. So then a file like fort.<unit> will appear containing the PYTHIA/JETSET events.

18

For CompHEP matrix elements, and for O’Mega without the ’c’ option.