Chapter 5 Detailed WHIZARD Steering: SINDARIN
5.1 Data and expressions5.1.1 Realvalued objectsReal literals have their usual form, mantissa and, optionally, exponent:
0. 3.14 .5 2.345e3 .890E023
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 hardcoded, 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: 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 machineprecision 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 realvalued intrinsic
variable tolerance. This variable is initially zero, but can be
set to any value (for instance, tolerance = 1.e13 by the user.
Note that for nonzero tolerance, operators like
5.1.2 Integervalued objectsInteger literals are obvious:
1 98765 0123
Integers are always signed. Their range is the defaultinteger 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 Complexvalued objectsComplex 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 Logicalvalued objectsThere 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 The names of logical variables begin with a question mark ?. Here is the declaration of a logical user variable:
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 Stringvalued objects and string operationsString literals are enclosed in double quotes: "This is a string."
The empty string is "". String variables begin with the dollar
sign:
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 formatstring
(arglist)
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 $(stringexpr). Example:
The related printf command with the same syntax prints the formatted string to standard output^{2}. 5.2 Particles and (sub)events5.2.1 Particle aliasesA particle species is denoted by its name as a string: An alias may either denote a single particle species or a class of particles species. A colon : concatenates particle names and aliases to yield multispecies aliases:
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 (singleparticle) aliases for all
particles it contains. Furthermore, it defines the class aliases
5.2.2 SubeventsSubevents 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:
These subevents evaluate to the set of all W^{+} bosons (to be precise, their fourmomenta), all u, d, or s quarks, and all colored particles, respectively. A subevent can contain pseudoparticles, i.e., particle combinations. That is, the fourmomenta of distinct particles are combined (added conmponentwise), and the results become subevent elements just like ordinary particles. The (pseudo)particles in a subevent are nonoverlapping. 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
In this expression, we first define 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 5.2.3 Subevent functionsThere are several functions that take a subevent (or an alias) as an argument and return a new subevent. Here we describe them: collectcollect [particles] First version: collect all particle momenta in the argument and combine them to a single fourmomentum. The particles argument may either be a subevt expression or an alias expression. The result is a oneentry 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. clusterNOTE: This is an experimental feature, available from version 2.2.1 on. cluster [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 oneentry 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 jetalgorithm values (integer constants): kt_algorithm and the variable jet_r to the desired R parameter value, as appropriate for the analysis and the jet algorithm. Example:
combinecombine [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:
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
selectselect if condition [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. extractextract [particles] Return a singleparticle subevent. In the first version, it contains the first particle in the subevent particles. In the second version, the particle with index indexvalue is returned, where indexvalue 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 welldefined, so you may wish to sort the subevent before applying the extract function to it. sortsort [particles] 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:
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 refparticle, if it has more than one. joinjoin [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
5.2.4 Calculating observablesObservables (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.: evaleval expr [particles] 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,
evaluates to the transverse momentum of the first colored particle,
evaluates to the invariant mass of the first distinct pair of jets (assuming
that
evaluates to the difference of energy and mass of the combination of the first electronneutrino 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 hardcoding them as new functions. 5.2.5 Cuts and event selectionInstead of a numeric value, we can use observables to compute a logical value. allall logical_expr [particles] The all construct expects a logical expression and one or two subevent arguments in square brackets.
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 twoargument case it must be true for all nonoverlapping combinations of particles in the two subevents. If one of the arguments is the empty subevent, the result is also true. anyany logical_expr [particles] The any construct is true if the logical expression is true for at least one particle or nonoverlapping particle combination:
This defines a trigger or selection condition. If a subevent argument is empty, it evaluates to false nono logical_expr [particles] The no construct is true if the logical expression is true for no single one particle or nonoverlapping particle combination:
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 functionscountcount [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 nonoverlapping particle pairs (that pass the test, if any). Predefined observablesThe following realvalued 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.
There is also an integervalued observable:
5.3 Physics ModelsA 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 phasespace 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 nomodel 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:
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 G_{F}m_{W}m_{Z} scheme. (While this scheme is unusual for loop calculations, it is natural for a treelevel event generator where the Z and W poles have to be at their experimentally determined location^{3}.) 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:
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
After the model is reassigned, you will see the MSSM value of m_{b} which still has its default value, not the one you have given. However, if you revert to the SM later,
you will see that your modification of the SM’s m_{b} 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 pureSM 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 ProcessesThe purpose of WHIZARD is the integration and simulation of highenergy 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 definitionA process object is defined in a straightforward notation. The definition syntax is straightforward:
process processid = incomingparticles Here are typical examples:
Throughout the program, the process will be identified by its processid, 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 uppercase 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 matrixelement 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 leadingorder 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 fieldtheoretical 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 namesThe 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
This tells that you can identify an electron either as 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.
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
unless they deal with processes where this simplification is phenomenologically unacceptable. Often m_{τ} and m_{b} 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 multiflavor aliases
and handle processes such as
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 massdegenerate 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 processesThe process definition may contain an optional argument:
process processid = incomingparticles 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 reassignmentIt 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. Restrictions on matrix elementsAnother useful option is the setting
This option allows to select particular classes of Feynman graphs for the
process when using the O’Mega matrix element generator. The
The complete process e^{−}e^{+} → ννH, summed over all neutrino
generations, contains both ZH pair production (Higgsstrahlung) and
W^{+}W^{−}→ H fusion. The restrictions string selects the Higgsstrahlung
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
which is redundant here, however. The restriction keeps the full energy dependence in the intermediate propagator, so the BreitWigner 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 restrictions that can be combined with the restrictions above on intermediate propagators allow to exclude certain particles from intermediate propagators, or to exclude certain vertices from the matrix elements. For example,
would exclude all photon propagators from the matrix element and
leaves only the Z exchange here. In the same way,
excludes all Z and H propagators from the matrix elements. Besides excluding certain particles as internal lines, it is also possible to exclude certain vertices using the restriction command
This would generate the matrix element for the production of two W bosons at LEP without the nonAbelian vertex W^{+}W^{−}Z. Again, these restrictions are able to work on lists, so
would exclude all triple gauge boson vertices from the above process and leave only the tchannel neutrino exchange. It is also possible to exlude vertices by their coupling constants, e.g. the photon exchange in the process e^{+} e^{−} → µ^{+} µ^{−} can also be removed by the following restriction:
Here, qlep is the Fortran variable for the coupling constant of the electronpositronphoton vertex.
The Tab. 5.2 gives a list of options that can be applied to the O’Mega matrix elements. Other optionsThere are some further options that the O’Mega matrixelement generator can take. If desired, any string of options that is contained in this variable
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 tchannel. This is modified by the It is well known that for some processes, e.g., single W production from photonW fusion, gauge invariance puts constraints on the treatment of the unstableparticle 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 onshell. An alternative is
which puts zero width in the matrix element, so that gauge cancellations hold, and reinstates the schannel width in the appropriate places by an overall factor that multiplies the whole matrix element. Another possibility is
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 OpenMPOn 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 componentsIt was mentioned above that processes with flavor sums (in the initial or final state or both) have to be massdegenerate (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 phasespace parameterization and integration for the flavorsummed 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, gluinosquark associated production), or maybe leptoninclusive 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 socalled 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 oneloop program, OLP) can be mixed. The very same infrastructure can also be used for nexttoleading 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.:
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 phasespace 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 MonteCarlo 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:
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:
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 CompilationOnce processes have been set up, to make them available for integration they have to be compiled. More precisely, the matrixelement 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
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
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
code will be regenerated and recompiled even if WHIZARD would think that this is unncessary. The same effect is achieved by calling WHIZARD with a commandline switch,
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
Setting this switch is always a good idea when starting a new project, just in
case some old files clutter the working directory. When rerunning the same
script, possibly modified, the 5.4.6 Process librariesProcesses 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 Again, there are two possibilities. You may start the script with the specification
to open a library Alternatively, you may call WHIZARD with the option
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,
to compile only a specific subset. The command
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 Standalone WHIZARD with precompiled processesOnce you have set up a process library, it is straightforward to make a special standalone 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:
which produces an executable As an example, the script
will make a new executable program 5.5 BeamsBefore 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 setupIf 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 centerofmomentum energy of the collider by setting the value of √s, for instance
The beams statement comes into play if
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:
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.
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 singlebeam 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 anglesWHIZARD 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 protonproton collisions are equivalent: beams = p, p => pdf_builtin and beams = p, p => pdf_builtin 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:
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+" 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 nonvanishing 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 It is important that when a crossing angle is being specified, and the collision system consequently never is the centerofmomentum 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 5.5.3 LHAPDFFor incoming hadron beams, the beams statement specifies which structure functions are used. The simplest example is the study of partonparton scattering processes at a hadronhadron 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:
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 structurefunction 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:
for the actual LHAPDF 6 series, and
for LHAPDF5.Similarly, a member within the set is selected by the
numeric variable In some cases, different structure functions have to be chosen for the two beams. For instance, we may look at ep collisions:
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
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
integervalued parameter
Finally, the scattering of elementary photons on partons is described by
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 = "<pathtolhapdf>". Usually, it is not necessary to set this because WHIZARD detects this path via the lhapdfconfig 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 Builtin PDFsIn addition to the possibility of linking against LHAPDF, WHIZARD
comes with a couple of builtin PDFs which are selected via the
The default PDF set is CTEQ6L, but other choices are also available by
setting the string variable
would select the proton PDF from the MRST2004QED set. A list of all currently available PDFs can be found in Table 5.3.
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 matchingWhen the HOPPET tool [71] for hadroncollider 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 4flavor and 5flavor schemes in bparton 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 g→ bb 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β: pp → H with a partonic subprocess bb → H, or directly take the gluon PDFs and use pp → bb H with a partonic subprocess gg → b 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:
The first two processes are single top production from b PDFs, the last two processes contain an explicit g→ bb 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 functionsInitial 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 higherorder effect for the SLC and LEP I colliders. The softcollinear and soft photon radiation can indeed be resummed/exponentiated to all orders in perturbation theory [7], while higher orders in hardcollinear 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 softcollinear photons as well as up to the third order in hardcollinear photons. It can be switched on by the following statement:
As the ISR structure function is a singlebeam structure function, this expression is synonymous for
The ISR structure function can again be applied to only one of the two beams, e.g. in a HERAlike setup:
Their are several options for the leptoncollider ISR structure function that are summarized in the following:
The maximal order of the hardcollinear 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. 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. While the ISR structure function is evaluated in the collinear limit, it is possible to generate transverse momentum for both the radiated photons and the recoiling partonic system. We recommend to stick to the collinear approximation for the integration step. Integration cuts should be set up such that they do not significantly depend on photon transverse momentum. In a subsequent simulation step, it is possible to transform the events with collinear ISR radiation into more realistic events with noncollinear radiation. To this end, WHIZARD provides a separate ISR photon handler which can be activated in the simulation step. The algorithm operates on the partonic event: it takes the radiated photons and the partons entering the hard process, and applies a p_{T} distribution to those particles and their interaction products, i.e., all outgoing particles. Cuts that depend on photon p_{T} may be applied to the modified events. For details on the ISR photon handler, cf. Sec. 10.4. The flag ?isr_recoil switches on p_{T} recoil of the emitting lepton against photon radiation during integration; per default it is off. The flag ?isr_keep_energy controls the mode of onshell projection for the splitting process with p_{T}. Note that this feature is kept for backwards compatibility, but should not be used for new simulations. The reason is as follows: For a fraction of events, p_{T} will become significant, and (i) energy/momentum nonconservation, applied to both beams separately, can lead to unexpected and unphysical effects, and (ii) the modified momenta enter the hard process, so the collinear approximation used in the ISR structure function computation does not hold. 5.5.7 Lepton Collider BeamstrahlungAt 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 beambeam interactions in leptoncollider physics is GuineaPig++ [10, 11, 12]. A direct interface between this tool GuineaPig++ and WHIZARD had existed as an inofficial addon to the legacy branch WHIZARD1, but is no longer applicable in WHIZARD2. A WHIZARDinternal 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:
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):
or
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
This is an overview over all options and flags for the CIRCE1 setup for lepton collider beamstrahlung:
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 semianalytic approximation formulae for the spectra, or a MonteCarlo 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 phasespace 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 lowenergy 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:
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 eventsAs 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 pregenerated 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:
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/beamsim. 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 beamenergy spreadReal 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.
(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 approximationThe equivalent photon approximation (EPA) uses an onshell approximation for the e → eγ collinear splitting to allow the simulation of photoninduced backgrounds in lepton collider physics. The original concept is that of the WeizsäckerWilliams approximation [21, 22, 23]. This is a singlebeam structure function that can be applied to both beams, or also to one beam only. Examples are:
or for a single beam:
The last process allows the reaction of (quasi) onshell photons with protons. In the following, we collect the parameters and flags that can be adjusted when using the EPA inside WHIZARD:
The adjustable parameters are partially similar to the parameters in the QED initialstate 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 momentumtransfer cutoff for the EPA structure function, epa_q_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 softphoton 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:
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^{+} e^{−} → e^{+} 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 tradeoff 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 approximationAn 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 highenergy weak vectorboson fusion and scattering processes at hadron colliders, particularly the Superconducting SuperCollider (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 pp → jjVV and subsequent decays of the bosons V ≡ W, 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 singlebeam 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:
And this is for LHC or a higherenergy followup collider (which also shows the concatenation of the singlebeam structure functions, applied to both beams consecutively, cf. Sec. 5.5.14:
Again, we list all the options, parameters and flags that can be adapted for the EWA:
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 f → f 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 p_{T} 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 noncollinear kinematics violates 4four momentum conservation, so there are two choices: either to conserve the energy (?ewa_keep_energy = true) or to conserve 3momentum (?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 functionsIn 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 falloff with the Z resonance in e^{+}e^{−} → µ^{+} µ^{−}, where the syntax is very easy:
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 finalstate 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 spectraOne option that has been discussed as an alternative possibility for a highenergy 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 photoninitiated process:
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:
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 functionsAs 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 singlebeam structure functions or pair spectra. One important thing is whether there is a phasespace mapping for these structure functions. Also, there are some combinations which do not make sense from the physics point of view, for example using leptoncollider 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 beambeam interactions) and electron QED initialstate radiation are both switched on:
Another possibility is the simulation of photoninduced backgrounds at ILC or CLIC, using beamstrahlung and equivalent photon approximation (EPA):
or with beam events from a data file:
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 highenergy vector boson scattering:
Note that this last case involves a flavor sum over the five active quark (and antiquark) 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:
5.6 Polarization5.6.1 Initial state polarizationWHIZARD 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 statement^{5}:
The command beams_pol_fraction gives the degree of polarization of the two beams:
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 % lefthanded electron polarization means that 80 % of the electron beam are polarized, 20 % are unpolarized, i.e. 20 % have half left and half righthanded polarization each. Hence, 90 % of the electron beam is lefthanded, 10 % is righthanded. 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%) lefthanded 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 nonzero, but might not be shown. E.g. ILClike, 80 % polarization of the electrons, 30 % positron polarization will be specified like this for lefthanded electrons and righthanded positrons:
The screen output will be like this:
But because of the fraction of unpolarized electrons and positrons, the spin density matrices for electrons and positrons are:
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):
Offdiagonal entries that are equal to one (up to the normalization) of the spindensity 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 offdiagonal 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.
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:
Although beam polarization should be straightforward to use, some pitfalls exist for the unwary:
5.6.2 Final state polarizationFinal 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
The effect of polarized can be reversed with the unpolarized statement which has the same syntax. For example,
will cause the polarization of all final state W and Z bosons to be retained, while
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.
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
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
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 sectionsIntegrating 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 multichannel MonteCarlo integration. The integration, in turn, requires a phasespace setup, i.e., a collection of phasespace channels, which are mappings of the unit hypercube onto the complete space of multiparticle kinematics. This phasespace information is encoded in the file xxx.phs, where xxx is the process tag. WHIZARD generates the phasespace file on the fly and can reuse it in later integrations. For each phasespace 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 phasespace points, so they adapt to the actual phasespace 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 IntegrationSince everything can be handled automatically using default parameters, it often suffices to write the command
for integrating the process with name tag proc1, and similarly
for integrating several processes consecutively. Options to the integrate command are specified, if not globally, by a local option string
(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 matrixelement 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^{+}e^{−}→ qq 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 p_{T} and energy of jets, and on the invariant mass of jet pairs. Here is the script:
After the run is finished, the integration output looks like
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 basis. The next two columns display the error in percent, and the accuracy which is the same error normalized by √n_{calls}. The accuracy value has the property that it is independent of n_{calls}, it describes the quality of adaptation of the current grids. Goodquality 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 phasespace for this process and set of cuts is wellbehaved. 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 suboptimally 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 randomnumber generator will differ by default. It is possible to get exactly reproducible results by setting the randomnumber seed explicitly, e.g.,
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 L^{A}T_{E}X 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 IDsA 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 rerun: all results that have been overwritten have to be recreated. To avoid this, the user may identify a specific run by a stringvalued ID, e.g.
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:
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 iterationsWHIZARD has some predefined numbers of iterations and calls for the first and second integration pass, respectively, which depend on the number of initial and finalstate particles. They are guesses for values that yield goodquality 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 shortcut 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:
This is a commaseparated 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
where the relative channel weights will not be adjusted (because this is the final pass). This is appropriate for wellbehaved 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 asis, 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:
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
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 spaceBefore integrate can start its work, it must have a phasespace configuration for the process at hand. The method for the phasespace 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 finetuned for electroweak processes and might not be the ideal for for pure jet cross sections. In future versions of WHIZARD, more options for phasespace parameterizations will be made available, e.g. the RAMBO algorithm and its massive cousin, and phasespace parameterizations that take care of the dipolelike emission structure in collinear QCD (or QED) splittings. For the standard method, the phasespace parameterization is laid out in an ASCII file <processname>_i<comp>.phs. Here, <processname> 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 pp → Z + nj and pp → W + 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 phasespace file and inspect and/or
modify it before proceeding further. To this end, there is the parameter
Things might go wrong with the default phasespace generation, or manual intervention might be necessary to improve later performance. There are a few parameters that control the algorithm of phasespace generation. To understand their meaning, you should realize that phasespace parameterizations are modeled after (dominant) Feynman graphs for the current process. The main phase space setup woodFor the main phasespace parameterization of WHIZARD, which is called "wood", there are many different parameters and flags that allow to tune and customize the phasespace setup for every certain process: The parameter There is a similar parameter There are two numerical parameters that describe whether particles are treated
like massless particles in particular situations. The value of
Such logarithmic mappings need a dimensionful scale as parameter. There are
three such scales, all with default value 10 GeV: There are more flags. These and more details about the phase space parameterization will be described in Sec. 8.2. 5.7.5 CutsWHIZARD 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
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 righthand 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 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 righthand 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
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:
Example:
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
processindependent 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
5.7.6 QCD scale and couplingWHIZARD treats all physical parameters of a model, the coefficients in the Lagrangian, as constants. As a leadingorder 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
This is controlled by the parameter One option is to set In the very same way, the α_{s} running from the PDFs implemented
intrinsically in WHIZARD can be taken by setting
If this is not appropriate, there are again two possibilities. If
Otherwise there is the option to set Note that for using one of the running options for α_{s}, always ?alphas_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 userdefined scale, the special object scale. This is assigned and used just like the cuts object. The righthand side is a realvalued expression. Here is an example:
This selects the p_{T} value of the first entry in the list of colored particles sorted by decreasing p_{T}, i.e., the p_{T} 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 factorIt is possible to reweight the integrand by a userdefined function of the event kinematics. This is done by specifying a weight expression. Syntax and usage is exactly analogous to the scale expression. Example:
We should note that the phasespace 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 EventsAfter the cross section integral of a scattering process is known (or the partialwidth integral of a decay process), WHIZARD can generate event samples. There are two limiting cases or modes of event generation:
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 MonteCarlo integration yields nonuniform 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 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 5.8.1 SimulationThe simulate command generates an event sample. The number of events
can be set either by specifying the integer variable As usual, both parameters can be set either as global or as local parameters:
In the second example, both 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
will produce an event file This event file is in a machinedependent binary format, so it is not of
immediate use. Its principal purpose is to serve as a cache: if you rerun
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
timeconsuming 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
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., Generating an event sample can serve several purposes. First of all, it can be analyzed directly, by WHIZARD’s builtin 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 eventfile 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.
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 DecaysNormally, 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 pp→ W^{+}W^{−}, the finalstate particles
are onshell 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., pp→ udūd. In
this case, the intermediate vector bosons, if reconstructed, are offshell 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
nonresonant intermediate states. (This can be restricted via the
Another approach is to factorize the process in production (of W bosons) and decays (W→ qq). This is actually the traditional approach, since it is much less computingintensive. The factorization neglects all offshell effects and irreducible background diagrams that do not have the decaying particles as an intermediate resonance. While WHIZARD is able to deal with multiparticle 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 pp → Wj final state:
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 narrowwidth 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 onshell, or otherwise gaugeinvariance issues could spoil the calculation. (There are plans for a future WHIZARD version to also include BreitWigner or Gaussian distributions when using the factorized approach.) Decays can be concatenated, e.g. for top pair production and decay, e^{+} e^{−} → t t with decay t → W^{+} b, and subsequent leptonic decay of the W as in W^{+} → µ^{+} ν_{µ}:
Note that in this case the final state in the event file will consist of t b µ^{+} ν_{µ} because the antitop is not decayed. If more than one decay process is being specified like in
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
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
constructor for the particle <particletag>. Resetting branching fractionsAs described above, decay processes that appear in a simulation must
first be integrated by the program, either explicitly via the
However, it may be desired to rather use different branchingfraction
values for the decays of a particle, for instance, NLOcorrected
values. In fact, after WHIZARD has integrated any process, the
integration result becomes available as an ordinary
SINDARIN variable. For instance, if a decay process has the ID
The new value will be used for all subsequent Higgs branchingratio calculations and decays, if an unstable Higgs appears in a process for simulation. Spin correlations in decaysBy 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 productiondecay 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:
Here, the isotropic option completely eliminates spin correlations. The diagonaldecays option eliminates just the offdiagonal entries of the W spindensity 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 definitehelicity 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 decaysA 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):
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:
Future shorter notation for decaysIn 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:
as well as
5.8.3 Event formatsAs mentioned above, the internal WHIZARD event format is a machinedependent event format. There are a series of humanreadable 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 VisualizationSINDARIN natively supports basic methods of data analysis and visualization which are frequently used in highenergy physics studies. Data generated during script execution, in particular simulated event samples, can be analyzed to evaluate further observables, fill histograms, and draw twodimensional 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 ObservablesAnalyses in highenergy 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 realvalued variable with a statistical distribution. It is declared by a command of the form observable analysistag where analysistag 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 analysistag (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 (analysistag) Here is an example:
The result is displayed on screen:
5.9.2 The analysis expressionThe most common application is the computation of event observables – for instance, a forwardbackward 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 = logicalexpr 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 forwardbackward symmetry in a process ee_mumu with final state µ^{+}µ^{−}:
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:
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:
Alternatively, we can specify a separate selection expression:
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:
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:
5.9.3 HistogramsIn SINDARIN, a histogram is declared by the command histogram analysistag (lowerbound, upperbound) This creates a histogram data structure for an (unspecified) observable. The entries are organized in bins between the real values lowerbound and upperbound. 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 (realexpr) puts the value of realexpr 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
blankseparated 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 5.9.4 PlotsWhile 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 analysistag Like observables and histograms, plots are filled by the record command or expression. To this end, it can take two arguments, record (xexpr, yexpr) or up to four: record (xexpr, yexpr, yerror) 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,y,Δ y, Δ x. 5.9.5 Analysis OutputThere 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.
5.10 Custom Input/OutputWHIZARD 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 processspecific 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 Filesopen_out
Open an external file for writing. If the file exists, it is overwritten without warning, otherwise it is created. Example:
close_out
Close an external file that is open for writing. Example:
5.10.2 Printing Dataprintf
Format ⟨dataobjects⟩ according to ⟨formatstringexpr⟩ 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 ⟨formatstringexpr⟩ 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:
For more details, consult the 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 Note that for simply printing a text string, you may call printf with just a format string and no data arguments. Examples:
There is a related sprintf function, cf. Sec. 5.1.5. 5.11 WHIZARD at nexttoleading order5.11.1 PrerequisitesA full NLO computation requires virtual matrix elements obtained from loop diagrams. Since O’Mega cannot calculate such diagrams, external programs are used. WHIZARD has a generic interface to matrixelement generators that are BLHAcompatible. Explicit implementations exist for Gosam, OpenLoops and Recola. Setting up GosamThe installation of Gosam is detailed on the HepForge page https://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 relevant 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 oneloop 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. A more convenient way to install Gosam, is the automatic installation script https://gosam.hepforge.org/gosam_installer.py. Setting up OpenLoopsThe installation of OpenLoops is explained in detail on the HepForge page https://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 withopenloops=<openloops_dir>. OpenLoops can be checked out with
The program can be build by running scons or ./scons, a local version that is included in the OpenLoops directory. This produces the script ./openloops, which is the main hook for the further usage of the program. OpenLoops works by downloading prebuild process libraries, which have to be installed for each individual process. This requires the file openloops.cfg, which should contain the following content:
The first line instructs OpenLoops to also look for process libraries in an additional lepton collider repository. The second line triggers the inclusion of N+1particle treelevel 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
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. The process library ppll, for example, can also be used to compute the matrix elements for e^{+} e^{−} → q q (massless quarks only). The massive case of the top quark is handled in eett. 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
where $OPENLOOPS_PATH is the directory the OpenLoops executable is located in. OpenLoops oneloop diagrams can then be used with the SINDARIN option
The functional tests which check the OpenLoops functionality 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. OpenLoops SINDARIN flagsSeveral SINDARIN options exist to control the behavior of OpenLoops.
Setting up RecolaThe installation of Recola is explained in detail on the HepForge page https://recola.hepforge.org. In the following the main steps for usage with WHIZARD are summarized. The minimal required version number of Recola is 1.2. Recola can be linked to WHIZARD during configuration by including
where $RECOLA_PATH is the directory the Recola library is located in. Recola can then be used with the SINDARIN option
or any other of the matrix element methods. Note that there might be a clash of the Collier libraries when you have Collier installed both via Recola and via OpenLoops, but have compiled them with different Fortran compilers. 5.11.2 NLO cross sectionsAn NLO computation can be switched on in SINDARIN with
where the nlo_calculation can be followed by a list of strings specifying the desired NLOcomponents to be integrated, i.e. born, real, virtual, dglap, (for hadron collisions) or mismatch (for the soft mismatch in resonanceaware computations) and full. The full option switches on all components and is required if the total NLO result is desired. For example, specifying
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 and a fine tuned control over the iterations in each component. In the combined integration mode, all components are added up during integration so that the sum of them is evaluated. Here, only one integration 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 is in this case achieved activating the resonanceaware FKS subtraction,
This mode comes with an additional integration component, the socalled soft mismatch. Note that you can modify the number of iterations in each component with the multipliers:
5.11.3 Fixedorder NLO eventsFixedorder NLO events can also be produced in three different modes: Combined weighted, combined unweighted and separated weighted.
Weighted fixedorder NLO events are supported by any output format that supports weights like the HepMC format and unweighted NLO events work with any format. The output can either be written to disk or put into a FIFO to interface it to an analysis program without writing events to file. The weights in the real event output, both in the combined and separate weighted mode, are divided by a factor N_{phs} + 1. This is to account for the fact that we artificially increase the number of events in the output file. Thus, the sum of all event weights correctly reproduces the total cross section. 5.11.4 POWHEG matchingTo match the NLO computation with a parton shower, WHIZARD supports the POWHEG matching. It generates a distribution according to The subscript s refers to the singular part of the real component, cf. to the next subsection. Eq. (1) produces either no or one additional emission. These events can then either be analyzed directly or passed on to the parton shower^{6} for the full simulation. You activate this with
The p_{T}^{min} of Eq. (1) can be set with powheg_pt_min. It sets the minimal scale for the POWHEG evolution and should be of order 1 GeV and set accordingly in the interfaced shower. The maximal scale is currently given by sqrts but should in the future be changeable with powheg_pt_min. Note that the POWHEG event generation needs an additional grid for efficient event generation that is automatically generated during integration. Further options that steer the efficiency of this grid are powheg_grid_size_xi and powheg_grid_size_y. 5.11.5 Separation of finite and singular contributionsFor both the pure NLO computations as well as the POWHEG event generation, WHIZARD supports the partitioning of the real into finite and singular contributions with the flag
The finite contributions, which by definition should not contain soft or collinear emissions, will then integrate like a ordinary LO integration with one additional particle. Similarly, the event generation will produce only real events without subtraction terms with Born kinematics for this additional finite component. The POWHEG event generation will also only use the singular parts. The current implementation uses the following parametrization
Thus, a point is singular (F=1), if any of the FKS tuples forms an invariant mass that is smaller than the hardness scale h. This parameter is controlled in SINDARIN with real_partition_scale. This simplifies in massless case to
