WHIZARD is hosted by Hepforge, IPPP Durham
Previous Up Next

Chapter 5  Detailed WHIZARD Steering: SINDARIN

5.1  Data and expressions

5.1.1  Real-valued objects

Real literals have their usual form, mantissa and, optionally, exponent:

0.  3.14  -.5  2.345e-3  .890E-023

Internally, real values are treated as double precision. The values are read by the Fortran library, so details depend on its implementation.

A special feature of SINDARIN is that numerics (real and integer) can be immediately followed by a physical unit. The supported units are presently hard-coded, they are

meV  eV  keV  MeV  GeV  TeV
nbarn  pbarn  fbarn  abarn
rad  mrad  degree
%

If a number is followed by a unit, it is automatically normalized to the corresponding default unit: 14.TeV is transformed into the real number 14000. Default units are GeV, fbarn, and rad. The % sign after a number has the effect that the number is multiplied by 0.01. Note that no checks for consistency of units are done, so you can add 1 meV + 3 abarn if you absolutely wish to. Omitting units is always allowed, in that case, the default unit is assumed.

Units are not treated as variables. In particular, you can’t write theta / degree, the correct form is theta / 1 degree.

There is a single predefined real constant, namely π which is referred to by the keyword pi. In addition, there is a single predefined complex constant, which is the complex unit i, being referred to by the keyword I.

The arithmetic operators are

+ - * / ^

with their obvious meaning and the usual precedence rules.

SINDARIN supports a bunch of standard numerical functions, mostly equivalent to their Fortran counterparts:

abs  conjg  sgn  mod  modulo
sqrt  exp  log  log10
sin  cos  tan  asin  acos  atan
sinh  cosh  tanh

(Unlike Fortran, the sgn function takes only one argument and returns 1., or −1.) The function argument is enclosed in brackets: sqrt (2.), tan (11.5 degree).

There are two functions with two real arguments:

max  min

Example: real lighter_mass = min (mZ, mH)

The following functions of a real convert to integer:

int  nint  floor  ceiling

and this converts to complex type:

complex

Real values can be compared by the following operators, the result is a logical value:

==  <>
>  <  >=  <=

In SINDARIN, it is possible to have more than two operands in a logical expressions. The comparisons are done from left to right. Hence,

115 GeV < mH < 180 GeV

is valid SINDARIN code and evaluates to true if the Higgs mass is in the given range.

Tests for equality and inequality with machine-precision real numbers are notoriously unreliable and should be avoided altogether. To deal with this problem, SINDARIN has the possibility to make the comparison operators “fuzzy” which should be read as “equal (unequal) up to an absolute tolerance”, where the tolerance is given by the real-valued intrinsic variable tolerance. This variable is initially zero, but can be set to any value (for instance, tolerance = 1.e-13 by the user. Note that for non-zero tolerance, operators like == and <> or < and > are not mutually exclusive1.

5.1.2  Integer-valued objects

Integer literals are obvious:

1  -98765  0123

Integers are always signed. Their range is the default-integer range as determined by the Fortran compiler.

Like real values, integer values can be followed by a physical unit: 1 TeV, 30 degree. This actually transforms the integer into a real.

Standard arithmetics is supported:

+ - * / ^

It is important to note that there is no fraction datatype, and pure integer arithmetics does not convert to real. Hence 3/4 evaluates to 0, but 3 GeV / 4 GeV evaluates to 0.75.

Since all arithmetics is handled by the underlying Fortran library, integer overflow is not detected. If in doubt, do real arithmetics.

Integer functions are more restricted than real functions. We support the following:

abs  sgn  mod  modulo
max  min

and the conversion functions

real  complex

Comparisons of integers among themselves and with reals are possible using the same set of comparison operators as for real values. This includes the operators with a finite tolerance.

5.1.3  Complex-valued objects

Complex variables and values are currently not yet used by the physics models implemented in WHIZARD. There complex input coupling constants are always split into their real and imaginary parts (or modulus and phase). They are exclusively available for arithmetic calculations.

There is no form for complex literals. Complex values must be created via an arithmetic expression,

complex c = 1 + 2 * I

where the imaginary unit I is predefined as a constant.

The standard arithmetic operations are supported (also mixed with real and integer). Support for functions is currently still incomplete, among the supported functions there are sqrt, log, exp.

5.1.4  Logical-valued objects

There are two predefined logical constants, true and false. Logicals are not equivalent to integers (like in C) or to strings (like in PERL), but they make up a type of their own. Only in printf output, they are treated as strings, that is, they require the %s conversion specifier.

The names of logical variables begin with a question mark ?. Here is the declaration of a logical user variable:

logical ?higgs_decays_into_tt = mH > 2 * mtop  

Logical expressions use the standard boolean operations

or  and  not

The results of comparisons (see above) are logicals.

There is also a special logical operator with lower priority, concatenation by a semicolon:

lexpr1 ; lexpr2

This evaluates lexpr1 and throws its result away, then evaluates lexpr2 and returns that result. This feature is to used with logical expressions that have a side effect, namely the record function within analysis expressions.

The primary use for intrinsic logicals are flags that change the behavior of commands. For instance, ?unweighted = true and ?unweighted = false switch the unweighting of simulated event samples on and off.

5.1.5  String-valued objects and string operations

String literals are enclosed in double quotes: "This is a string." The empty string is "". String variables begin with the dollar sign: $. There is only one string operation, concatenation

string $foo = "abc" & "def"

However, it is possible to transform variables and values to a string using the sprintf function. This function is an interface to the system’s C function sprintf with some restrictions and modifications. The allowed conversion specifiers are

%d  %i (integer)
%e  %f  %g  %E  %F  %G (real)
%s (string and logical)

The conversions can use flag parameter, field width, and precision, but length modifiers are not supported since they have no meaning for the application. (See also Sec. 5.10.)

The sprintf function has the syntax

sprintf format-string (arg-list)

This is an expression that evaluates to a string. The format string contains the mentioned conversion specifiers. The argument list is optional. The arguments are separated by commas. Allowed arguments are integer, real, logical, and string variables, and numeric expressions. Logical and string expressions can also be printed, but they have to be dressed as anonymous variables. A logical anonymous variable has the form ?(logical_expr) (example: ?(mH > 115 GeV)). A string anonymous variable has the form $(string-expr).

Example:

string $unit = "GeV"
string $str = sprintf "mW = %f %s" (mW, $unit)

The related printf command with the same syntax prints the formatted string to standard output2.

5.2  Particles and (sub)events

5.2.1  Particle aliases

A particle species is denoted by its name as a string: "W+". Alternatively, it can be addressed by an alias. For instance, the W+ boson has the alias Wp. Aliases are used like variables in a context where a particle species is expected, and the user can specify his/her own aliases.

An alias may either denote a single particle species or a class of particles species. A colon : concatenates particle names and aliases to yield multi-species aliases:

alias quark = u:d:s
alias wboson = "W+":"W-"

Such aliases are used for defining processes with summation over flavors, and for defining classes of particles for analysis.

Each model files define both names and (single-particle) aliases for all particles it contains. Furthermore, it defines the class aliases colored and charged which are particularly useful for event analysis.

5.2.2  Subevents

Subevents are sets of particles, extracted from an event. The sets are unordered by default, but may be ordered by appropriate functions. Obviously, subevents are meaningful only in a context where an event is available. The possible context may be the specification of a cut, weight, scale, or analysis expression.

To construct a simple subevent, we put a particle alias or an expression of type particle alias into square brackets:

["W+"]  [u:d:s]  [colored]

These subevents evaluate to the set of all W+ bosons (to be precise, their four-momenta), all u, d, or s quarks, and all colored particles, respectively.

A subevent can contain pseudoparticles, i.e., particle combinations. That is, the four-momenta of distinct particles are combined (added conmponent-wise), and the results become subevent elements just like ordinary particles.

The (pseudo)particles in a subevent are non-overlapping. That is, for any of the particles in the original event, there is at most one (pseudo)particle in the subevent in which it is contained.

Sometimes, variables (actually, named constants) of type subevent are useful. Subevent variables are declared by the subevt keyword, and their names carry the prefix @. Subevent variables exist only within the scope of a cuts (or scale, analysis, etc.) macro, which is evaluated in the presence of an actual event. In the macro body, they are assigned via the let construct:

cuts =
  let subevt @jets = select if Pt > 10 GeV [colored]
  in
  all Theta > 10 degree [@jets, @jets]

In this expression, we first define @jets to stand for the set of all colored partons with pT>10 GeV. This abbreviation is then used in a logical expression, which evaluates to true if all relative angles between distinct jets are greater than 10 degree.

We note that the example also introduces pairs of subevents: the square bracket with two entries evaluates to the list of all possible pairs which do not overlap. The objects within square brackets can be either subevents or alias expressions. The latter are transformed into subevents before they are used.

As a special case, the original event is always available as the predefined subevent @evt.

5.2.3  Subevent functions

There are several functions that take a subevent (or an alias) as an argument and return a new subevent. Here we describe them:

collect

collect [particles]
collect if condition [particles]
collect if condition [particles, ref_particles]

First version: collect all particle momenta in the argument and combine them to a single four-momentum. The particles argument may either be a subevt expression or an alias expression. The result is a one-entry subevt. In the second form, only those particles are collected which satisfy the condition, a logical expression. Example: collect if Pt > 10 GeV [colored]

The third version is useful if you want to put binary observables (i.e., observables constructed from two different particles) in the condition. The ref_particles provide the second argument for binary observables in the condition. A particle is taken into account if the condition is true with respect to all reference particles that do not overlap with this particle. Example: collect if Theta > 5 degree [photon, charged]: combine all photons that are separated by 5 degrees from all charged particles.

cluster

NOTE: This is an experimental feature, available from version 2.2.1 on.

cluster [particles]
cluster if condition [particles]

First version: collect all particle momenta in the argument and cluster them to a set of jets. The particles argument may either be a subevt expression or an alias expression. The result is a one-entry subevt. In the second form, only those particles are clustered which satisfy the condition, a logical expression. Example: cluster if Pt > 10 GeV [colored]

This command is available from WHIZARD version 2.2.1 on, and only if the FastJet package has been installed and linked with WHIZARD (cf. Sec.2.2.9); in a future version of WHIZARD it is foreseen to have also an intrinsic clustering package inside WHIZARD which will be able to support some of the clustering algorithms below. To use it in an analysis, you have to set the variable jet_algorithm to one of the predefined jet-algorithm values (integer constants):

kt_algorithm
cambridge_algorithm
antikt_algorithm
genkt_algorithm
cambridge_for_passive_algorithm
genkt_for_passive_algorithm
ee_kt_algorithm
ee_genkt_algorithm
plugin_algorithm

and the variable jet_r to the desired R parameter value, as appropriate for the analysis and the jet algorithm. Example:

  jet_algorithm = antikt_algorithm
  jet_r = 0.7
  cuts = all Pt > 15 GeV [cluster if Pt > 5 GeV [colored]]

combine

combine [particles_1, particles_2]
combine if condition [particles_1, particles_2]

Make a new subevent of composite particles. The composites are generated by combining all particles from subevent particles_1 with all particles from subevent particles_2 in all possible combinations. Overlapping combinations are excluded, however: if a (composite) particle in the first argument has a constituent in common with a composite particle in the second argument, the combination is dropped. In particular, this applies if the particles are identical.

If a condition is provided, the combination is done only when the logical expression, applied to the particle pair in question, returns true. For instance, here we reconstruct intermediate W bosons:

let @W_candidates = combine if 70 GeV < M < 80 GeV ["mu-", "numubar"] 
in ...

Note that the combination may fail, so the resulting subevent could be empty.

operator +

If there is no condition, the + operator provides a convenient shorthand for the combine command. In particular, it can be used if there are several particles to combine. Example:

cuts = any 170 GeV < M < 180 GeV [b + lepton + invisible]

select

select if condition [particles]
select if condition [particles, ref_particles]

One argument: select all particles in the argument that satisfy the condition and drop the rest. Two arguments: the ref_particles provide a second argument for binary observables. Select particles if the condition is satisfied for all reference particles.

extract

extract [particles]
extract index index-value [particles]

Return a single-particle subevent. In the first version, it contains the first particle in the subevent particles. In the second version, the particle with index index-value is returned, where index-value is an integer expression. If its value is negative, the index is counted from the end of the subevent.

The order of particles in an event or subevent is not always well-defined, so you may wish to sort the subevent before applying the extract function to it.

sort

sort [particles]
sort by observable [particles]
sort by observable [particles, ref_particle]

Sort the subevent according to some criterion. If no criterion is supplied (first version), the subevent is sorted by increasing PDG code (first particles, then antiparticles). In the second version, the observable is a real expression which is evaluated for each particle of the subevent in turn. The subevent is sorted by increasing value of this expression, for instance:

let @sorted_evt = sort by Pt [@evt]
in ...

In the third version, a reference particle is provided as second argument, so the sorting can be done for binary observables. It doesn’t make much sense to have several reference particles at once, so the sort function uses only the first entry in the subevent ref-particle, if it has more than one.

join

join [particles, new_particles]
join if condition [particles, new_particles]

This commands appends the particles in subevent new_particles to the subevent particles, i.e., it joins the two particle sets. To be precise, a (pseudo)particle from new_particles is only appended if it does not overlap with any of the (pseudo)particles present in particles, so the function will not produce overlapping entries.

In the second version, each particle from new_particles is also checked with all particles in the first set whether condition is fulfilled. If yes, and there is no overlap, it is appended, otherwise it is dropped.

operator &

Subevents can also be concatenated by the operator &. This effectively applies join to all operands in turn. Example:

let @visible = 
    select if Pt > 10 GeV and E > 5 GeV [photon]
  & select if Pt > 20 GeV and E > 10 GeV [colored]
  & select if Pt > 10 GeV [lepton]
in ...

5.2.4  Calculating observables

Observables (invariant mass M, energy E, …) are used in expressions just like ordinary numeric variables. By convention, their names start with a capital letter. They are computed using a particle momentum (or two particle momenta) which are taken from a subsequent subevent argument.

We can extract the value of an observable for an event and make it available for computing the scale value, or for histogramming etc.:

eval

eval expr [particles]
eval expr [particles_1, particles_2]

The function eval takes an expression involving observables and evaluates it for the first momentum (or momentum pair) of the subevent (or subevent pair) in square brackets that follows the expression. For example,

  eval Pt [colored]

evaluates to the transverse momentum of the first colored particle,

  eval M [@jets, @jets]

evaluates to the invariant mass of the first distinct pair of jets (assuming that @jets has been defined in a let construct), and

  eval E - M [combine [e1, N1]]

evaluates to the difference of energy and mass of the combination of the first electron-neutrino pair in the event.

The last example illustrates why observables are treated like variables, even though they are functions of particles: the eval construct with the particle reference in square brackets after the expression allows to compute derived observables – observables which are functions of new observables – without the need for hard-coding them as new functions.

5.2.5  Cuts and event selection

Instead of a numeric value, we can use observables to compute a logical value.

all

all logical_expr [particles]
all logical_expr [particles_1, particles_2]

The all construct expects a logical expression and one or two subevent arguments in square brackets.

  all Pt > 10 GeV [charged]
  all 80 GeV < M < 100 GeV [lepton, antilepton]

In the second example, lepton and antilepton should be aliases defined in a let construct. (Recall that aliases are promoted to subevents if they occur within square brackets.)

This construction defines a cut. The result value is true if the logical expression evaluates to true for all particles in the subevent in square brackets. In the two-argument case it must be true for all non-overlapping combinations of particles in the two subevents. If one of the arguments is the empty subevent, the result is also true.

any

any logical_expr [particles]
any logical_expr [particles_1, particles_2]

The any construct is true if the logical expression is true for at least one particle or non-overlapping particle combination:

  any E > 100 GeV [photon]

This defines a trigger or selection condition. If a subevent argument is empty, it evaluates to false

no

no logical_expr [particles]
no logical_expr [particles_1, particles_2]

The no construct is true if the logical expression is true for no single one particle or non-overlapping particle combination:

  no 5 degree < Theta < 175 degree ["e-":"e+"]

This defines a veto condition. If a subevent argument is empty, it evaluates to true. It is equivalent to not any…, but included for notational convenience.

5.2.6  More particle functions

count

count [particles]
count [particles_1, particles_2]
count if logical-expr [particles]
count if logical-expr [particles, ref_particles]

This counts the number of events in a subevent, the result is of type int. If there is a conditional expression, it counts the number of particle in the subevent that pass the test. If there are two arguments, it counts the number of non-overlapping particle pairs (that pass the test, if any).

Predefined observables

The following real-valued observables are available in SINDARIN for use in eval, all, any, no, and count constructs. The argument is always the subevent or alias enclosed in square brackets.

  • M2
    • One argument: Invariant mass squared of the (composite) particle in the argument.
    • Two arguments: Invariant mass squared of the sum of the two momenta.
  • M
    • Signed square root of M2: positive if M2>0, negative if M2<0.
  • E
    • One argument: Energy of the (composite) particle in the argument.
    • Two arguments: Sum of the energies of the two momenta.
  • Px, Py, Pz
    • Like E, but returning the spatial momentum components.
  • P
    • Like E, returning the absolute value of the spatial momentum.
  • Pt, Pl
    • Like E, returning the transversal and longitudinal momentum, respectively.
  • Theta
    • One argument: Absolute polar angle in the lab frame
    • Two arguments: Angular distance of two particles in the lab frame.
  • Theta_star Only with two arguments, gives the relative polar angle of the two momenta in the rest system of the momentum sum (i.e. mother particle).
  • Phi
    • One argument: Absolute azimuthal angle in the lab frame
    • Two arguments: Azimuthal distance of two particles in the lab frame
  • Rap, Eta
    • One argument: rapidity / pseudorapidity
    • Two arguments: rapidity / pseudorapidity difference
  • Dist
    • Two arguments: Distance on the η-φ cylinder, i.e., √Δη2 + Δφ2
  • kT
    • Two arguments: kT jet clustering variable: 2 min(Ej12, Ej22) / Q2 × (1 − cosθj1,j2). At the moment, Q2 = 1 GeV2.

There is also an integer-valued observable:

  • PDG
    • One argument: PDG code of the particle. For a composite particle, the code is undefined (value 0).

5.3  Physics Models

A physics model is a combination of particles, numerical parameters (masses, couplings, widths), and Feynman rules. Many physics analyses are done in the context of the Standard Model (SM). The SM is also the default model for WHIZARD. Alternatively, you can choose a subset of the SM (QED or QCD), variants of the SM (e.g., with or without nontrivial CKM matrix), or various extensions of the SM. The complete list is displayed in Table 10.1.

The model definitions are contained in text files with filename extension .mdl, e.g., SM.mdl, which are located in the share/models subdirectory of the WHIZARD installation. These files are easily readable, so if you need details of a model implementation, inspect their contents. The model file contains the complete particle and parameter definitions as well as their default values. It also contains a list of vertices. This is used only for phase-space setup; the vertices used for generating amplitudes and the corresponding Feynman rules are stored in different files within the O’Mega source tree.

In a SINDARIN script, a model is a special object of type model. There is always a current model. Initially, this is the SM, so on startup WHIZARD reads the SM.mdl model file and assigns its content to the current model object. (You can change the default model by the --model option on the command line. Also the preloading of a model can be switched off with the --no-model option) Once the model has been loaded, you can define processes for the model, and you have all independent model parameters at your disposal. As noted before, these are intrinsic parameters which need not be declared when you assign them a value, for instance:

  mW = 80.33 GeV
  wH = 243.1 MeV  

Other parameters are derived. They can be used in expressions like any other parameter, they are also intrinsic, but they cannot be modified directly at all. For instance, the electromagnetic coupling ee is a derived parameter. If you change either GF (the Fermi constant), mW (the W mass), or mZ (the Z mass), this parameter will reflect the change, but setting it directly is an error. In other words, the SM is defined within WHIZARD in the GF-mW-mZ scheme. (While this scheme is unusual for loop calculations, it is natural for a tree-level event generator where the Z and W poles have to be at their experimentally determined location3.)

The model also defines the particle names and aliases that you can use for defining processes, cuts, or analyses.

If you would like to generate a SUSY process instead, for instance, you can assign a different model (cf. Table 10.1) to the current model object:

  model = MSSM  

This assignment has the consequence that the list of SM parameters and particles is replaced by the corresponding MSSM list (which is much longer). The MSSM contains essentially all SM parameters by the same name, but in fact they are different parameters. This is revealed when you say

  model = SM  
  mb = 5.0 GeV  
  model = MSSM
  show (mb)  

After the model is reassigned, you will see the MSSM value of mb which still has its default value, not the one you have given. However, if you revert to the SM later,

  model = SM  
  show (mb)  

you will see that your modification of the SM’s mb value has been remembered. If you want both mass values to agree, you have to set them separately in the context of their respective model. Although this might seem cumbersome at first, it is nevertheless a sensible procedure since the parameters defined by the user might anyhow not be defined or available for all chosen models.

When using two different models which need an SLHA input file, these have to be provided for both models.

Within a given scope, there is only one current model. The current model can be reset permanently as above. It can also be temporarily be reset in a local scope, i.e., the option body of a command or the body of a scan loop. It is thus possible to use several models within the same script. For instance, you may define a SUSY signal process and a pure-SM background process. Each process depends only on the respective model’s parameter set, and a change to a parameter in one of the models affects only the corresponding process.

5.4  Processes

The purpose of WHIZARD is the integration and simulation of high-energy physics processes: scatterings and decays. Hence, process objects play the central role in SINDARIN scripts.

A SINDARIN script may contain an arbitrary number of process definitions. The initial states need not agree, and the processes may belong to different physics models.

5.4.1  Process definition

A process object is defined in a straightforward notation. The definition syntax is straightforward:

process process-id = incoming-particles => outgoing-particles

Here are typical examples:

  process w_pair_production = e1, E1 => "W+", "W-"
  process zdecay = Z => u, ubar

Throughout the program, the process will be identified by its process-id, so this is the name of the process object. This identifier is arbitrary, chosen by the user. It follows the rules for variable names, so it consists of alphanumeric characters and underscores, where the first character is not numeric. As a special rule, it must not contain upper-case characters. The reason is that this name is used for identifying the process not just within the script, but also within the Fortran code that the matrix-element generator produces for this process.

After the equals sign, there follow the lists of incoming and outgoing particles. The number of incoming particles is either one or two: scattering processes and decay processes. The number of outgoing particles should be two or larger (as 2→ 1 processes are proportional to a δ function they can only be sensibly integrated when using a structure function like a hadron collider PDF or a beamstrahlung spectrum.). There is no hard upper limit; the complexity of processes that WHIZARD can handle depends only on the practical computing limitations (CPU time and memory). Roughly speaking, one can assume that processes up to 2→ 6 particles are safe, 2→ 8 processes are feasible given sufficient time for reaching a stable integration, while more complicated processes are largely unexplored.

We emphasize that in the default setup, the matrix element of a physics process is computed exactly in leading-order perturbation theory, i.e., at tree level. There is no restriction of intermediate states, the result always contains the complete set of Feynman graphs that connect the initial with the final state. If the result would actually be expanded in Feynman graphs (which is not done by the O’Mega matrix element generator that WHIZARD uses), the number of graphs can easily reach several thousands, depending on the complexity of the process and on the physics model.

More details about the different methods for quantum field-theoretical matrix elements can be found in Chap. 9. In the following, we will discuss particle names, options for processes like restrictions on intermediate states, parallelization, flavor sums and process components for inclusive event samples (process containers).

5.4.2  Particle names

The particle names are taken from the particle definition in the current model file. Looking at the SM, for instance, the electron entry in share/models/SM.mdl reads

particle E_LEPTON 11
  spin 1/2  charge  -1   isospin -1/2  
  name "e-" e1 electron e
  anti "e+" E1 positron 
  tex_name "e^-"
  tex_anti "e^+"
  mass me

This tells that you can identify an electron either as "e-", e1, electron, or simply e. The first version is used for output, but needs to be quoted, because otherwise SINDARIN would interpret the minus sign as an operator. (Technically, unquoted particle identifiers are aliases, while the quoted versions – you can say either e1 or "e1" – are names. On input, this makes no difference.) The alternative version e1 follows a convention, inherited from CompHEP [57], that particles are indicated by lower case, antiparticles by upper case, and for leptons, the generation index is appended: e2 is the muon, e3 the tau. These alternative names need not be quoted because they contain no special characters.

In Table 5.1, we list the recommended names as well as mass and width parameters for all SM particles. For other models, you may look up the names in the corresponding model file.


 ParticleOutput nameAlternative namesMassWidth
Leptonsee-e1  electronme 
 e+e+E1  positronme 
 µmu-e2  muonmmu 
 µ+mu+E2mmu 
 τtau-e3  tauonmtau 
 τ+tau+E3mtau 
Neutrinosνenuen1  
 νenuebarN1  
 νµnumun2  
 νµnumubarN2  
 ντnutaun3  
 ντnutaubarN3  
Quarksdddown  
 ddbarD  
 uuup  
 ūubarU  
 ssstrangems 
 ssbarSms 
 cccharmmc 
 ccbarCmc 
 bbbottommb 
 bbbarBmb 
 tttopmtopwtop
 ttbarTmtopwtop
Vector bosonsgglg  G  gluon  
 γAgamma  photon  
 ZZ mZwZ
 W+W+WpmWwW
 WW-WmmWwW
Scalar bosonsHHh  HiggsmHwH
Table 5.1: Names that can be used for SM particles. Also shown are the intrinsic variables that can be used to set mass and width, if applicable.

Where no mass or width parameters are listed in the table, the particle is assumed to be massless or stable, respectively. This is obvious for particles such as the photon. For neutrinos, the mass is meaningless to particle physics collider experiments, so it is zero. For quarks, the u or d quark mass is unobservable directly, so we also set it zero. For the heavier quarks, the mass may play a role, so it is kept. (The s quark is borderline; one may argue that its mass is also unobservable directly.) On the other hand, the electron mass is relevant, e.g., in photon radiation without cuts, so it is not zero by default.

It pays off to set particle masses to zero, if the approximation is justified, since fewer helicity states will contribute to the matrix element. Switching off one of the helicity states of an external fermion speeds up the calculation by a factor of two. Therefore, script files will usually contain the assignments

  me = 0  mmu = 0  ms = 0  mc = 0

unless they deal with processes where this simplification is phenomenologically unacceptable. Often mτ and mb can also be neglected, but this excludes processes where the Higgs couplings of τ or b are relevant.

Setting fermion masses to zero enables, furthermore, the possibility to define multi-flavor aliases

  alias q = d:u:s:c
  alias Q = D:U:S:C  

and handle processes such as

  process two_jets_at_ilc = e1, E1 => q, Q
  process w_pairs_at_lhc = q, Q => Wp, Wm  

where a sum over all allowed flavor combination is automatically included. For technical reasons, such flavor sums are possible only for massless particles (or more general for mass-degenerate particles). If you want to generate inclusive processes with sums over particles of different masses (e.g. summing over W/Z in the final state etc.), confer below the section about process components, Sec. 5.4.4.

Assignments of masses, widths and other parameters are actually in effect when a process is integrated, not when it is defined. So, these assignments may come before or after the process definition, with no significant difference. However, since flavor summation requires masses to be zero, the assignments may be put before the alias definition which is used in the process.

The muon, tau, and the heavier quarks are actually unstable. However, the width is set to zero because their decay is a macroscopic effect and, except for the muon, affected by hadron physics, so it is not described by WHIZARD. (In the current WHIZARD setup, all decays occur at the production vertex. A future version may describe hadronic physics and/or macroscopic particle propagation, and this restriction may be eventually removed.)

5.4.3  Options for processes

The process definition may contain an optional argument:

process process-id = incoming-particles => outgoing-particles {options…}

The options are a SINDARIN script that is executed in a context local to the process command. The assignments it contains apply only to the process that is defined. In the following, we describe the set of potentially useful options (which all can be also set globally):

Model reassignment

It is possible to locally reassign the model via a model = statment, permitting the definition of process using a model other than the globally selected model. The process will retain this association during integration and event generation.

Restriction on intermediate states

Another useful option is the setting

$restrictions = string

This option allows to select particular classes of Feynman graphs for the process when using the O’Mega matrix element generator. The $restrictions string specifies propagators that the graph must contain. Here is an example:

  process zh_invis = e1, E1 => n1:n2:n3, N1:N2:N3, H { $restrictions = "1+2 ~ Z" }

The complete process ee+ → ννH, summed over all neutrino generations, contains both ZH pair production (Higgs-strahlung) and W+WH fusion. The restrictions string selects the Higgs-strahlung graph where the initial electrons combine to a Z boson. Here, the particles in the process are consecutively numbered, starting with the initial particles. An alternative for the same selection would be $restrictions = "3+4 ~ Z". Restrictions can be combined using &&, for instance

  $restrictions = "1+2 ~ Z && 3 + 4 ~ Z"

which is redundant here, however.

The restriction keeps the full energy dependence in the intermediate propagator, so the Breit-Wigner shape can be observed in distributions. This breaks gauge invariance, in particular if the intermediate state is off shell, so you should use the feature only if you know the implications. For more details, cf. the Chap. 9 and the O’Mega manual.

Other options

There are some further options that the O’Mega matrix-element generator can take. If desired, any string of options that is contained in this variable

$omega_flags = string

will be copied verbatim to the O’Mega call, after all other options.

One important application is the scheme of treating the width of unstable particles in the t-channel. This is modified by the model: class of O’Mega options.

It is well known that for some processes, e.g., single W production from photon-W fusion, gauge invariance puts constraints on the treatment of the unstable-particle width. By default, O’Mega puts a nonzero width in the s channel only. This correctly represents the resummed Dyson series for the propagator, but it violates QED gauge invariance, although the effect is only visible if the cuts permit the photon to be almost on-shell.

An alternative is

$omega_flags = "-model:fudged_width"

which puts zero width in the matrix element, so that gauge cancellations hold, and reinstates the s-channel width in the appropriate places by an overall factor that multiplies the whole matrix element.

Another possibility is

$omega_flags = "-model:constant_width"

which puts the width both in the s and in the t channel everywhere.

Note that both options apply only to charged unstable particles, such as the W boson.

Multithreaded calculation of helicity sums via OpenMP

On multicore and / or multiprocessor systems, it is possible to speed up the calculation by using multiple threads to perform the helicity sum in the matrix element calculation. As the processing time used by WHIZARD is not used up solely in the matrix element, the speedup thus achieved varies greatly depending on the process under consideration; while simple processes without flavor sums do not profit significantly from this parallelization, the computation time for processes involving flavor sums with four or more particles in the final state is typically reduced by a factor between two and three when utilizing four parallel threads.

The parallization is implemented using OpenMP and requires WHIZARD to be compiled with an OpenMP aware compiler and the appropiate compiler flags This is done in the configuration step, cf. Sec. 2.3.

As with all OpenMP programs, the default number of threads used at runtime is up to the compiler runtime support and typically set to the number of independent hardware threads (cores / processors / hyperthreads) available in the system. This default can be adjusted by setting the OMP_NUM_THREADS environment variable prior to calling WHIZARD. Alternatively, the available number of threads can be reset anytime by the SINDARIN parameter openmp_num_threads. Note however that the total number of threads that can be sensibly used is limited by the number of nonvanishing helicity combinations.

5.4.4  Process components

It was mentioned above that processes with flavor sums (in the initial or final state or both) have to be mass-degenerate (in most cases massless) in all particles that are summed over at a certain position. This condition is necessary in order to use the same phase-space parameterization and integration for the flavor-summed process. However, in many applications the user wants to handle inclusive process definitions, e.g. by defining inclusive decays, inclusive SUSY samples at hadron colliders (gluino pairs, squark pairs, gluino-squark associated production), or maybe lepton-inclusive samples where the tau and muon mass should be kept at different values. In WHIZARD  from version v2.2.0 on, there is the possibility to define such inclusive process containers. The infrastructure for this feature is realized via so-called process components: processes are allowed to contain several process components. Those components need not be provided by the same matrix element generator, e.g. internal matrix elements, O’Mega matrix elements, external matrix element (e.g. from a one-loop program, OLP) can be mixed. The very same infrastructure can also be used for next-to-leading order (NLO) calculations, containing the born with real emission, possible subtraction terms to make the several components infrared- and collinear finite, as well as the virtual corrections.

Here, we want to discuss the use for inclusive particle samples. There are several options, the simplest of which to add up different final states by just using the + operator in SINDARIN, e.g.:

   process multi_comp = e1, E1 => (e2, E2) + (e3, E3) + (A, A)

The brackets are not only used for a better grouping of the expressions, they are not mandatory for WHIZARD to interpret the sum correctly. When integrating, WHIZARD tells you that this a process with three different components:

| Initializing integration for process multi_comp_1_p1:
| ------------------------------------------------------------------------
| Process [scattering]: 'multi_comp'
|   Library name  = 'default_lib'
|   Process index = 1
|   Process components:
|     1: 'multi_comp_i1':   e-, e+ => m-, m+ [omega]
|     2: 'multi_comp_i2':   e-, e+ => t-, t+ [omega]
|     3: 'multi_comp_i3':   e-, e+ => A, A [omega]
| ------------------------------------------------------------------------  

A different phase-space setup is used for each different component. The integration for each different component is performed separately, and displayed on screen. At the end, a sum of all components is shown. All files that depend on the components are being attached an _i<n> where <n> is the number of the process component that appears in the list above: the Fortran code for the matrix element, the .phs file for the phase space parameterization, and the grid files for the VAMP Monte-Carlo integration (or any other integration method). However, there will be only one event file for the inclusive process, into which a mixture of events according to the size of the individual process component cross section enter.

More options are to specify additive lists of particles. WHIZARD then expands the final states according to tensor product algebra:

   process multi_tensor = e1, E1 => e2 + e3 + A, E2 + E3 + A

This gives the same three process components as above, but WHIZARD recognized that e.g. e e+ → µ γ is a vanishing process, hence the numbering is different:

| Process component 'multi_tensor_i2': matrix element vanishes
| Process component 'multi_tensor_i3': matrix element vanishes
| Process component 'multi_tensor_i4': matrix element vanishes
| Process component 'multi_tensor_i6': matrix element vanishes
| Process component 'multi_tensor_i7': matrix element vanishes
| Process component 'multi_tensor_i8': matrix element vanishes
| ------------------------------------------------------------------------
| Process [scattering]: 'multi_tensor'
|   Library name  = 'default_lib'
|   Process index = 1
|   Process components:
|     1: 'multi_tensor_i1':   e-, e+ => m-, m+ [omega]
|     5: 'multi_tensor_i5':   e-, e+ => t-, t+ [omega]
|     9: 'multi_tensor_i9':   e-, e+ => A, A [omega]
| ------------------------------------------------------------------------  

Identical copies of the same process that would be created by expanding the tensor product of final states are eliminated and appear only once in the final sum of process components.

Naturally, inclusive process definitions are also available for decays:

process multi_dec = Wp => E2 + E3, n2 + n3  

This yields:

| Process component 'multi_dec_i2': matrix element vanishes
| Process component 'multi_dec_i3': matrix element vanishes
| ------------------------------------------------------------------------
| Process [decay]: 'multi_dec'
|   Library name  = 'default_lib'
|   Process index = 2
|   Process components:
|     1: 'multi_dec_i1':   W+ => mu+, numu [omega]
|     4: 'multi_dec_i4':   W+ => tau+, nutau [omega]
| ------------------------------------------------------------------------  

5.4.5  Compilation

Once processes have been set up, to make them available for integration they have to be compiled. More precisely, the matrix-element generator O’Mega (and it works similarly if a different matrix element method is chosen) is called to generate matrix element code, the compiler is called to transform this Fortran code into object files, and the linker is called to collect this in a dynamically loadable library. Finally, this library is linked to the program. From version v2.2.0 of WHIZARD this is no longer done by system calls of the OS but steered via process library Makefiles. Hence, the user can execute and manipulate those Makefiles in order to manually intervene in the particular steps, if he/she wants to do so.

All this is done automatically when an integrate, unstable, or simulate command is encountered for the first time. You may also force compilation explicitly by the command

  compile  

which performs all steps as listed above, including loading the generated library.

The Fortran part of the compilation will be done using the Fortran compiler specified by the string variable $fc and the compiler flags specified as $fcflags. The default settings are those that have been used for compiling WHIZARD itself during installation. For library compatibility, you should stick to the compiler. The flags may be set differently. They are applied in the compilation and loading steps, and they are processed by libtool, so libtool-specific flags can also be given.

WHIZARD has some precautions against unnecessary repetitions. Hence, when a compile command is executed (explicitly, or implicitly by the first integration), the program checks first whether the library is already loaded, and whether source code already exists for the requested processes. If yes, this code is used and no calls to O’Mega (or another matrix element method) or to the compiler are issued. Otherwise, it will detect any modification to the process configuration and regenerate the matrix element or recompile accordingly. Thus, a SINDARIN script can be executed repeatedly without rebuilding everything from scratch, and you can safely add more processes to a script in a subsequent run without having to worry about the processes that have already been treated.

This default behavior can be changed. By setting

  ?rebuild_library = true

code will be re-generated and re-compiled even if WHIZARD would think that this is unncessary. The same effect is achieved by calling WHIZARD with a command-line switch,

  /home/user$ whizard --rebuild_library

There are further rebuild switches which are described below. If everything is to be rebuilt, you can set a master switch ?rebuild or the command line option --rebuild. The latter can be abbreviated as a short command-line option:

  /home/user$ whizard -r

Setting this switch is always a good idea when starting a new project, just in case some old files clutter the working directory. When re-running the same script, possibly modified, the -r switch should be omitted, so the existing files can be reused.

5.4.6  Process libraries

Processes are collected in libraries. A script may use more than one library, although for most applications a single library will probably be sufficient.

The default library is default_lib. If you do not specify anything else, the processes you compile will be collected by a driver file default_lib.f90 which is compiled together with the process code and combined as a libtool archive default_lib.la, which is dynamically linked to the running WHIZARD process.

Once in a while, you work on several projects at once, and you didn’t care about opening a new working directory for each. If the -r option is given, a new run will erase the existing library, which may contain processes needed for the other project. You could omit -r, so all processes will be collected in the same library (this does not hurt), but you may wish to cleanly separate the projects. In that case, you should open a separate library for each project.

Again, there are two possibilities. You may start the script with the specification

  library = "my_lhc_proc"  

to open a library my_lhc_proc in place of the default library. Repeating the command with different arguments, you may introduce several libraries in the script. The active library is always the one specified last. It is possible to issue this command locally, so a particular process goes into its own library.

Alternatively, you may call WHIZARD with the option

  /home/user$ whizard --library=my_lhc_proc  

If several libraries are open simultaneously, the compile command will compile all libraries that the script has referenced so far. If this is not intended, you may give the command an argument,

  compile ("my_lhc_proc", "my_other_proc")

to compile only a specific subset.

The command

  show (library)

will display the contents of the actually loaded library together with a status code which indicates the status of the library and the processes within.

5.4.7  Stand-alone WHIZARD with precompiled processes

Once you have set up a process library, it is straightforward to make a special stand-alone WHIZARD executable which will have this library preloaded on startup. This is a matter of convenience, and it is also useful if you need a statically linked executable for reasons of profiling, batch processing, etc.

For this task, there is a variant of the compile command:

  compile as "my_whizard" ()

which produces an executable my_whizard. You can omit the library argument if you simply want to include everything. (Note that this command will not load a library into the current process, it is intended for creating a separate program that will be started independently.)

As an example, the script

  process proc1 = e1, E1 => e1, E1
  process proc2 = e1, E1 => e2, E2
  process proc3 = e1, E1 => e3, E3
  compile as "whizard-leptons" ()

will make a new executable program whizard-leptons. This program behaves completely identical to vanilla WHIZARD, except for the fact that the processes proc1, proc2, and proc3 are available without configuring them or loading any library.

5.5  Beams

Before processes can be integrated and simulated, the program has to know about the collider properties. They can be specified by the beams statement.

In the command script, it is irrelevant whether a beams statement comes before or after process specification. The integrate or simulate commands will use the beams statement that was issued last.

5.5.1  Beam setup

If the beams have no special properties, and the colliding particles are the incoming particles in the process themselves, there is no need for a beams statement at all. You only must specify the center-of-momentum energy of the collider by setting the value of √s, for instance

  sqrts = 14 TeV

The beams statement comes into play if

  • the beams have nontrivial structure, e.g., parton structure in hadron collision or photon radiation in lepton collision, or
  • the beams have non-standard properties: polarization, asymmetry, crossing angle.

Note that some of the abovementioned beam properties had not yet been reimplemented in the WHIZARD2 release series. From version v2.2.0 on all options of the legacy series WHIZARD1 are available again. From version v2.1 to version v2.2 of WHIZARD there has also been a change in possible options to the beams statement: in the early versions of WHIZARD2 (v2.0/v2.1), local options could be specified within the beam settings, e.g. beams = p, p sqrts = 14 TeV => pdf_builtin. These possibility has been abandoned from version v2.2 on, and the beams command does not allow for any optional arguments any more.

Hence, beam parameters can – with the exception of the specification of structure functions – be specified only globally:

  sqrts = 14 TeV   
  beams = p, p => lhapdf

It does not make any difference whether the value of sqrts is set before or after the beams statement, the last value found before an integrate or simulate is the relevant one. This in particularly allows to specify the beam structure, and then after that perform a loop or scan over beam energies, beam parameters, or structure function settings.

The beams statement also applies to particle decay processes, where there is only a single beam. Here, it is usually redundant because no structure functions are possible, and the energy is fixed to the decaying particle’s mass. However, it is needed for computing polarized decay, e.g.

  beams = Z 
  beams_pol_density = @(0)

where for a boson at rest, the polarization axis is defined to be the z axis.

Beam polarization is described in detail below in Sec. 5.6.

Note also that future versions of WHIZARD might give support for single-beam events, where structure functions for single particles indeed do make sense.

In the following sections we list the available options for structure functions or spectra inside WHIZARD and explain their usage. More about the physics of the implemented structure functions can be found in Chap. 9.

5.5.2  Asymmetric beams and Crossing angles

WHIZARD not only allows symmetric beam collisions, but basically arbitrary collider setups. In the case there are two different beam energies, the command

beams_momentum = <beam_mom1>, <beam_mom2>

allows to specify the momentum (or as well energies for massless particles) for the beams. Note that for scattering processes both values for the beams must be present. So the following to setups for 14 TeV LHC proton-proton collisions are equivalent:

beams = p, p => pdf_builtin
sqrts = 14 TeV

and

beams = p, p => pdf_builtin
beams_momentum = 7 TeV, 7 TeV

Asymmetric setups can be set by using different values for the two beam momenta, e.g. in a HERA setup:

beams = e, p => none, pdf_builtin beams_momentum = 27.5 GeV, 920 GeV

or for the BELLE experiment at the KEKB accelerator:

beams = e1, E1 beams_momentum = 8 GeV, 3.5 GeV

WHIZARD lets you know about the beam structure and calculates for you that the center of mass energy corresponds to 10.58 GeV:

| Beam structure: e-, e+
|   momentum = 8.000000000000E+00, 3.500000000000E+00
| Beam data (collision):
|   e-  (mass = 5.1099700E-04 GeV)
|   e+  (mass = 5.1099700E-04 GeV)
|   sqrts = 1.058300530253E+01 GeV
| Beam structure: lab and c.m. frame differ  

It is also possible to specify beams for decaying particles, where beams_momentum then only has a single argument, e.g.:

process zee = Z => "e-", "e+"
beams = Z
beams_momentum = 500 GeV
simulate (zee) { n_events = 100 }

This would corresponds to a beam of Z bosons with a momentum of 500 GeV. Note, however, that WHIZARD will always do the integration of the particle width in the particle’s rest frame, while the moving beam is then only taken into account for the frame of reference for the simulation.

Further options then simply having different beam energies describe a non-vanishing between the two incoming beams. Such concepts are quite common e.g. for linear colliders to improve the beam properties in the collimation region at the beam interaction points. Such crossing angles can be specified in the beam setup, too, using the beams_theta command:

beams = e1, E1
beams_momentum = 500 GeV, 500 GeV
beams_theta = 0, 10 degree

It is important that when a crossing angle is being specified, and the collision system consequently never is the center-of-momentum system, the beam momenta have to explicitly set. Besides a planar crossing angle, one is even able to rotate an azimuthal distance:

beams = e1, E1
beams_momentum = 500 GeV, 500 GeV
beams_theta = 0, 10 degree
beams_phi = 0, 45 degree

5.5.3  LHAPDF

For incoming hadron beams, the beams statement specifies which structure functions are used. The simplest example is the study of parton-parton scattering processes at a hadron-hadron collider such as LHC or Tevatron. The LHAPDF structure function set is selected by a syntax similar to the process setup, namely the example already shown above:

  beams = p, p => lhapdf

Note that there are slight differences in using the LHAPDF release series 6 and the older Fortran LHAPDF release series 5, at least concerning the naming conventions for the PDF sets 4. The above beams statement selects a default LHAPDF structure-function set for both proton beams (which is the CT10 central set for LHAPDF 6, and cteq6ll.LHpdf central set for LHAPDF5). The structure function will apply for all quarks, antiquarks, and the gluon as far as supported by the particular LHAPDF set. Choosing a different set is done by adding the filename as a local option to the lhapdf keyword:

  beams = p, p => lhapdf 
  $lhapdf_file = "MSTW2008lo68cl"

for the actual LHAPDF 6 series, and

  beams = p, p => lhapdf 
  $lhapdf_file = "MSTW2008lo68cl.LHgrid"

for LHAPDF5.Similarly, a member within the set is selected by the numeric variable lhapdf_member (for both release series of LHAPDF).

In some cases, different structure functions have to be chosen for the two beams. For instance, we may look at ep collisions:

  beams = "e-", p => none, lhapdf

Here, there is a list of two independent structure functions (each with its own option set, if applicable) which applies to the two beams.

Another mixed case is pγ collisions, where the photon is to be resolved as a hadron. The simple assignment

  beams = p, gamma => lhapdf, lhapdf_photon

will be understood as follows: WHIZARD selects the appropriate default structure functions (here we are using LHAPDF 5 as an example as the support of photon and pion PDFs in LHAPDF 6 has been dropped), cteq6ll.LHpdf for the proton and GSG960.LHgrid for the photon. The photon case has an additional integer-valued parameter lhapdf_photon_scheme. (There are also pion structure functions available.) For modifying the default, you have to specify separate structure functions

  beams = p, gamma => lhapdf, lhapdf_photon
  $lhapdf_file = ...
  $lhapdf_photon_file = ...

Finally, the scattering of elementary photons on partons is described by

  beams = p, gamma => lhapdf, none

Note that for LHAPDF version 5.7.1 or higher and for PDF sets which support it, photons can be used as partons.

There is one more option for the LHAPDF PDFs, namely to specify the path where the LHAPDF PDF sets reside: this is done with the string variable $lhapdf_dir = "<path-to-lhapdf>". Usually, it is not necessary to set this because WHIZARD detects this path via the lhapdf-config script during configuration, but in the case paths have been moved, or special files/special locations are to be used, the user can specify this location explicitly.

5.5.4  Built-in PDFs

In addition to the possibility of linking against LHAPDF, WHIZARD comes with a couple of built-in PDFs which are selected via the pdf_builtin keyword

  beams = p, p => pdf_builtin

The default PDF set is CTEQ6L, but other choices are also available by setting the string variable $pdf_builtin_set to an appropiate value. E.g, modifying the above setup to

  beams = p, p => pdf_builtin 
  $pdf_builtin_set = "mrst2004qedp"

would select the proton PDF from the MRST2004QED set. A list of all currently available PDFs can be found in Table 5.2.


TagNameNotesReferences
cteq6lCTEQ6L [58]
cteq6l1CTEQ6L1 [58]
cteq6dCTEQ6D [58]
cteq6mCTEQ6M [58]
mrst2004qedpMRST2004QED (proton)includes photon [59]
mrst2004qednMRST2004QED (neutron)includes photon [59]
mstw2008loMSTW2008LO [60]
mstw2008nloMSTW2008NLO [60]
mstw2008nnloMSTW2008NNLO [60]
ct10CT10 [61]
CJ12_maxCJ12_max [62]
CJ12_midCJ12_mid [62]
CJ12_minCJ12_min [62]
CJ15LOCJ15LO [63]
CJ15NLOCJ15NLO [63]
mmht2014loMMHT2014LO [64]
mmht2014nloMMHT2014NLO [64]
mmht2014nnloMMHT2014NNLO [64]
CT14LLCT14LLO [65]
CT14LCT14LO [65]
CT14NCT1414NLO [65]
CT14NNCT14NNLO [65]
Table 5.2: All PDF sets available as builtin sets. The two MRST2004QED sets also contain a photon.

The two MRST2004QED sets also contain the photon as a parton, which can be used in the same way as for LHAPDF from v5.7.1 on. Note, however, that there is no builtin PDF that contains a photon structure function. There is a beams structure function specifier pdf_builtin_photon, but at the moment this throws an error. It just has been implemented for the case that in future versions of WHIZARD a photon structure function might be included.

Note that in general only the data sets for the central values of the different PDFs ship with WHIZARD. Using the error sets is possible, i.e. it is supported in the syntax of the code, but you have to download the corresponding data sets from the web pages of the PDF fitting collaborations.

5.5.5  HOPPET b parton matching

When the HOPPET tool [66] for hadron-collider PDF structure functions and their manipulations are correctly linked to WHIZARD, it can be used for advanced calculations and simulations of hadron collider physics. Its main usage inside WHIZARD is for matching schemes between 4-flavor and 5-flavor schemes in b-parton initiated processes at hadron colliders. Note that in versions 2.2.0 and 2.2.1 it only worked together with LHAPDF version 5, while with the LHAPDF version 6 interface from version 2.2.2 on it can be used also with the modern version of PDFs from LHAPDF. Furthermore, from version 2.2.2, the HOPPET b parton matching also works for the builtin PDFs.

It depends on the corresponding process and the energy scales involved whether it is a better description to use the gbb splitting from the DGLAP evolution inside the PDF and just take the b parton content of a PDF, e.g. in BSM Higgs production for large tanβ: ppH with a partonic subprocess bbH, or directly take the gluon PDFs and use ppbb H with a partonic subprocess ggb b H. Elaborate schemes for a proper matching between the two prescriptions have been developed and have been incorporated into the HOPPET interface.

Another prime example for using these matching schemes is single top production at hadron colliders. Let us consider the following setup:

process proc1 = b, u => t, d
process proc2 = u, b => t, d
process proc3 = g, u => t, d, B { $restrictions = "2+4 ~ W+" }
process proc4 = u, g => t, d, B { $restrictions = "1+4 ~ W+" }

beams = p,p => pdf_builtin
sqrts = 14 TeV
?hoppet_b_matching = true

$sample = "single_top_matched"
luminosity = 1 / 1 fbarn
simulate (proc1, proc2, proc3, proc4)

The first two processes are single top production from b PDFs, the last two processes contain an explicit gbb splitting (the restriction, cf. Sec. 5.4.3 has been placed in order to single out the single top production signal process). PDFs are then chosen from the default builtin PDF (which is CTEQ6L), and the HOPPET matching routines are switched on by the flag ?hoppet_b_matching.

5.5.6  Lepton Collider ISR structure functions

Initial state QED radiation off leptons is an important feature at all kinds of lepton colliders: the radiative return to the Z resonance by ISR radiation was in fact the largest higher-order effect for the SLC and LEP I colliders. The soft-collinear and soft photon radiation can indeed be resummed/exponentiated to all orders in perturbation theory [7], while higher orders in hard-collinear photons have to be explicitly calculated order by order [8, 9]. WHIZARD has an intrinsic implementation of the lepton ISR structure function that includes all orders of soft and soft-collinear photons as well as up to the third order in hard-collinear photons. It can be switched on by the following statement:

beams = e1, E1 => isr  

As the ISR structure function is a single-beam structure function, this expression is synonymous for

beams = e1, E1 => isr, isr

The ISR structure function can again be applied to only one of the two beams, e.g. in a HERA-like setup:

beams = e1, p => isr, pdf_builtin  

Their are several options for the lepton-collider ISR structure function that are summarized in the following:


ParameterDefaultMeaning
isr_alpha0/intrinsicvalue of αQED for ISR
isr_order3max. order of hard-collinear photon emission
isr_mass0/intrinsicmass of the radiating lepton
isr_q_max0/√supper cutoff for ISR
?isr_recoilfalseflag to switch on recoil/pT
?isr_keep_energyfalserecoil flag to conserve energy in splitting

The maximal order of the hard-collinear photon emission taken into account by WHIZARD is set by the integer variable isr_order; the default is the maximally available order of three. With the variable isr_alpha, the value of the QED coupling constant αQED used in the ISR structure function can be set. The default is taken from the active physics model. The mass of the radiating lepton (in most cases the electron) is set by isr_mass; again the default is taken from the active physics model. Furthermore, the upper integration border for the ISR structure function which acts roughly as an upper hardness cutoff for the emitted photons, can be set through isr_q_max; if not set, the collider energy (possibly after beamstrahlung, cf. Sec. 5.5.7) √s (or √s) is taken. Finally, with the flag ?isr_recoil, the pT recoil of the emitting lepton against the photon radiation can be switched on; per default it is off. Note that WHIZARD accounts for the exclusive effects of ISR radiation at the moment by a single (hard, resolved) photon in the event; a more realistic treatment of exclusive ISR photons in simulation is foreseen for a future version.

For more information on the underlying physics, see Chap. 9.

5.5.7  Lepton Collider Beamstrahlung

At linear lepton colliders, the macroscopic electromagnetic interaction of the bunches leads to a distortion of the spectrum of the bunches that is important for an exact simulation of the beam spectrum. There are several methods to account for these effects. The most important tool to simulate classical beam-beam interactions in lepton-collider physics is GuineaPig++ [10, 11, 12]. A direct interface between this tool GuineaPig++ and WHIZARD had existed as an inofficial add-on to the legacy branch WHIZARD1, but is no longer applicable in WHIZARD2. A WHIZARD-internal interface is foreseen for the very near future, most probably within this v2.2 release. Other options are to use parameterizations of the beam spectrum that have been included in the package CIRCE1 [6] which has been interfaced to WHIZARD since version v1.20 and been included in the WHIZARD2 release series. Another option is to generate a beam spectrum externally and then read it in as an ASCII data file, cf. Sec. 5.5.8. More about this can be found in a dedicated section on lepton collider spectra, Sec. 10.3.

In this section, we discuss the usage of beamstrahlung spectra by means of the CIRCE1 package. The beamstrahlung spectra are true spectra, so they have to be applied to pairs of beams, and an application to only one beam is meaningless. They are switched on by this beams statement including structure functions:

beams = e1, E1 => circe1

It is important to note that the parameterization of the beamstrahlung spectra within CIRCE1 contain also processes where e→γ conversions have been taking place, i.e. also hard processes with one (or two) initial photons can be simulated with beamstrahlung switched on. In that case, the explicit photon flags, ?circe1_photon1 and ?circe1_photon2, for the two beams have to be properly set, e.g. (ordering in the final state does not play a role):

process proc1 = A, e1 => A, e1
sqrts = 500 GeV
beams = e1, E1 => circe1
?circe1_photon1 = true
integrate (proc1)

process proc2 = e1, A => A, e1
sqrts = 1000 GeV
beams = e1, A => circe1
?circe1_photon2 = true

or

process proc1 = A, A => Wp, Wm
sqrts = 200 GeV
beams = e1, E1 => circe1
?circe1_photon1 = true
?circe1_photon2 = true
?circe1_generate = false

In all cases (one or both beams with photon conversion) the beam spectrum applies to both beams simultaneously.

In the last example (γγ→ W+W) the default CIRCE1 generator mode was turned off by unsetting ?circe1_generate. In the other examples this flag is set, by default. For standard use cases, CIRCE1 implements a beam-event generator inside the WHIZARD generator, which provides beam-event samples with correctly distributed probability. For electrons, the beamstrahlung spectrum sharply peaks near maximum energy. This distribution is most efficiently handled by the generator mode. By contrast, in the γγ mode, the beam-event c.m. energy is concentrated at low values. For final states with low invariant mass, which are typically produced by beamstrahlung photons, the generator mode is appropriate. However, the W+W system requires substantial energy, and such events will be very rare in the beam-event sample. Switching off the CIRCE1 generator mode solves this problem.

This is an overview over all options and flags for the CIRCE1 setup for lepton collider beamstrahlung:


ParameterDefaultMeaning
?circe1_photon1falsee→γ conversion for beam 1
?circe1_photon2falsee→γ conversion for beam 2
circe1_sqrtsscollider energy for the beam spectrum
?circe1_generatetrueflag for the CIRCE1 generator mode
?circe1_maptrueflag to apply special phase-space mapping
circe1_mapping_slope2.value of PS mapping exponent
circe1_eps1E-5parameter for mapping of spectrum peak position
circe1_ver0internal version of CIRCE1 package
circe1_rev0/most recentinternal revision of CIRCE1
$circe1_accSBANDaccelerator type
circe1_chat0chattiness/verbosity of CIRCE1

The collider energy relevant for the beamstrahlung spectrum is set by circe1_sqrts. As a default, this is always the value of sqrts set in the SINDARIN script. However, sometimes these values do not match, e.g. the user wants to simulate tt h at sqrts = 550 GeV, but the only available beam spectrum is for 500 GeV. In that case, circe1_sqrts = 500 GeV has to be set to use the closest possible available beam spectrum.

As mentioned in the discussion of the examples above, in CIRCE1 there are two options to use the beam spectra for beamstrahlung: intrinsic semi-analytic approximation formulae for the spectra, or a Monte-Carlo sampling of the sampling. The second possibility always give a better description of the spectra, and is the default for WHIZARD. It can, however, be switched off by setting the flag ?circe1_generate to false.

As the beamstrahlung spectra are sharply peaked at the collider energy, but still having long tails, a mapping of the spectra for an efficient phase-space sampling is almost mandatory. This is the default in WHIZARD, which can be changed by the flag ?circe1_map. Also, the default exponent for the mapping can be changed from its default value 2. with the variable circe1_mapping_slope. It is important to efficiently sample the peak position of the spectrum; the effective ratio of the peak to the whole sampling interval can be set by the parameter circe1_eps. The integer parameter circe1_chat sets the chattiness or verbosity of the CIRCE1 package, i.e. how many messages and warnings from the beamstrahlung generation/sampling will be issued.

The actual internal version and revision of the CIRCE1 package are set by the two integer parameters circe1_ver and circe1_rev. The default is in any case always the newest version and revision, while older versions are still kept for backwards compatibility and regression testing.

Finally, the geometry and design of the accelerator type is set with the string variable $circe1_acc: it contains the possible options for the old "SBAND" and "XBAND" setups, as well as the "TESLA" and JLC/NLC SLAC design "JLCNLC". The setups for the most important energies of the ILC as they are summarized in the ILC TDR [13, 14, 15, 16] are available as ILC. Beam spectra for the CLIC [18, 19, 20] linear collider are much more demanding to correctly simulate (due to the drive beam concept; only the low-energy modes where the drive beam is off can be simulated with the same setup as the abovementioned machines). Their setup will be supported soon in one of the upcoming WHIZARD versions within the CIRCE2 package.

An example of how to generate beamstrahlung spectra with the help of the package CIRCE2 (that is also a part of WHIZARD) is this:

process eemm = e1, E1 => e2, E2
sqrts = 500 GeV
beams = e1, E1 => circe2
$circe2_file = "ilc500.circe"
$circe2_design = "ILC"
?circe_polarized = false

Here, the ILC design is used for a beamstrahlung spectrum at 500 GeV nominal energy, with polarization averaged (hence, the setting of polarization to false). A list of all available options can be found in Sec. 5.5.13.

More technical details about the simulation of beamstrahlung spectra see the documented source code of the CIRCE1 package, as well as Chap. 9. In the next section, we discuss how to read in beam spectra from external files.

5.5.8  Beam events

As mentioned in the previous section, beamstrahlung is one of the crucial ingredients for a realistic simulation of linear lepton colliders. One option is to take a pre-generated beam spectrum for such a machine, and make it available for simulation within WHIZARD as an external ASCII data file. Such files basically contain only pairs of energy fractions of the nominal collider energy √s (x values). In WHIZARD they can be used in simulation with the following beams statement:

beams = e1, E1 => beam_events
$beam_events_file = "<beam_spectrum_file>"   

Note that beam spectra must always be pair spectra, i.e. they are automatically applied to both beam simultaneously. Beam spectra via external files are expected to reside in the current working directory. Alternatively, WHIZARD searches for them in the install directory of WHIZARD in share/beam-sim. There you can find an example file, uniform_spread_2.5%.dat for such a beam spectrum. The only possible parameter that can be set is the flag ?beam_events_warn_eof whose default is true. This triggers the issuing of a warning when the end of file of an external beam spectrum file is reached. In such a case, WHIZARD starts to reuse the same file again from the beginning. If the available data points in the beam events file are not big enough, this could result in an insufficient sampling of the beam spectrum.

5.5.9  Gaussian beam-energy spread

Real beams have a small energy spread. If beamstrahlung is small, the spread may be approximately described as Gaussian. As a replacement for the full simulation that underlies CIRCE2 spectra, it is possible to impose a Gaussian distributed beam energy, separately for each beam.

beams = e1, E1 => gaussian
gaussian_spread_1 = 0.1\%
gaussian_spread_2 = 0.2\%

(Note that the % sign means multiplication by 0.01, as it should.) The spread values are defined as the σ value of the Gaussian distribution, i.e., 2/3 of the events are within ± 1σ for each beam, respectively.

5.5.10  Equivalent photon approximation

The equivalent photon approximation (EPA) uses an on-shell approximation for the eeγ collinear splitting to allow the simulation of photon-induced backgrounds in lepton collider physics. The original concept is that of the Weizsäcker-Williams approximation [21, 22, 23]. This is a single-beam structure function that can be applied to both beams, or also to one beam only. Examples are:

beams = e1, E1 => epa  

or for a single beam:

beams = e1, p => epa, pdf_builtin

The last process allows the reaction of (quasi-) on-shell photons with protons.

In the following, we collect the parameters and flags that can be adjusted when using the EPA inside WHIZARD:


ParameterDefaultMeaning
epa_alpha0/intrinsicvalue of αQED for EPA
epa_x_min0.soft photon cutoff in x (mandatory)
epa_q_min0.minimal γ momentum transfer
epa_mass0/intrinsicmass of the radiating fermion (mandatory)
epa_e_max0/√supper cutoff for EPA
?epa_recoilfalseflag to switch on recoil/pT
?epa_keep_energyfalserecoil flag to conserve energy in splitting

The adjustable parameters are partially similar to the parameters in the QED initial-state radiation (ISR), cf. Sec. 5.5.6: the parameter epa_alpha sets the value of the electromagnetic coupling constant, αQED used in the EPA structure function. If not set, this is taken from the value inside the active physics model. The same is true for the mass of the particle that radiates the photon of the hard interaction, which can be reset by the user with the variable epa_mass. There are two dimensionful scale parameters, the minimal momentum transfer to the photon, epa_q_min, which must not be zero, and the upper energy cutoff for the EPA structure function, epa_e_max. The default for the latter value is the collider energy, √s, or the energy reduced by another structure function like e.g. beamstrahlung, √ŝ. Furthermore, there is a soft-photon regulator for the splitting function in x space, epa_x_min, which also has to be explicitly set different from zero. Hence, a minimal viable scenario that will be accepted by WHIZARD looks like this:

beams = e1, E1 => epa  
epa_q_min = 5 GeV
epa_x_min = 0.01

Finally, like the ISR case in Sec. 5.5.6, there is a flag to consider the recoil of the photon against the radiating electron by setting ?epa_recoil to true (default: false).

Though in principle processes like e+ ee+ e γ γ where the two photons have been created almost collinearly and then initiate a hard process could be described by exact matrix elements and exact kinematics. However, the numerical stability in the very far collinear kinematics is rather challenging, such that the use of the EPA is very often an acceptable trade-off between quality of the description on the one hand and numerical stability and speed on the other hand.

In the case, the EPA is set after a second structure function like a hadron collider PDF, there is a flavor summation over the quark constituents inside the proton, which are then the radiating fermions for the EPA. Here, the masses of all fermions have to be identical.

More about the physics of the equivalent photon approximation can be found in Chap. 9.

5.5.11  Effective W approximation

An approach similar to the equivalent photon approximation (EPA) discussed in the previous section Sec. 5.5.10, is the usage of a collinear splitting function for the radiation of massive electroweak vector bosons W/Z, the effective W approximation (EWA). It has been developed for the description of high-energy weak vector-boson fusion and scattering processes at hadron colliders, particularly the Superconducting Super-Collider (SSC). This was at a time when the simulation of 2→ 4 processes war still very challenging and 2→ 6 processes almost impossible, such that this approximation was the only viable solution for the simulation of processes like ppjjVV and subsequent decays of the bosons VW, Z.

Unlike the EPA, the EWA is much more involved as the structure functions do depend on the isospin of the radiating fermions, and are also different for transversal and longitudinal polarizations. Also, a truely collinear kinematics is never possible due to the finite W and Z boson masses, which start becoming more and more negligible for energies larger than the nominal LHC energy of 14 TeV.

Though in principle all processes for which the EWA might be applicable are technically feasible in WHIZARD to be generated also via full matrix elements, the EWA has been implemented in WHIZARD for testing purposes, backwards compatibility and comparison with older simulations. Like the EPA, it is a single-beam structure function that can be applied to one or both beams. We only give an example for both beams here, this is for a 3 TeV CLIC collider:

sqrts = 3 TeV
beams = e1, E1 => ewa

And this is for LHC or a higher-energy follow-up collider (which also shows the concatenation of the single-beam structure functions, applied to both beams consecutively, cf. Sec. 5.5.14:

sqrts = 14 TeV
beams = p, p => pdf_builtin => ewa

Again, we list all the options, parameters and flags that can be adapted for the EWA:


ParameterDefaultMeaning
ewa_x_min0.soft W/Z cutoff in x (mandatory)
ewa_mass0/intrinsicmass of the radiating fermion
ewa_pt_max0/√ŝupper cutoff for EWA
?ewa_recoilfalserecoil switch
?ewa_keep_energyfalseenergy conservation for recoil in splitting

First of all, all coupling constants are taken from the active physics model as they have to be consistent with electroweak gauge invariance. Like for EPA, there is a soft x cutoff for the ff V splitting, ewa_x_min, that has to be set different from zero by the user. Again, the mass of the radiating fermion can be set explicitly by the user; and, also again, the masses for the flavor sum of quarks after a PDF as radiators of the electroweak bosons have to be identical. Also for the EWA, there is an upper cutoff for the pT of the electroweak boson, that can be set via eta_pt_max. Indeed, the transversal W/Z structure function is logarithmically divergent in that variable. If it is not set by the user, it is estimated from √s and the splitting kinematics.

For the EWA, there is a flag to switch on a recoil for the electroweak boson against the radiating fermion, ?ewa_recoil. Note that this is an experimental feature that is not completely tested. In any case, the non-collinear kinematics violates 4-four momentum conservation, so there are two choices: either to conserve the energy (?ewa_keep_energy = true) or to conserve 3-momentum (?ewa_keep_energy = false). Momentum conservation for the kinematics is the default. This is due to the fact that for energy conservation, there will be a net total momentum in the event including the beam remnants (ISR/EPA/EWA radiated particles) that leeds to unexpected or unphysical features in the energy distributions of the beam remnants recoiling against the rest of the event.

More details about the physics can be found in Chap. 9.

5.5.12  Energy scans using structure functions

In WHIZARD, there is an implementation of a pair spectrum, energy_scan, that allows to scan the energy dependence of a cross section without actually scanning over the collider energies. Instead, only a single integration at the upper end of the scan interval over the process with an additional pair spectrum structure function performed. The structure function is chosen in such a way, that the distribution of x values of the energy scan pair spectrum translates in a plot over the energy of the final state in an energy scan from 0 to sqrts for the process under consideration.

The simplest example is the 1/s fall-off with the Z resonance in e+e → µ+ µ, where the syntax is very easy:

process eemm = e1, E1 => e2, E2
sqrts = 500 GeV
cuts = sqrts_hat > 50
beams = e1, E1 => energy_scan
integrate (eemm)

The value of sqrts = 500 GeV gives the upper limit for the scan, while the cut effectively let the scan start at 50 GeV. There are no adjustable parameters for this structure function. How to plot the invariant mass distribution of the final-state muon pair to show the energy scan over the cross section, will be explained in Sec. 5.9.

More details can be found in Chap. 9.

5.5.13  Photon collider spectra

One option that has been discussed as an alternative possibility for a high-energy linear lepton collider is to convert the electron and positron beam via Compton backscattering off intense laser beams into photon beams [24, 25, 26]. Naturally, due to the production of the photon beams and the inherent electron spectrum, the photon beams have a characteristic spectrum. The simulation of such spectra is possible within WHIZARD by means of the subpackage CIRCE2, which have been mentioned already in Sec. 5.5.7. It allows to give a much more elaborate description of a linear lepton collider environment than CIRCE1 (which, however, is not in all cases necessary, as the ILC beamspectra for electron/positrons can be perfectly well described with CIRCE1).

Here is a typical photon collider setup where we take a photon-initiated process:

process aaww = A, A => Wp, Wm  

beams = A, A => circe2
$circe2_file = "teslagg_500_polavg.circe"
$circe2_design = "TESLA/GG"
?circe2_polarized = false

Here, the photons are the initial states initiating the hard scattering. The structure function is circe2 which always is a pair spectrum. The list of available options are:


ParameterDefaultMeaning
?circe2_polarizedtruespectrum respects polarization info
$circe2_filename of beam spectrum data file
$circe2_design"*"collider design

The only logical flag ?circe2_polarized let WHIZARD know whether it should keep polarization information in the beam spectra or average over polarizations. Naturally, because of the Compton backscattering generation of the photons, photon spectra are always polarized. The collider design can be specified by the string variable $circe2_design, where the default setting "*" corresponds to the default of CIRCE2 (which is the TESLA 500 GeV machine as discussed in the TESLA Technical Design Report [27, 28]). Note that up to now there have not been any setups for a photon collider option for the modern linear collider concepts like ILC and CLIC. The string variable $circe2_file then allows to give the name of the file containing the actual beam spectrum; all files that ship with WHIZARD are stored in the directory circe2/share/data.

More details about the subpackage CIRCE2 and the physics it covers, can be found in its own manual and the chapter Chap. 9.

5.5.14  Concatenation of several structure functions

As has been shown already in Sec. 5.5.10 and Sec. 5.5.11, it is possible within WHIZARD to concatenate more than one structure function, irrespective of the fact, whether the structure functions are single-beam structure functions or pair spectra. One important thing is whether there is a phase-space mapping for these structure functions. Also, there are some combinations which do not make sense from the physics point of view, for example using lepton-collider ISR for protons, and then afterwards switching on PDFs. Such combinations will be vetoed by WHIZARD, and you will find an error message like (cf. also Sec. 4.3):

******************************************************************************
******************************************************************************
*** FATAL ERROR: Beam structure: [....] not supported 
******************************************************************************
******************************************************************************  

Common examples for the concatenation of structure functions are linear collider applications, where beamstrahlung (macroscopic electromagnetic beam-beam interactions) and electron QED initial-state radiation are both switched on:

beams = e1, E1 => circe1 => isr

Another possibility is the simulation of photon-induced backgrounds at ILC or CLIC, using beamstrahlung and equivalent photon approximation (EPA):

beams = e1, E1 => circe1 => epa  

or with beam events from a data file:

beams = e1, E1 => beam_events => isr

In hadron collider physics, parton distribution functions (PDFs) are basically always switched on, while afterwards the user could specify to use the effective W approximation (EWA) to simulate high-energy vector boson scattering:

sqrts = 100 TeV
beams = p, p => pdf_builtin => ewa  

Note that this last case involves a flavor sum over the five active quark (and anti-quark) species u, d, c, s, b in the proton, all of which act as radiators for the electroweak vector bosons in the EWA.

This would be an example with three structure functions:

beams = e1, E1 => circe1 => isr => epa

5.6  Polarization

5.6.1  Initial state polarization

WHIZARD supports polarizing the inital state fully or partially by assigning a nontrivial density matrix in helicity space. Initial state polarization requires a beam setup and is initialized by means of the beams_pol_density statement5:

beams_pol_density = @([<spin entries>]), @([<spin entries>])

The command beams_pol_fraction gives the degree of polarization of the two beams:

beams_pol_fraction = <degree beam 1>, <degree beam 2>

Both commands in the form written above apply to scattering processes, where the polarization of both beams must be specified. The beams_pol_density and beams_pol_fraction are possible with a single beam declaration if a decay process is considered, but only then.

While the syntax for the command beams_pol_fraction is pretty obvious, the syntax for the actual specification of the beam polarization is more intricate. We start with the polarization fraction: for each beam there is a real number between zero (unpolarized) and one (complete polarization) that can be specified either as a floating point number like 0.4 or with a percentage: 40 %. Note that the actual arithmetics is sometimes counterintuitive: 80 % left-handed electron polarization means that 80 % of the electron beam are polarized, 20 % are unpolarized, i.e. 20 % have half left- and half right-handed polarization each. Hence, 90 % of the electron beam is left-handed, 10 % is right-handed.

How does the specification of the polarization work? If there are no entries at all in the polarization constructor, @(), the beam is unpolarized, and the spin density matrix is proportional to the unit/identity matrix. Placing entries into the @() constructor follows the concept of sparse matrices, i.e. the entries that have been specified will be present, while the rest remains zero. Single numbers do specify entries for that particular helicity on the main diagonal of the spin density matrix, e.g. for an electron @(-1) means (100%) left-handed polarization. Different entries are separated by commas: @(1,-1) sets the two diagonal entries at positions (1,1) and (−1,−1) in the density matrix both equal to one. Two remarks are in order already here. First, note that you do not have to worry about the correct normalization of the spin density matrix, WHIZARD is taking care of this automatically. Second, in the screen output for the beam data, only those entries of the spin density matrix that have been specified by the user, will be displayed. If a beams_pol_fraction statement appears, other components will be non-zero, but might not be shown. E.g. ILC-like, 80 % polarization of the electrons, 30 % positron polarization will be specified like this for left-handed electrons and right-handed positrons:

beams = e1, E1
beams_pol_density = @(-1), @(+1)
beams_pol_fraction = 80%, 30%  

The screen output will be like this:

| ------------------------------------------------------------------------
| Beam structure: e-, e+
|   polarization (beam 1):
|     @(-1: -1: ( 1.000000000000E+00, 0.000000000000E+00))
|   polarization (beam 2):
|     @(+1: +1: ( 1.000000000000E+00, 0.000000000000E+00))
|   polarization degree = 0.8000000, 0.3000000
| Beam data (collision):
|   e-   (mass = 0.0000000E+00 GeV)  polarized
|   e+   (mass = 0.0000000E+00 GeV)  polarized  

But because of the fraction of unpolarized electrons and positrons, the spin density matrices for electrons and positrons are:

ρ(e) = diag 
0.10, 0.90 
      ρ(e+) = diag 
0.65, 0.35 
   ,

respectively. So, in general, only the entries due to the polarized fraction will be displayed on screen. We will come back to more examples below.

Again, the setting of a single entry, e.g. @(± m), which always sets the diagonal component (± m, ± m) of the spin density matrix equal to one. Here m can have the following values for the different spins (in parentheses are entries that exist only for massive particles):

Spin jParticle typepossible m values
0Scalar boson0
1/2Spinor+1, -1
1(Massive) Vector boson+1, (0), -1
3/2(Massive) Vectorspinor+2, (+1), (-1), -2
2(Massive) Tensor+2, (+1), (0), (-1), -2

Off-diagonal entries that are equal to one (up to the normalization) of the spin-density matrix can be specified simply by the position, namely: @(m:m, m). This would result in a spin density matrix with diagonal entry 1 for the position (m″, m″), and an entry of 1 for the off-diagonal position (m,m′).

Furthermore, entries in the density matrix different from 1 with a numerical value <val> can be specified, separated by another colon: @(m:m:<val>). Here, it does not matter whether m and m′ are different or not. For m = m′ also diagonal spin density matrix entries different from one can be specified. Note that because spin density matrices have to be Hermitian, only the entry (m,m′) has to be set, while the complex conjugate entry at the transposed position (m′,m) is set automatically by WHIZARD.

We will give some general density matrices now, and after that a few more definite examples. In the general setups below, we always give the expression for the spin density matrix only for one single beam.

  • Unpolarized:
    beams_pol_density = @()
    This has the same effect as not specifying any polarization at all and is the only constructor available for scalars and fermions declared as left- or right-handed (like the neutrino). Density matrix:
    ρ = 
    1
    |m|
    (|m|: particle multiplicity which is 2 for massless, 2j + 1 for massive particles).
  • Circular polarization:
    beams_pol_density = @(± j)     beams_pol_fraction = f
    A fraction f (parameter range f ∈ [0 ; 1]) of the particles are in the maximum / minimum helicity eigenstate ± j, the remainder is unpolarized. For spin 1/2 and massless particles of spin >0, only the maximal / minimal entries of the density matrix are populated, and the density matrix looks like this:
    ρ = diag


    1± f
    2
     , 0 , … , 0 ,
    1∓ f
    2



  • Longitudinal polarization (massive):
    beams_pol_density = @(0)     beams_pol_fraction = f
    We consider massive particles with maximal spin component j, a fraction f of which having longitudinal polarization, the remainder is unpolarized. Longitudinal polarization is (obviously) only available for massive bosons of spin >0. Again, the parameter range for the fraction is: f ∈ [0 ; 1]. The density matrix has the form:
    ρ = diag


    1−f
    |m|
     , … , 
    1−f
    |m|
     , 
    1+f 
    |m| − 1
    |m|
     , 
    1−f
    |m|
     ,  … , 
    1−f
    |m|



    (|m| = 2j+1 : particle multiplicity)
  • Transverse polarization (along an axis):
    beams_pol_density = @(j, -j, j:-j:exp(-I*phi))     beams_pol_fraction = f
    This so called transverse polarization is a polarization along an arbitrary direction in the xy plane, with φ=0 being the positive x direction and φ=90 the positive y direction. Note that the value of phi has either to be set inside the beam polarization expression explicitly or by a statement real phi = val degree before. A fraction f of the particles are polarized, the remainder is unpolarized. Note that, although this yields a valid density matrix for all particles with multiplicity >1 (in which the only the highest and lowest helicity states are populated), it is meaningful only for spin 1/2 particles and massless bosons of spin >0. The range of the parameters are: f ∈ [0 ; 1] and φ ∈ ℝ. This yields a density matrix:
    ρ =












      10
    f
    2
     eiφ     
      00 0      
      ⋮⋮ 
      0 00      
      
    f
    2
     eiφ
    01












    (for antiparticles, the matrix is conjugated).
  • Polarization along arbitrary axis (θ, φ):
    beams_pol_density = @(j:j:1-cos(theta), j:-j:sin(theta)*exp(-I*phi), -j:-j:1+cos(theta))           beams_pol_fraction = f
    This example describes polarization along an arbitrary axis in polar coordinates (polar axis in positive z direction, polar angle θ, azimuthal angle φ). A fraction f of the particles are polarized, the remainder is unpolarized. Note that, although axis polarization defines a valid density matrix for all particles with multiplicity >1, it is meaningful only for particles with spin 1/2. Valid ranges for the parameters are f ∈ [0 ; 1], θ ∈ ℝ, φ ∈ ℝ. The density matrix then has the form:
    ρ = 
    1
    2
    ·






      1 − fcosθ0fsinθ  eiφ     
      00 0      
      ⋮⋮ 
      0 00      
      fsinθ  eiφ01 + fcosθ






  • Diagonal density matrix:
    beams_pol_density = @(j:j:hj, j-1:j-1:hj−1, , -j:-j:hj)
    This defines an arbitrary diagonal density matrix with entries ρj,j , … , ρj,−j.
  • Arbitrary density matrix:
    beams_pol_density = @({m:m′:xm,m}):
    Here, {m:m′:xm,m} denotes a selection of entries at various positions somewhere in the spin density matrix. WHIZARD will check whether this is a valid spin density matrix, but it does e.g. not have to correspond to a pure state.

The beam polarization statements can be used both globally directly with the beams specification, or locally inside the integrate or simulate command. Some more specific examples are in order to show how initial state polarization works:

  • beams = A, A 
    beams_pol_density = @(+1),  @(1, -1, 1:-1:-I)
    
    This declares the initial state to be composed of two incoming photons, where the first photon is right-handed, and the second photon has transverse polarization in y direction.
  • beams = A, A 
    beams_pol_density = @(+1),  @(1, -1, 1:-1:-1)
    
    Same as before, but this time the second photon has transverse polarization in x direction.
  • beams = "W+"
    beams_pol\_density = @(0)
    
    This example sets up the decay of a longitudinal vector boson.
  • beams = E1, e1
    scan int hel_ep = (-1, 1) {
       scan int hel_em = (-1, 1) {
          beams_pol_density = @(hel_ep), @(hel_em)
          integrate (eeww)
       }
    }
    integrate (eeww)
    
    This example loops over the different positron and electron helicity combinations and calculates the respective integrals. The beams_pol_density statement is local to the scan loop(s) and, therefore, the last integrate calculates the unpolarized integral.

Although beam polarization should be straightforward to use, some pitfalls exist for the unwary:

  • Once beams_pol_density is set globally, it persists and is applied every time beams is executed (unless it is reset). In particular, this means that code like
    process wwaa = Wp, Wm => A, A
    process zee = Z => e1, E1
    
    sqrts = 200 GeV
    beams_pol_density = @(1, -1, 1:-1:-1), @()
    beams = Wp, Wm
    integrate (wwaa)
    beams = Z
    integrate (zee)
    beams_pol_density = @(0)
    
    will throw an error, because WHIZARD complains that the spin density matrix has the wrong dimensionality for the second (the decay) process. This kind of trap can be avoided be using beams_pol_density only locally in integrate or simulate statements.
  • On-the-fly integrations executed by simulate use the beam setup found at the point of execution. This implies that any polarization settings you have previously done affect the result of the integration.
  • The unstable command also requires integrals of the selected decay processes, and will compute them on-the-fly if they are unavailable. Here, a polarized integral is not meaningful at all. Therefore, this command ignores the current beam setting and issues a warning if a previous polarized integral is available; this will be discarded.

5.6.2  Final state polarization

Final state polarization is available in WHIZARD in the sense that the polarization of real final state particles can be retained when generating simulated events. In order for the polarization of a particle to be retained, it must be declared as polarized via the polarized statement

polarized particle [, particle, ...]

The effect of polarized can be reversed with the unpolarized statement which has the same syntax. For example,

polarized "W+", "W-", Z

will cause the polarization of all final state W and Z bosons to be retained, while

unpolarized "W+", "W-", Z

will reverse the effect and cause the polarization to be summed over again. Note that polarized and unpolarized are global statements which cannot be used locally as command arguments and if you use them e.g. in a loop, the effects will persist beyond the loop body. Also, a particle cannot be polarized and unstable at the same time (this restriction might be loosened in future versions of WHIZARD).

After toggling the polarization flag, the generation of polarized events can be requested by using the ?polarized_events option of the simulate command, e.g.

simulate (eeww) { ?polarized_events = true }

When simulate is run in this mode, helicity information for final state particles that have been toggled as polarized is written to the event file(s) (provided that polarization is supported by the selected event file format(s) ) and can also be accessed in the analysis by means of the Hel observable. For example, an analysis definition like

analysis =
  if (all Hel == -1 ["W+"] and all Hel == -1 ["W-"] ) then
    record cta_nn (eval cos (Theta) ["W+"]) endif;
  if (all Hel == -1 ["W+"] and all Hel ==  0 ["W-"] )
    then record cta_nl (eval cos (Theta) ["W+"]) endif

can be used to histogram the angular distribution for the production of polarized W pairs (obviously, the example would have to be extended to cover all possible helicity combinations). Note, however, that helicity information is not available in the integration step; therefore, it is not possible to use Hel as a cut observable.

While final state polarization is straightforward to use, there is a caveat when used in combination with flavor products. If a particle in a flavor product is defined as polarized, then all particles “originating” from the product will act as if they had been declared as polarized — their polarization will be recorded in the generated events. E.g., the example

process test = u:d, ubar:dbar => d:u, dbar:ubar, u, ubar

! insert compilation, cuts and integration here

polarized d, dbar
simulate (test) {?polarized_events = true}

will generate events including helicity information for all final state d and d quarks, but only for part of the final state u and u quarks. In this case, if you had wanted to keep the helicity information also for all u and u, you would have had to explicitely include them into the polarized statement.

5.7  Cross sections

Integrating matrix elements over phase space is the core of WHIZARD’s activities. For any process where we want the cross section, distributions, or event samples, the cross section has to be determined first. This is done by a doubly adaptive multi-channel Monte-Carlo integration. The integration, in turn, requires a phase-space setup, i.e., a collection of phase-space channels, which are mappings of the unit hypercube onto the complete space of multi-particle kinematics. This phase-space information is encoded in the file xxx.phs, where xxx is the process tag. WHIZARD generates the phase-space file on the fly and can reuse it in later integrations.

For each phase-space channel, the unit hypercube is binned in each dimension. The bin boundaries are allowed to move during a sequence of iterations, each with a fixed number of sampled phase-space points, so they adapt to the actual phase-space density as far as possible. In addition to this intrinsic adaptation, the relative channel weights are also allowed to vary.

All these steps are done automatically when the integrate command is executed. At the end of the iterative adaptation procedure, the program has obtained an estimate for the integral of the matrix element over phase space, together with an error estimate, and a set of integration grids which contains all information on channel weights and bin boundaries. This information is stored in a file xxx.vg, where xxx is the process tag, and is used for event generation by the simulate command.

5.7.1  Integration

Since everything can be handled automatically using default parameters, it often suffices to write the command

  integrate (proc1)

for integrating the process with name tag proc1, and similarly

  integrate (proc1, proc2, proc3)

for integrating several processes consecutively. Options to the integrate command are specified, if not globally, by a local option string

  integrate (proc1, proc2, proc3) { mH = 200 GeV }

(It is possible to place a beams statement inside the option string, if desired.)

If the process is configured but not compiled, compilation will be done automatically. If it is not available at all, integration will fail.

The integration method can be specified by the string variable

$integration_method = "<method>"

The default method is called "vamp" and uses the VAMP algorithm and code. (At the moment, there is only a single simplistic alternative, using the midpoint rule or rectangle method for integration, "midpoint". This is mainly for testing purposes. In future versions of WHIZARD, more methods like e.g. Gauss integration will be made available). VAMP, however, is clearly the main integration method. It is done in several passes (usually two), and each pass consists of several iterations. An iteration consists of a definite number of calls to the matrix-element function.

For each iteration, WHIZARD computes an estimate of the integral and an estimate of the error, based on the binned sums of matrix element values and squares. It also computes an estimate of the rejection efficiency for generating unweighted events, i.e., the ratio of the average sampling function value over the maximum value of this function.

After each iteration, both the integration grids (the binnings) and the relative weights of the integration channels can be adapted to minimize the variance estimate of the integral. After each pass of several iterations, WHIZARD computes an average of the iterations within the pass, the corresponding error estimate, and a χ2 value. The integral, error, efficiency and χ2 value computed for the most recent integration pass, together with the most recent integration grid, are used for any subsequent calculation that involves this process, in particular for event generation.

In the default setup, during the first pass(es) both grid binnings and channel weights are adapted. In the final (usually second) pass, only binnings are further adapted. Roughly speaking, the final pass is the actual calculation, while the previous pass(es) are used for “warming up” the integration grids, without using the numerical results. Below, in the section about the specification of the iterations, Sec. 5.7.3, we will explain how it is possible to change the behavior of adapting grids and weights.

Here is an example of the integration output, which illustrates these properties. The SINDARIN script describes the process e+eqq qq with q being any light quark, i.e., W+W and ZZ production and hadronic decay together will any irreducible background. We cut on pT and energy of jets, and on the invariant mass of jet pairs. Here is the script:

alias q = d:u:s:c
alias Q = D:U:S:C
process proc_4f = e1, E1 => q, Q, q, Q

ms = 0  mc = 0
sqrts = 500 GeV
cuts = all (Pt > 10 GeV and E > 10 GeV) [q:Q] 
   and all M > 10 GeV [q:Q, q:Q]

integrate (proc_4f)

After the run is finished, the integration output looks like

| Process library 'default_lib': loading
| Process library 'default_lib': ... success.
| Integrate: compilation done
| RNG: Initializing TAO random-number generator
| RNG: Setting seed for random-number generator to 12511
| Initializing integration for process proc_4f:
| ------------------------------------------------------------------------
| Process [scattering]: 'proc_4f'
|   Library name  = 'default_lib'
|   Process index = 1
|   Process components:
|     1: 'proc_4f_i1':   e-, e+ => d:u:s:c, dbar:ubar:sbar:cbar,
|                                  d:u:s:c, dbar:ubar:sbar:cbar [omega]
| ------------------------------------------------------------------------
| Beam structure: [any particles]
| Beam data (collision):
|   e-  (mass = 5.1099700E-04 GeV)
|   e+  (mass = 5.1099700E-04 GeV)
|   sqrts = 5.000000000000E+02 GeV
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'proc_4f_i1.phs'
| Phase space: 123 channels, 8 dimensions
| Phase space: found 123 channels, collected in 15 groves.
| Phase space: Using 195 equivalences between channels.
| Phase space: wood
| Applying user-defined cuts.
| OpenMP: Using 8 threads
| Starting integration for process 'proc_4f'
| Integrate: iterations not specified, using default
| Integrate: iterations = 10:10000:"gw", 5:20000:""
| Integrator: 15 chains, 123 channels, 8 dimensions
| Integrator: Using VAMP channel equivalences
| Integrator: 10000 initial calls, 20 bins, stratified = T
| Integrator: VAMP
|=============================================================================|
| It      Calls  Integral[fb]  Error[fb]   Err[%]    Acc  Eff[%]   Chi2 N[It] |
|=============================================================================|
   1       9963  2.3797857E+03  3.37E+02   14.15   14.13*   4.02
   2       9887  2.8307603E+03  9.58E+01    3.39    3.37*   4.31
   3       9815  3.0132091E+03  5.10E+01    1.69    1.68*   8.37
   4       9754  2.9314937E+03  3.64E+01    1.24    1.23*  10.65
   5       9704  2.9088284E+03  3.40E+01    1.17    1.15*  12.99
   6       9639  2.9725788E+03  3.53E+01    1.19    1.17   15.34
   7       9583  2.9812484E+03  3.10E+01    1.04    1.02*  17.97
   8       9521  2.9295139E+03  2.88E+01    0.98    0.96*  22.27
   9       9435  2.9749262E+03  2.94E+01    0.99    0.96   20.25
  10       9376  2.9563369E+03  3.01E+01    1.02    0.99   21.10
|-----------------------------------------------------------------------------|
  10      96677  2.9525019E+03  1.16E+01    0.39    1.22   21.10    1.15  10
|-----------------------------------------------------------------------------|
  11      19945  2.9599072E+03  2.13E+01    0.72    1.02   15.03
  12      19945  2.9367733E+03  1.99E+01    0.68    0.96*  12.68
  13      19945  2.9487747E+03  2.03E+01    0.69    0.97   11.63
  14      19945  2.9777794E+03  2.03E+01    0.68    0.96*  11.19
  15      19945  2.9246612E+03  1.95E+01    0.67    0.94*  10.34
|-----------------------------------------------------------------------------|
  15      99725  2.9488622E+03  9.04E+00    0.31    0.97   10.34    1.05   5
|=============================================================================|
| Time estimate for generating 10000 events: 0d:00h:00m:51s
| Creating integration history display proc_4f-history.ps and proc_4f-history.pdf

Each row shows the index of a single iteration, the number of matrix element calls for that iteration, and the integral and error estimate. Note that the number of calls displayed are the real calls to the matrix elements after all cuts and possible rejections. The error should be viewed as the 1σ uncertainty, computed on a statistical


Figure 5.1: Graphical output of the convergence of the adaptation during the integration of a WHIZARD process.

basis. The next two columns display the error in percent, and the accuracy which is the same error normalized by √ncalls. The accuracy value has the property that it is independent of ncalls, it describes the quality of adaptation of the current grids. Good-quality grids have a number of order one, the smaller the better. The next column is the estimate for the rejection efficiency in percent. Here, the value should be as high as possible, with 100 % being the possible maximum.

In the example, the grids are adapted over ten iterations, after which the accuracy and efficiency have saturated at about 1.0 and 10 %, respectively. The asterisk in the accuracy column marks those iterations where an improvement over the previous iteration is seen. The average over these iterations exhibits an accuracy of 1.22, corresponding to 0.39 % error, and a χ2 value of 1.15, which is just right: apparently, the phase-space for this process and set of cuts is well-behaved. The subsequent five iterations are used for obtaining the final integral, which has an accuracy below one (error 0.3 %), while the efficiency settles at about 10 %. In this example, the final χ2 value happens to be quite small, i.e., the individual results are closer together than the error estimates would suggest. One should nevertheless not scale down the error, but rather scale it up if the χ2 result happens to be much larger than unity: this often indicates sub-optimally adapted grids, which insufficiently map some corner of phase space.

One should note that all values are subject to statistical fluctuations, since the number of calls within each iterations is finite. Typically, fluctuations in the efficiency estimate are considerably larger than fluctuations in the error/accuracy estimate. Two subsequent runs of the same script should yield statistically independent results which may differ in all quantities, within the error estimates, since the seed of the random-number generator will differ by default.

It is possible to get exactly reproducible results by setting the random-number seed explicitly, e.g.,

  seed = 12345  

at any point in the SINDARIN script. seed is a predefined intrinsic variable. The value can be any 32bit integer. Two runs with different seeds can be safely taken as statistically independent. In the example above, no seed has been set, and the seed has therefore been determined internally by WHIZARD from the system clock.

The concluding line with the time estimate applies to a subsequent simulation step with unweighted events, which is not actually requested in the current example. It is based on the timing and efficiency estimate of the most recent iteration.

As a default, a graphical output of the integration history will be produced (if both LATEX and MetaPost have been available during configuration). Fig. 5.1 shows how this looks like, and demonstrates how a proper convergence of the integral during the adaptation looks like. The generation of these graphical history files can be switched off using the command ?vis_history = false.

5.7.2  Integration run IDs

A single SINDARIN script may contain multiple calls to the integrate command with different parameters. By default, files generated for the same process in a subsequent integration will overwrite the previous ones. This is undesirable when the script is re-run: all results that have been overwritten have to be recreated.

To avoid this, the user may identify a specific run by a string-valued ID, e.g.

  integrate (foo) { $run_id = "first" }

This ID will become part of the file name for all files that are created specifically for this run. Often it is useful to create a run ID from a numerical value using sprintf, e.g., in this scan:

  scan real mh = (100 => 200 /+ 10) {
    $run_id = sprintf "%e" (mh)
    integrate (h_production)
  }

With unique run IDs, a subsequent run of the same SINDARIN script will be able to reuse all previous results, even if there is more than a single integration per process.

5.7.3  Controlling iterations

WHIZARD has some predefined numbers of iterations and calls for the first and second integration pass, respectively, which depend on the number of initial and final-state particles. They are guesses for values that yield good-quality grids and error values in standard situations, where no exceptionally strong peaks or loose cuts are present in the integrand. Actually, the large number of warmup iterations in the previous example indicates some safety margin in that respect.

It is possible, and often advisable, to adjust the iteration and call numbers to the particular situation. One may reduce the default numbers to short-cut the integration, if either less accuracy is needed, or CPU time is to be saved. Otherwise, if convergence is bad, the number of iterations or calls might be increased.

To set iterations manually, there is the iterations command:

  iterations = 5:50000, 3:100000  

This is a comma-separated list. Each pair of values corresponds to an integration pass. The value before the colon is the number of iterations for this pass, the other number is the number of calls per iteration.

While the default number of passes is two (one for warmup, one for the final result), you may specify a single pass

  iterations = 5:100000  

where the relative channel weights will not be adjusted (because this is the final pass). This is appropriate for well-behaved integrands where weight adaptation is not necessary.

You can also define more than two passes. That might be useful when reusing a previous grid file with insufficient quality: specify the previous passes as-is, so the previous results will be read in, and then a new pass for further adaptation.

In the final pass, the default behavior is to not adapt grids and weights anymore. Otherwise, different iterations would be correlated, and a final reliable error estimate would not be possible. For all but the final passes, the user can decide whether to adapt grids and weights by attaching a string specifier to the number of iterations: "g" does adapt grids, but not weights, "w" the other way round. "gw" or "wg" does adapt both. By the setting "", all adaptations are switched off. An example looks like this:

  iterations = 2:10000:"gw", 3:5000

Since it is often not known beforehand how many iterations the grid adaptation will need, it is generally a good idea to give the first pass a large number of iterations. However, in many cases these turn out to be not necessary. To shortcut iterations, you can set any of

accuracy_goal
error_goal
relative_error_goal

to a positive value. If this is done, WHIZARD will skip warmup iterations once all of the specified goals are reached by the current iteration. The final iterations (without weight adaptation) are always performed.

5.7.4  Phase space

Before integrate can start its work, it must have a phase-space configuration for the process at hand. The method for the phase-space parameterization is determined by the string variable $phs_method. At the moment there are only two options, "single", for testing purposes, that is mainly used internally, and WHIZARD’s traditional method, "wood". This parameterization is particularly adapted and fine-tuned for electroweak processes and might not be the ideal for for pure jet cross sections. In future versions of WHIZARD, more options for phase-space parameterizations will be made available, e.g. the RAMBO algorithm and its massive cousin, and phase-space parameterizations that take care of the dipole-like emission structure in collinear QCD (or QED) splittings. For the standard method, the phase-space parameterization is laid out in an ASCII file <process-name>_i<comp>.phs. Here, <process-name> is the process name chosen by the user while <comp> is the number of the process component of the corresponding process. This immediately shows that different components of processes are getting different phase space setups. This is necessary for inclusive processes, e.g. the sum of ppZ + nj and ppW + nj, or in future versions of WHIZARD for NLO processes, where one component is the interference between the virtual and the Born matrix element, and another one is the subtraction terms. Normally, you do not have to deal with this file, since WHIZARD will generate one automatically if it does not find one. (WHIZARD is careful to check for consistency of process definition and parameters before using an existing file.)

Experts might find it useful to generate a phase-space file and inspect and/or modify it before proceeding further. To this end, there is the parameter ?phs_only. If you set this true, WHIZARD skips the actual integration after the phase-space file has been generated. There is also a parameter ?vis_channels which can be set independently; if this is true, WHIZARD will generate a graphical visualization of the phase-space parameterizations encoded in the phase-space file. This file has to be taken with a grain of salt because phase space channels are represented by sample Feynman diagrams for the corresponding channel. This does however not mean that in the matrix element other Feynman diagrams are missing (the default matrix element method, O’Mega, is not using Feynman-diagrammatic amplitudes at all).

Things might go wrong with the default phase-space generation, or manual intervention might be necessary to improve later performance. There are a few parameters that control the algorithm of phase-space generation. To understand their meaning, you should realize that phase-space parameterizations are modeled after (dominant) Feynman graphs for the current process.

The main phase space setup wood

For the main phase-space parameterization of WHIZARD, which is called "wood", there are many different parameters and flags that allow to tune and customize the phase-space setup for every certain process:

The parameter phs_off_shell controls the number of off-shell lines in those graphs, not counting s-channel resonances and logarithmically enhanced s- and t-channel lines. The default value is 2. Setting it to zero will drop everything that is not resonant or logarithmically enhanced. Increasing it will include more subdominant graphs. (WHIZARD increases the value automatically if the default value does not work.)

There is a similar parameter phs_t_channel which controls multiperipheral graphs in the parameterizations. The default value is 6, so graphs with up to 6 t/u-channel lines are considered. In particular cases, such as e+enγ, all graphs are multiperipheral, and for n>7 WHIZARD would find no parameterizations in the default setup. Increasing the value of phs_t_channel solves this problem. (This is presently not done automatically.)

There are two numerical parameters that describe whether particles are treated like massless particles in particular situations. The value of phs_threshold_s has the default value 50 GeV. Hence, W and Z are considered massive, while b quarks are considered massless. This categorization is used for deciding whether radiation of b quarks can lead to (nearly) singular behavior, i.e., logarithmic enhancement, in the infrared and collinear regions. If yes, logarithmic mappings are applied to phase space. Analogously, phs_threshold_t decides about potential t-channel singularities. Here, the default value is 100 GeV, so amplitudes with W and Z in the t-channel are considered as logarithmically enhanced. For a high-energy hadron collider of 40 or 100 TeV energy, also W and Z in s-channel like situations might be necessary to be considered massless.

Such logarithmic mappings need a dimensionful scale as parameter. There are three such scales, all with default value 10 GeV: phs_e_scale (energy), phs_m_scale (invariant mass), and phs_q_scale (momentum transfer). If cuts and/or masses are such that energies, invariant masses of particle pairs, and momentum transfer values below 10 GeV are excluded or suppressed, the values can be kept. In special cases they should be changed: for instance, if you want to describe γ*→µ+µ splitting well down to the muon mass, no cuts, you may set phs_m_scale = mmu. The convergence of the Monte-Carlo integration result will be considerably faster.

There are more flags. These and more details about the phase space parameterization will be described in Sec. 8.2.

5.7.5  Cuts

WHIZARD 2 does not apply default cuts to the integrand. Therefore, processes with massless particles in the initial, intermediate, or final states may not have a finite cross section. This fact will manifest itself in an integration that does not converge, or is unstable, or does not yield a reasonable error or reweighting efficiency even for very large numbers of iterations or calls per iterations. When doing any calculation, you should verify first that the result that you are going to compute is finite on physical grounds. If not, you have to apply cuts that make it finite.

A set of cuts is defined by the cuts statement. Here is an example

cuts = all Pt > 20 GeV [colored]  

This implies that events are kept only (for integration and simulation) if the transverse momenta of all colored particles are above 20 GeV.

Technically, cuts is a special object, which is unique within a given scope, and is defined by the logical expression on the right-hand side of the assignment. It may be defined in global scope, so it is applied to all subsequent processes. It may be redefined by another cuts statement. This overrides the first cuts setting: the cuts statement is not cumulative. Multiple cuts should be specified by the logical operators of SINDARIN, for instance

cuts = all Pt > 20 GeV [colored]  
  and all E > 5 GeV [photon]

Cuts may also be defined local to an integrate command, i.e., in the options in braces. They will apply only to the processes being integrated, overriding any global cuts.

The right-hand side expression in the cuts statement is evaluated at the point where it is used by an integrate command (which could be an implicit one called by simulate). Hence, if the logical expression contains parameters, such as

mH = 120 GeV
cuts = all M > mH [b, bbar]
mH = 150 GeV
integrate (myproc)

the Higgs mass value that is inserted is the value in place when integrate is evaluated, 150 GeV in this example. This same value will also be used when the process is called by a subsequent simulate; it is integrate which compiles the cut expression and stores it among the process data. This behavior allows for scanning over parameters without redefining the cuts every time.

The cut expression can make use of all variables and constructs that are defined at the point where it is evaluated. In particular, it can make use of the particle content and kinematics of the hard process, as in the example above. In addition to the predefined variables and those defined by the user, there are the following variables which depend on the hard process:

integer:n_in, n_out, n_tot
real:sqrts, sqrts_hat

Example:

cuts = sqrts_hat > 150 GeV

The constants n_in etc. are sometimes useful if a generic set of cuts is defined, which applies to various processes simultaneously.

The user is encouraged to define his/her own set of cuts, if possible in a process-independent manner, even if it is not required. The include command allows for storing a set of cuts in a separate SINDARIN script which may be read in anywhere. As an example, the system directories contain a file default_cuts.sin which may be invoked by

include ("default_cuts.sin")

5.7.6  QCD scale and coupling

WHIZARD treats all physical parameters of a model, the coefficients in the Lagrangian, as constants. As a leading-order program, WHIZARD does not make use of running parameters as they are described by renormalization theory. For electroweak interactions where the perturbative expansion is sufficiently well behaved, this is a consistent approach.

As far as QCD is concerned, this approach does not yield numerically reliable results, even on the validity scale of the tree approximation. In WHIZARD2, it is therefore possible to replace the fixed value of αs (which is accessible as the intrinsic model variable alphas), by a function of an energy scale µ.

This is controlled by the parameter ?alpha_s_is_fixed, which is true by default. Setting it to false enables running αs. The user has then to decide how αs is calculated.

One option is to set ?alpha_s_from_lhapdf (default false). This is recommended if the LHAPDF library is used for including structure functions, but it may also be set if LHAPDF is not invoked. WHIZARD will then use the αs formula and value that matches the active LHAPDF structure function set and member.

In the very same way, the αs running from the PDFs implemented intrinsically in WHIZARD can be taken by setting ?alpha_s_from_pdf_builtin to true. This is the same running then the one from LHAPDF, if the intrinsic PDF coincides with a PDF chosen from LHAPDF.

If this is not appropriate, there are again two possibilities. If ?alpha_s_from_mz is true, the user input value alphas is interpreted as the running value αs(mZ), and for the particular event, the coupling is evolved to the appropriate scale µ. The formula is controlled by the further parameters alpha_s_order (default 0, meaning leading-log; maximum 2) and alpha_s_nf (default 5).

Otherwise there is the option to set ?alpha_s_from_lambda_qcd = true in order to evaluate αs from the scale ΛQCD, represented by the intrinsic variable lambda_qcd. The reference value for the QCD scale is Λ_QCD = 200 MeV. alpha_s_order and alpha_s_nf apply analogously.

Note that for using one of the running options for αs, always ?alpha_s_is_fixed = false has to be invoked.

In any case, if αs is not fixed, each event has to be assigned an energy scale. By default, this is √ŝ, the partonic invariant mass of the event. This can be replaced by a user-defined scale, the special object scale. This is assigned and used just like the cuts object. The right-hand side is a real-valued expression. Here is an example:

scale = eval Pt [sort by -Pt [colored]]

This selects the pT value of the first entry in the list of colored particles sorted by decreasing pT, i.e., the pT of the hardest jet.

The scale definition is used not just for running αs (if enabled), but it is also the factorization scale for the LHAPDF structure functions.

These two values can be set differently by specifying factorization_scale for the scale at which the PDFs are evaluated. Analogously, there is a variable renormalization_scale that sets the scale value for the running αs. Whenever any of these two values is set, it supersedes the scale value.

Just like the cuts expression, the expressions for scale, factorization_scale and also renormalization_scale are evaluated at the point where it is read by an explicit or implicit integrate command.

5.7.7  Reweighting factor

It is possible to reweight the integrand by a user-defined function of the event kinematics. This is done by specifying a weight expression. Syntax and usage is exactly analogous to the scale expression. Example:

weight = eval (1 + cos (Theta) ^ 2) [lepton]

We should note that the phase-space setup is not aware of this reweighting, so in complicated cases you should not expect adaptation to achieve as accurate results as for plain cross sections.

Needless to say, the default weight is unity.

5.8  Events

After the cross section integral of a scattering process is known (or the partial-width integral of a decay process), WHIZARD can generate event samples. There are two limiting cases or modes of event generation:

  1. For a physics simulation, one needs unweighted events, so the probability of a process and a kinematical configuration in the event sample is given by its squared matrix element.
  2. Monte-Carlo integration yields weighted events, where the probability (without any grid adaptation) is uniformly distributed over phase space, while the weight of the event is given by its squared matrix element.

The choice of parameterizations and the iterative adaptation of the integration grids gradually shift the generation mode from option 2 to option 1, which obviously is preferred since it simulates the actual outcome of an experiment. Unfortunately, this adaptation is perfect only in trivial cases, such that the Monte-Carlo integration yields non-uniform probability still with weighted events. Unweighted events are obtained by rejection, i.e., accepting an event with a probability equal to its own weight divided by the maximal possible weight. Furthermore, the maximal weight is never precisely known, so this probability can only be estimated.

The default generation mode of WHIZARD is unweighted. This is controlled by the parameter ?unweighted with default value true. Unweighted events are easy to interpret and can be directly compared with experiment, if properly interfaced with detector simulation and analysis.

However, when applying rejection to generate unweighted events, the generator discards information, and for a single event it needs, on the average, 1/є calls, where the efficiency є is the ratio of the average weight over the maximal weight. If ?unweighted is false, all events are kept and assigned their respective weights in histograms or event files.

5.8.1  Simulation

The simulate command generates an event sample. The number of events can be set either by specifying the integer variable n_events, or by the real variable luminosity. (This holds for unweighted events. If weighted events are requested, the luminosity value is ignored.) The luminosity is measured in femtobarns, but other units can be used, too. Since the cross sections for the processes are known at that point, the number of events is determined as the luminosity multiplied by the cross section.

As usual, both parameters can be set either as global or as local parameters:

  n_events = 10000
  simulate (proc1)
  simulate (proc2, proc3) { luminosity = 100 / 1 pbarn }

In the second example, both n_events and luminosity are set. In that case, WHIZARD chooses whatever produces the larger number of events.

If more than one process is specified in the argument of simulate, events are distributed among the processes with fractions proportional to their cross section values. The processes are mixed randomly, as it would be the case for real data.

The raw event sample is written to a file which is named after the first process in the argument of simulate. If the process name is proc1, the file will be named proc1.evx. You can choose another basename by the string variable $sample. For instance,

  simulate (proc1) { n_events = 4000  $sample = "my_events" }

will produce an event file my_events.evx which contains 4000 events.

This event file is in a machine-dependent binary format, so it is not of immediate use. Its principal purpose is to serve as a cache: if you re-run the same script, before starting simulation, it will look for an existing event file that matches the input. If nothing has changed, it will find the file previously generated and read in the events, instead of generating them. Thus you can modify the analysis or any further steps without repeating the time-consuming task of generating a large event sample. If you change the number of events to generate, the program will make use of the existing event sample and generate further events only when it is used up. If necessary, you can suppress the writing/reading of the raw event file by the parameters ?write_raw and ?read_raw.

If you try to reuse an event file that has been written by a previous version of WHIZARD, you may run into an incompatibility, which will be detected as an error. If this happens, you may enforce a compatibility mode (also for writing) by setting $event_file_version to the appropriate version string, e.g., "2.0". Be aware that this may break some more recent features in the event analysis.

Generating an event sample can serve several purposes. First of all, it can be analyzed directly, by WHIZARD’s built-in capabilities, to produce tables, histograms, or calculate inclusive observables. The basic analysis features of WHIZARD are described below in Sec. 5.9. It can be written to an external file in a standard format that a human or an external program can understand. In Chap. 11, you will find a more thorough discussion of event generation with WHIZARD, which also covers in detail the available event-file formats. Finally, WHIZARD can rescan an existing event sample. The event sample may either be the result of a previous simulate run or, under certain conditions, an external event sample produced by another generator or reconstructed from data.

rescan "my_events" (proc1) { $pdf_builtin_set = "MSTW2008LO" }

The rescanning may apply different parameters and recalculate the matrix element, it may apply a different event selection, it may reweight the events by a different PDF set (as above). The modified event sample can again be analyzed or written to file. For more details, cf. Sec. 11.7.

5.8.2  Decays

Normally, the events generated by the simulate command will be identical in structure to the events that the integrate command generates. This implies that for a process such as ppW+W, the final-state particles are on-shell and stable, so they appear explicitly in the generated event files. If events are desired where the decay products of the W bosons appear, one has to generate another process, e.g., ppudūd. In this case, the intermediate vector bosons, if reconstructed, are off-shell as dictated by physics, and the process contains all intermediate states that are possible. In this example, the matrix element contains also ZZ, photon, and non-resonant intermediate states. (This can be restricted via the $restrictions option, cf. 5.4.3.

Another approach is to factorize the process in production (of W bosons) and decays (Wqq). This is actually the traditional approach, since it is much less computing-intensive. The factorization neglects all off-shell effects and irreducible background diagrams that do not have the decaying particles as an intermediate resonance. While WHIZARD is able to deal with multi-particle processes without factorization, the needed computing resources rapidly increase with the number of external particles. Particularly, it is the phase space integration that becomes the true bottleneck for a high multiplicity of final state particles.

In order to use the factorized approach, one has to specify particles as unstable. (Also, the ?allow_decays switch must be true; this is however its default value.) We give an example for a ppWj final state:

process wj  = u, gl => d, Wp
process wen = Wp => E1, n1

integrate (wen)

sqrts = 7 TeV
beams = p, p => pdf_builtin
unstable Wp (wen)
simulate (wj) { n_events = 1 }

This defines a 2 → 2 hard scattering process of W + j production at the 7 TeV LHC 2011 run. The W+ is marked as unstable, with its decay process being W+e+ νe. In the simulate command both processes, the production process wj and the decay process wen will be integrated, while the W decays become effective only in the final event sample. This event sample will contain final states with multiplicity 3, namely e+ νe d. Note that here only one decay process is given, hence the branching ratio for the decay will be taken to be 100 % by WHIZARD.

A natural restriction of the factorized approach is the implied narrow-width approximation. Theoretically, this restriction is necessary since whenever the width plays an important role, the usage of the factorized approach will not be fully justified. In particular, all involved matrix elements must be evaluated on-shell, or otherwise gauge-invariance issues could spoil the calculation. (There are plans for a future WHIZARD version to also include Breit-Wigner or Gaussian distributions when using the factorized approach.)

Decays can be concatenated, e.g. for top pair production and decay, e+ et t with decay tW+ b, and subsequent leptonic decay of the W as in W+ → µ+ νµ:

process eett = e1, E1 => t, tbar
process t_dec = t => Wp, b
process W_dec = Wp => E2, n2

unstable t (t_dec)
unstable Wp (W_dec)

sqrts = 500
simulate (eett) { n_events = 1 }

Note that in this case the final state in the event file will consist of t b µ+ νµ because the anti-top is not decayed.

If more than one decay process is being specified like in

  process eeww = e1, E1 => Wp, Wm
  process w_dec1 = Wp => E2, n2
  process w_dec2 = Wp => E3, n3

  unstable Wp (w_dec1, w_dec2)
  
  sqrts = 500
  simulate (eeww) { n_events = 100 }

then WHIZARD takes the integrals of the specified decay processes and distributes the decays statistically according to the calculated branching ratio. Note that this might not be the true branching ratios if decay processes are missing, or loop corrections to partial widths give large(r) deviations. In the calculation of the code above, WHIZARD will issue an output like

| Unstable particle W+: computed branching ratios:
|   w_dec1: 5.0018253E-01   mu+, numu
|   w_dec2: 4.9981747E-01   tau+, nutau
|   Total width = 4.5496085E-01 GeV (computed)
|               = 2.0490000E+00 GeV (preset)
|   Decay options: helicity treated exactly  

So in this case, WHIZARD uses 50 % muonic and 50 % tauonic decays of the positively charged W, while the W appears directly in the event file. WHIZARD shows the difference between the preset W width from the physics model file and the value computed from the two decay channels.

Note that a particle in a SINDARIN input script can be also explictly marked as being stable, using the

stable <particle-tag>  

constructor for the particle <particle-tag>.

Spin correlations in decays

By default, WHIZARD applies full spin and color correlations to the factorized processes, so it keeps both color and spin coherence between productions and decays. Correlations between decay products of distinct unstable particles in the same event are also fully retained. The program sums over all intermediate quantum numbers.

Although this approach obviously yields the optimal description with the limits of production-decay factorization, there is support for a simplified handling of particle decays. Essentially, there are four options, taking a decay W_ud: W→ ū d as an example:

  1. Full spin correlations: unstable Wp (W_ud)
  2. Isotropic decay: unstable Wp (W_ud) { ?isotropic_decay = true }
  3. Diagonal decay matrix: unstable Wp (W_ud) { ?diagonal_decay = true }
  4. Project onto specific helicity: unstable Wp (W_ud) { decay_helicity = -1 }

Here, the isotropic option completely eliminates spin correlations. The diagonal-decays option eliminates just the off-diagonal entries of the W spin-density matrix. This is equivalent to a measurement of spin before the decay. As a result, spin correlations are still present in the classical sense, while quantum coherence is lost. The definite-helicity option is similar and additional selects only the specified helicity component for the decaying particle, so its decay distribution assumes the shape for an accordingly polarized particle. All options apply in the rest frame of the decaying particle, with the particle’s momentum as the quantization axis.

Automatic decays

A convenient option is if the user did not have to specify the decay mode by hand, but if they were generated automatically. WHIZARD does have this option: the flag ?auto_decays can be set to true, and is taking care of that. In that case the list for the decay processes of the particle marked as unstable is left empty (we take a W again as example):

unstable Wm () { ?auto_decays = true }  

WHIZARD then inspects at the local position within the SINDARIN input file where that unstable statement appears the masses of all the particles of the active physics model in order to determine which decays are possible. It then calculates their partial widths. There are a few options to customize the decays. The integer variable auto_decays_multiplicity allows to set the maximal multiplicity of the final states considered in the auto decay option. The defaul value of that variable is 2; please be quite careful when setting this to values larger than that. If you do so, the flag ?auto_decays_radiative allows to specify whether final states simply containing additional resolved gluons or photons are taken into account or not. For the example above, you almost hit the PDG value for the W total width:

| Unstable particle W-: computed branching ratios:
|   decay_a24_1: 3.3337068E-01   d, ubar
|   decay_a24_2: 3.3325864E-01   s, cbar
|   decay_a24_3: 1.1112356E-01   e-, nuebar
|   decay_a24_4: 1.1112356E-01   mu-, numubar
|   decay_a24_5: 1.1112356E-01   tau-, nutaubar
|   Total width = 2.0478471E+00 GeV (computed)
|               = 2.0490000E+00 GeV (preset)
|   Decay options: helicity treated exactly  

Future shorter notation for decays

In an upcoming WHIZARD version there will be a shorter and more concise notation already in the process definition for such decays, which, however, is current not yet implemented. The two first examples above will then be shorter and have this form:

    process wj = u, gl => (Wp => E1, n1), d
  

as well as

    process eett = e1, E1 => (t => (Wp => E2, n2), b), tbar
  

5.8.3  Event formats

As mentioned above, the internal WHIZARD event format is a machine-dependent event format. There are a series of human-readable ASCII event formats that are supported: very verbose formats intended for debugging, formats that have been agreed upon during the Les Houches workshops like LHA and LHEF, or formats that are steered through external packages like HepMC. More details about event formats can be found in Sec. 11.5.

5.9  Analysis and Visualization

SINDARIN natively supports basic methods of data analysis and visualization which are frequently used in high-energy physics studies. Data generated during script execution, in particular simulated event samples, can be analyzed to evaluate further observables, fill histograms, and draw two-dimensional plots.

So the user does not have to rely on his/her own external graphical analysis method (like e.g. gnuplot or ROOT etc.), but can use methods that automatically ship with WHIZARD. In many cases, the user, however, clearly will use his/her own analysis machinery, especially experimental collaborations.

In the following sections, we first summarize the available data structures, before we consider their graphical display.

5.9.1  Observables

Analyses in high-energy physics often involve averages of quantities other than a total cross section. SINDARIN supports this by its observable objects. An observable is a container that collects a single real-valued variable with a statistical distribution. It is declared by a command of the form

observable analysis-tag

where analysis-tag is an identifier that follows the same rules as a variable name.

Once the observable has been declared, it can be filled with values. This is done via the record command:

record analysis-tag (value)

To make use of this, after values have been filled, we want to perform the actual analysis and display the results. For an observable, these are the mean value and the standard deviation. There is the command write_analysis:

write_analysis (analysis-tag)

Here is an example:

observable obs
record obs (1.2)  record obs (1.3)  record obs (2.1)  record obs (1.4)
write_analysis (obs)

The result is displayed on screen:

###############################################################################
# Observable: obs
average     =  1.500000000000E+00
error[abs]  =  2.041241452319E-01
error[rel]  =  1.360827634880E-01
n_entries   = 4

5.9.2  The analysis expression

The most common application is the computation of event observables – for instance, a forward-backward asymmetry – during simulation. To this end, there is an analysis expression, which behaves very similar to the cuts expression. It is defined either globally

analysis = logical-expr

or as a local option to the simulate or rescan commands which generate and handle event samples. If this expression is defined, it is not evaluated immediately, but it is evaluated once for each event in the sample.

In contrast to the cuts expression, the logical value of the analysis expression is discarded; the expression form has been chosen just by analogy. To make this useful, there is a variant of the record command, namely a record function with exactly the same syntax. As an example, here is a calculation of the forward-backward symmetry in a process ee_mumu with final state µ+µ:

  observable a_fb
  analysis = record a_fb (eval sgn (Pz) ["mu-"])
  simulate (ee_mumu) { luminosity = 1 / 1 fbarn }

The logical return value of record – which is discarded here – is true if the recording was successful. In case of histograms (see below) it is true if the value falls within bounds, false otherwise.

Note that the function version of record can be used anywhere in expressions, not just in the analysis expression.

When record is called for an observable or histogram in simulation mode, the recorded value is weighted appropriately. If ?unweighted is true, the weight is unity, otherwise it is the event weight.

The analysis expression can involve any other construct that can be expressed as an expression in SINDARIN. For instance, this records the energy of the 4th hardest jet in a histogram pt_dist, if it is in the central region:

  analysis = 
    record pt_dist (eval E [extract index 4 
                             [sort by - Pt 
                               [select if -2.5 < Eta < 2.5 [colored]]]])

Here, if there is no 4th jet in the event which satisfies the criterion, the result will be an undefined value which is not recorded. In that case, record evaluates to false.

Selection cuts can be part of the analysis expression:

  analysis =
    if any Pt > 50 GeV [lepton] then
      record jet_energy (eval E [collect [jet]])
    endif

Alternatively, we can specify a separate selection expression:

  selection = any Pt > 50 GeV [lepton]
  analysis = record jet_energy (eval E [collect [jet]])

The former version writes all events to file (if requested), but applies the analysis expression only to the selected events. This allows for the simultaneous application of different selections to a single event sample. The latter version applies the selection to all events before they are analyzed or written to file.

The analysis expression can make use of all variables and constructs that are defined at the point where it is evaluated. In particular, it can make use of the particle content and kinematics of the hard process, as in the example above. In addition to the predefined variables and those defined by the user, there are the following variables which depend on the hard process. Some of them are constants, some vary event by event:

integer:event_index
integer:process_num_id
string:$process_id
integer:n_in, n_out, n_tot
real:sqrts, sqrts_hat
real:sqme, sqme_ref
real:event_weight, event_excess

The process_num_id is the numeric ID as used by external programs, while the process index refers to the current library. By default, the two are identical. The process index itself is not available as a predefined observable. The sqme and sqme_ref values indicate the squared matrix element and the reference squared matrix element, respectively. The latter applies when comparing with a reference sample (the rescan command).

record evaluates to a logical, so several record functions may be concatenated by the logical operators and or or. However, since usually the further evaluation should not depend on the return value of record, it is more advisable to concatenate them by the semicolon (;) operator. This is an operator (not a statement separator or terminator) that connects two logical expressions and evaluates both of them in order. The lhs result is discarded, the result is the value of the rhs:

  analysis = 
    record hist_pt (eval Pt [lepton]) ; record hist_ct (eval cos (Theta) [lepton])

5.9.3  Histograms

In SINDARIN, a histogram is declared by the command

histogram analysis-tag (lower-bound, upper-bound)

This creates a histogram data structure for an (unspecified) observable. The entries are organized in bins between the real values lower-bound and upper-bound. The number of bins is given by the value of the intrinsic integer variable n_bins, the default value is 20.

The histogram declaration supports an optional argument, so the number of bins can be set locally, for instance

histogram pt_distribution (0 GeV, 500 GeV) { n_bins = 50 }

Sometimes it is more convenient to set the bin width directly. This can be done in a third argument to the histogram command.

histogram pt_distribution (0 GeV, 500 GeV, 10 GeV)

If the bin width is specified this way, it overrides the setting of n_bins.

The record command or function fills histograms. A single call

record (real-expr)

puts the value of real-expr into the appropriate bin. If the call is issued during a simulation where unweighted is false, the entry is weighted appropriately.

If the value is outside the range specified in the histogram declaration, it is put into one of the special underflow and overflow bins.

The write_analysis command prints the histogram contents as a table in blank-separated fixed columns. The columns are: x (bin midpoint), y (bin contents), Δ y (error), excess weight, and n (number of entries). The output also contains comments initiated by a # sign, and following the histogram proper, information about underflow and overflow as well as overall contents is added.

5.9.4  Plots

While a histogram stores only summary information about a data set, a plot stores all data as (x,y) pairs, optionally with errors. A plot declaration is as simple as

plot analysis-tag

Like observables and histograms, plots are filled by the record command or expression. To this end, it can take two arguments,

record (x-expr, y-expr)

or up to four:

record (x-expr, y-expr, y-error)
record (x-expr, y-expr, y-error-expr, x-error-expr)

Note that the y error comes first. This is because applications will demand errors for the y value much more often than x errors.

The plot output, again written by write_analysis contains the four values for each point, again in the ordering x,yy, Δ x.

5.9.5  Analysis Output

There is a default format for piping information into observables, histograms, and plots. In older versions of WHIZARD there was a first version of a custom format, which was however rather limited. A more versatile custom output format will be coming soon.

  1. By default, the write_analysis command prints all data to the standard output. The data are also written to a default file with the name whizard_analysis.dat. Output is redirected to a file with a different name if the variable $out_file has a nonempty value. If the file is already open, the output will be appended to the file, and it will be kept open. If the file is not open, write_analysis will open the output file by itself, overwriting any previous file with the same name, and close it again after data have been written.

    The command is able to print more than one dataset, following the syntax

    write_analysis (analysis-tag1, analysis-tag2, …) { options }

    The argument in brackets may also be empty or absent; in this case, all currently existing datasets are printed.

    The default data format is suitable for compiling analysis data by WHIZARD’s built-in gamelan graphics driver (see below and particularly Chap. 13). Data are written in blank-separated fixed columns, headlines and comments are initiated by the # sign, and each data set is terminated by a blank line. However, external programs often require special formatting.

    The internal graphics driver gamelan of WHIZARD is initiated by the compile_analysis command. Its syntax is the same, and it contains the write_analysis if that has not been separately called (which is unnecessary). For more details about the gamelan graphics driver and data visualization within WHIZARD, confer Chap. 13.

  2. Custom format. Not yet (re-)implemented in a general form.

5.10  Custom Input/Output

WHIZARD is rather chatty. When you run examples or your own scripts, you will observe that the program echoes most operations (assignments, commands, etc.) on the standard output channel, i.e., on screen. Furthermore, all screen output is copied to a log file which by default is named whizard.log.

For each integration run, WHIZARD writes additional process-specific information to a file ⟨tag.log, where ⟨tag⟩ is the process name. Furthermore, the write_analysis command dumps analysis data – tables for histograms and plots – to its own set of files, cf. Sec. 5.9.

However, there is the occasional need to write data to extra files in a custom format. SINDARIN deals with that in terms of the following commands:

5.10.1  Output Files

open_out

open_out (filename)
open_out (
filename) { options }

Open an external file for writing. If the file exists, it is overwritten without warning, otherwise it is created. Example:

open_out ("my_output.dat")

close_out

close_out (filename)
close_out (
filename) { options }

Close an external file that is open for writing. Example:

close_out ("my_output.dat")

5.10.2  Printing Data

printf

printf format-string-expr
printf
format-string-expr (data-objects)

Format ⟨data-objects⟩ according to ⟨format-string-expr⟩ and print the resulting string to standard output if the string variable $out_file is undefined. If $out_file is defined and the file with this name is open for writing, print to this file instead.

Print a newline at the end if ?out_advance is true, otherwise don’t finish the line.

The ⟨format-string-expr⟩ must evaluate to a string. Formatting follows a subset of the rules for the printf(3) command in the C language. The supported rules are:

  • All characters are printed as-is, with the exception of embedded conversion specifications.
  • Conversion specifications are initiated by a percent (%) sign and followed by an optional prefix flag, an optional integer value, an optional dot followed by another integer, and a mandatory letter as the conversion specifier.
  • A percent sign immediately followed by another percent sign is interpreted as a single percent sign, not as a conversion specification.
  • The number of conversion specifiers must be equal to the number of data objects. The data types must also match.
  • The first integer indicates the minimum field width, the second one the precision. The field is expanded as needed.
  • The conversion specifiers d and i are equivalent, they indicate an integer value.
  • The conversion specifier e indicates a real value that should be printed in exponential notation.
  • The conversion specifier f indicates a real value that should be printed in decimal notation without exponent.
  • The conversion specifier g indicates a real value that should be printed either in exponential or in decimal notation, depending on its value.
  • The conversion specifier s indicates a logical or string value that should be printed as a string.
  • Possible prefixes are # (alternate form, mandatory decimal point for reals), 0 (zero padding), - (left adjusted), + (always print sign), ‘ ’ (print space before a positive number).

For more details, consult the printf(3) manpage. Note that other conversions are not supported and will be rejected by WHIZARD.

The data arguments are numeric, logical or string variables or expressions. Numeric expressions must be enclosed in parantheses. Logical expressions must be enclosed in parantheses prefixed by a question mark ?. String expressions must be enclosed in parantheses prefixed by a dollar sign $. These forms behave as anonymous variables.

Note that for simply printing a text string, you may call printf with just a format string and no data arguments.

Examples:

printf "The W mass is %8f GeV" (mW)

int i = 2
int j = 3
printf "%i + %i = %i" (i, j, (i+j))

string $directory = "/usr/local/share"
string $file = "foo.dat"
printf "File path: %s/%s" ($directory, $file)

There is a related sprintf function, cf. Sec. 5.1.5.

5.11  WHIZARD next-to-leading order mode

5.11.1  Prerequisites

A full NLO computation requires virtual matrix elements obtained from loop diagrams. Since O’Mega cannot calculate such diagrams, external programs have to be used. WHIZARD has a generic interface to each matrix-element generator which is BLHA-compatible. Explicit implementations exist for Gosam  OpenLoopsand Recola.

Setting up Gosam

The installation of Gosam is detailed on the HepForge page http://gosam/hepforge.org. We mention here some of the steps necessary to get it to be linked with WHIZARD.

Bug in Gosam installation scripts: In many versions of Gosam there is a bug in the installation scripts that is only effective if Gosam is installed with superuser privileges. Then all files in $installdir/share/golem do not have read privileges for normal users. These privileges must be given manually to all files in that directory.

Prerequisites for Gosam to produce code for one-loop matrix elements are the scientific algebra program form and the generator of loop topologies and diagrams, qgraf. These can be accessed via their respective webpages http://www.nikhef.nl/~form/ and http://cfif.ist.utl.pt/~paulo/qgraf.html. Note also that both Java and the Java runtime environment have to be installed in order for Gosam to properly work. Furthermore, libtool needs to be installed.

Setting up OpenLoops

The installation of OpenLoops is explained in detail on the HepForge page http://openloops.hepforge.org. In the following the main steps for usage with WHIZARD are summarized.

Please note that at the moment, OpenLoops cannot be installed such that in almost all cases the explicit OpenLoops package directory has to be set via --with-openloops=<openloops_dir>.

OpenLoops can be checked out using the URL

   \url{http://openloops.hepforge.org/svn/OpenLoops/branches/public}.

The program can simply be installed by running scons, of which a local version is included in the OpenLoopsdirectory. This produces the binary file OpenLoops, which is the main hook for all other usage of the program.

OpenLoops works by loading pre-compiled external libraries, which have to be installed for each individual process. This requires the file openloops.cfg, which should contain the following content:

   [OpenLoops]
   process\_repositories=public, whizard\\
   compile\_extra=1

The first line instructs OpenLoopsto also look for process libraries in an additional lepton collider repository. The second line triggers the inclusion of N+1-particle tree-level matrix elements in the process directory, so that a complete NLO calculation including real amplitudes can be performed only with OpenLoops.
The libraries can then be installed via

   ./openloops libinstall proc_name

A list of supported library names can be found on the OpenLoops web page. Note that a process library also includes all possible permutated processes, i.e. the process library ppll, for example, can also be used to compute the matrix elements for e+ eq q (massless quarks only). Additionally, there are process libraries for top and gauge boson decays, tbw, vjj, tbln and tbqq.

Finally, OpenLoops can be linked to WHIZARD during configuration by including

  --enable-openloops --with-openloops=$OPENLOOPS_PATH,

where $OPENLOOPS_PATH is the directory the OpenLoops binary is located in. OpenLoops can then be used with the SINDARINoption

  $loop_me_method = "openloops".

The functional tests which check the OpenLoopsfunctionality require the libraries ppll, eett and tbw to be installed (note that eett is not contained in ppll). During the configuration of WHIZARD, it is automatically checked that these two libraries, as well as the option compile_extra=1, are present.

OpenLoopsSindarin flags

Several sindarin options exist to control the behaviour of OpenLoops.

  • openloops_verbosity:
    Decide how much OpenLoopsoutput is printed. Can have values 0, 1 and 2.
  • ?openloops_use_cms:
    Activates the complex mass scheme
  • openloops_phs_tolerance:
  • openloops_switch_off_muon_yukawa:
    Sets the Yukawa coupling of muons to zero in order to assure agreement with O’Mega, which is possibly used for other components and per default does not take Hµµ couplings into account.
  • openloops_stability_log:
    Creates the directory stability_log containing information about the performance of the matrix elements. Possible values are
    • 0: No output (default),
    • 1: On finish() call,
    • 2: Adaptive,
    • 3: Always
  • openloops_top_signal: Only possible if the process is e+ et t. Loads a special library which only contains the signal diagram.
  • ?openloops_use_collier: Use Collier as the reduction method (default true).

Setting up Recola

The installation of Recola is explained in detail on the HepForge page http://recola.hepforge.org. In the following the main steps for usage with WHIZARD are summarized. The minimal required version number of Recolais 1.2.

Recola can be linked to WHIZARD during configuration by including

  --enable-recola --with-recola=$RECOLA_PATH,

where $RECOLA_PATH is the directory the Recola binary is located in. Recola can then be used with the SINDARINoption

  $loop_me_method = "recola".

5.11.2  NLO cross sections

NLO computation can be switched on in SINDARIN with

  process proc_nlo = in1, in2 => out1, ..., outN { nlo_calculation = "Components" },

where after the nlo_calculation can be followed by a list of strings specifying the desired NLO-components to be integrated, i.e. "Born", "Real", "Virtual", "Pdf", (for hadron collisions) or "Mismatch" (for the soft mismatch in resonance-aware computations) and "Full". The "Full"-option switches on all components and is required if the total NLO result is desired. For example, specifying

  nlo_calculation = "Born", "Virtual"

will result in the computation of the Born and virtual component.

The integration can be carried out in two different modes: Combined and separate integration. In the separate integration mode, each component is integrated individually, allowing for a good overview of their contributions to the total cross section. In the combined integration mode, all components are added up during integration so that the sum of them is evaluated. Here, only one integration pass will be displayed. The default method is the separate integration.

The convergence of the integration can crucially be influenced by the presence of resonances. A better convergence in this case is achieved activating the resonance-aware FKS subtraction,

  $fks_method = "resonances".

This mode comes with an additional integration component, the so-called soft mismatch.

5.11.3  Fixed-order NLO events

Fixed-order NLO events can be produced in two different modes: Combined and separated.

  • Combined:
    In the combined mode, one single integration grid is produced, from which events are generated with the total NLO weight. The corresponding event file contains N events with born-like kinematics and weight equal to B + V + ∑αr Cαr, where B is the Born matrix element, V is the virtual matrix element and Cαr are the subtraction terms in each singular region. For resonance-aware processes, also the mismatch value is added. Each born-like event is followed by Nphs associated events with real kinematics, i.e. events where one additional QCD particle is present. Nphs the number of distinct phase spaces. Two phase spaces are distinct if they share the same resonance history but have different emitters. So, two αr can share the same phase space index.
    The combined event mode is activated by
      ?combined_nlo_integration
      ?unweighted = false
      ?fixed_order_nlo_events = true
    
    Moreover, the process must be specified at next-to-leading-order in its definition using nlo_calculation = "Full". WHIZARDthen proceeds as in the usual simulation mode. I.e. it first checks whether integration grids are already present and uses them if they fit. Otherwise, it starts an integration.
  • Separete:
    In the separate mode, grids and events are generated for each individual component of the NLO process. This method is preferable for complicated processes, since it allows to individually tune each grid generation. Moreover, the grid generation then occurs trivially parallelized. The event files either contain only born-like kinematics with weight B or V (and mismatch in case of a resonance-aware process) or mixed born-like and real kinematics like in the combined mode. However, the born-like events have weight ∑αr Cαr.
    The separate event mode is activated by
      ?unweighted = false
      ?negative_weights = true
      ?fixed_order_nlo_events = true
    
    Note that negative weights have to be switched on because, in contrast to the combined mode, the total cross sections of the individual components can be negative.
    Also, the desired component has to appear in the process NLO specification, e.g. using nlo_calculation = "Real".

Fixed-order NLO events are so far only supported in the HepMC-format. The output can either be written to disk or put into a FIFO to interface it to an analysis program.
The weights in the real event output, both in the combined and separate mode, are multiplied with a factor Nphs + 1. This is to account for the fact that the number of events in the output file is artificially enhanced by this factor. Thus, if wi are the event weights, 1/Nreali wi, with Nreal = Nborn · (1 + Nphs) does not correctly reproduce the total cross section.

5.11.4  Powheg matching


1
In older versions of WHIZARD, until v2.1.1, there used to be separate comparators for the comparisons up to a tolerance, namely ==~ and <>~. These have been discarded from v2.2.0 on in order to simplify the syntax.
2
In older versions of WHIZARD, until v2.1.1, there also used to be a sprintd function and a printd command for default formats without a format string. They have been discarded in order to simplify the syntax from version v2.2.0 on.
3
In future versions of WHIZARD it is foreseen to implement other electroweak schemes.
4
Until WHIZARD version 2.2.1 including, only the LHAPDF series 5 was supported, while from version 2.2.2 on also the LHAPDF release series 6 has been supported.
5
Note that the syntax for the specification of beam polarization has changed from version v2.1 to v2.2 and is incompatible between the two release series. The old syntax beam_polarization with its different polarization constructors has been discarded in favor of a unified syntax.

Previous Up Next