Chapter 5 SINDARIN in Details
5.1 Data and expressions5.1.1 Real-valued objectsReal literals have their usual form, mantissa and, optionally, exponent:
0. 3.14 -.5 2.345e-3 .890E-023
Internally, real values are treated as double precision. The values are read by the Fortran library, so details depend on its implementation. A special feature of SINDARIN is that numerics (real and integer) can be immediately followed by a physical unit. The supported units are presently hard-coded, they are
meV eV keV MeV GeV TeV
nbarn pbarn fbarn abarn rad mrad degree % If a number is followed by a unit, it is automatically normalized to the corresponding default unit: 14.TeV is transformed into the real number 14000. Default units are GeV, fbarn, and rad. The % sign after a number has the effect that the number is multiplied by 0.01. Note that no checks for consistency of units are done, so you can add 1 meV + 3 abarn if you absolutely wish to. Omitting units is always allowed, in that case, the default unit is assumed. Units are not treated as variables. In particular, you can’t write theta / degree, the correct form is theta / 1 degree. There is a single predefined real constant, namely π which is referred to by the keyword pi. In addition, there is a single predefined complex constant, which is the complex unit i, being referred to by the keyword I. The arithmetic operators are + - * / ^
with their obvious meaning and the usual precedence rules. SINDARIN supports a bunch of standard numerical functions, mostly equivalent to their Fortran counterparts:
abs conjg sgn mod modulo
sqrt exp log log10 sin cos tan asin acos atan sinh cosh tanh (Unlike Fortran, the sgn function takes only one argument and returns 1., or −1.) The function argument is enclosed in brackets: sqrt (2.), tan (11.5 degree). There are two functions with two real arguments:
max min
Example: The following functions of a real convert to integer:
int nint floor ceiling and this converts to complex type:
complex
Real values can be compared by the following operators, the result is a logical value: == <>
> < >= <=
In SINDARIN, it is possible to have more than two operands in a logical expressions. The comparisons are done from left to right. Hence, 115 GeV < mH < 180 GeV
is valid SINDARIN code and evaluates to true if the Higgs mass is in the given range. Tests for equality and inequality with machine-precision real numbers are
notoriously unreliable and should be avoided altogether. To deal with this
problem, SINDARIN has the possibility to make the comparison operators
“fuzzy” which should be read as “equal (unequal) up to an absolute
tolerance”, where the tolerance is given by the real-valued intrinsic
variable tolerance. This variable is initially zero, but can be
set to any value (for instance, tolerance = 1.e-13 by the user.
Note that for non-zero tolerance, operators like
5.1.2 Integer-valued objectsInteger literals are obvious:
1 -98765 0123
Integers are always signed. Their range is the default-integer range as determined by the Fortran compiler. Like real values, integer values can be followed by a physical unit: 1 TeV, 30 degree. This actually transforms the integer into a real. Standard arithmetics is supported: + - * / ^
It is important to note that there is no fraction datatype, and pure integer arithmetics does not convert to real. Hence 3/4 evaluates to 0, but 3 GeV / 4 GeV evaluates to 0.75. Since all arithmetics is handled by the underlying Fortran library, integer overflow is not detected. If in doubt, do real arithmetics. Integer functions are more restricted than real functions. We support the following:
abs sgn mod modulo
max min and the conversion functions
real complex
Comparisons of integers among themselves and with reals are possible using the same set of comparison operators as for real values. This includes the operators with a finite tolerance. 5.1.3 Complex-valued 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 Logical-valued 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 String-valued 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 format-string
(arg-list)
This is an expression that evaluates to a string. The format string contains the mentioned conversion specifiers. The argument list is optional. The arguments are separated by commas. Allowed arguments are integer, real, logical, and string variables, and numeric expressions. Logical and string expressions can also be printed, but they have to be dressed as anonymous variables. A logical anonymous variable has the form ?(logical_expr) (example: ?(mH > 115 GeV)). A string anonymous variable has the form $(string-expr). Example:
The related printf command with the same syntax prints the formatted string to standard output2. 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 multi-species 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 (single-particle) 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 four-momenta), all u, d, or s quarks, and all colored particles, respectively. A subevent can contain pseudoparticles, i.e., particle combinations. That is, the four-momenta of distinct particles are combined (added conmponent-wise), and the results become subevent elements just like ordinary particles. The (pseudo)particles in a subevent are non-overlapping. That is, for any of the particles in the original event, there is at most one (pseudo)particle in the subevent in which it is contained. Sometimes, variables (actually, named constants) of type subevent are useful.
Subevent variables are declared by the subevt keyword, and their
names carry the prefix
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 four-momentum. The particles argument may either be a subevt expression or an alias expression. The result is a one-entry subevt. In the second form, only those particles are collected which satisfy the condition, a logical expression. Example: collect if Pt > 10 GeV [colored] The third version is useful if you want to put binary observables (i.e., observables constructed from two different particles) in the condition. The ref_particles provide the second argument for binary observables in the condition. A particle is taken into account if the condition is true with respect to all reference particles that do not overlap with this particle. Example: collect if Theta > 5 degree [photon, charged]: combine all photons that are separated by 5 degrees from all charged particles. clustercluster [particles] First version: collect all particle momenta in the argument and cluster them to a set of jets. The particles argument may either be a subevt expression or an alias expression. The result is a one-entry subevt. In the second form, only those particles are clustered which satisfy the condition, a logical expression. Example: cluster if Pt > 10 GeV [colored] This command is available from WHIZARD version 2.2.1 on, and only if the FastJet package has been installed and linked with WHIZARD (cf. Sec.2.2.12); in a future version of WHIZARD it is foreseen to have also an intrinsic clustering package inside WHIZARD which will be able to support some of the clustering algorithms below. To use it in an analysis, you have to set the variable jet_algorithm to one of the predefined jet-algorithm values (integer constants): kt_algorithm and the variable jet_r to the desired R parameter value, as appropriate for the analysis and the jet algorithm. Example:
select_b_jet, select_non_b_jet, select_c_jet, select_light_jetThis command is available from WHIZARD version 2.8.1 on, and it only generates anything non-trivial if the FastJet package has been installed and linked with WHIZARD (cf. Sec.2.2.12). It only returns sensible results when it is applied to subevents after the cluster command (cf. the paragraph before). It is similar to the select command, and accepts a logical expression as a possible condition. The four commands select_b_jet, select_non_b_jet, select_c_jet, and select_light_jet select b jets, non-b jets (anything lighter than bs), c jets (neither b nor light) and light jets (anything besides b and c), respectively. An example looks like this:
photon_isolationThis command is available from WHIZARD version 2.8.1 on. It provides isolation of photons from hadronic (and possibly electromagnetic) activity in the event to define a (especially) NLO cross section that is completely perturbative. The isolation criterion according to Frixione, cf. [61], removes the non-perturbative contribution from the photon fragmentation function. This command can in principle be applied to elementary hard process partons (and leptons), but generates something sensible only if the FastJet package has been installed and linked with WHIZARD (cf. Sec.2.2.12). There are three parameters which allow to tune the isolation, photon_iso_r0, which is the radius Rγ0 of the isolation cone, photon_iso_eps, which is the fraction єγ of the photon (transverse) energy that enters the isolation criterion, and the exponent of the isolation cone, photon_iso_n, nγ. For more information cf. [61]. The command allows also a conditional cut on the photon which is applied before the isolation takes place. The first argument are the photons in the event, the second the particles from which they should be isolated. If also the electromagnetic activity is to be isolated, photons need to be isolated from themselves and must be included in the second argument. This is mandatory if leptons appear in the second argument. Two examples look like this:
photon_recombinationphoton_recombination [particles] This function, which maps a subevent into another subevent, is used for electroweak (and mixed coupling) higher order calculations. It takes the selection of photons in particles (for the moment, WHIZARD restricts this to one explicit photon in the final state) and recombines it with the closest non-photon particle from particles in R-distance, if the R-distance is smaller than the parameter set by photon_rec_r0. Otherwise the particles subevent is left unchanged so that it may contain possibly non-recombined photons. The logical variable ?keep_flavors_when_recombining determines whether WHIZARD keeps the flavor of the particle with which the photon is recombined into the pseudoparticle, the default being true. An example for photon recombination is shown here:
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 single-particle subevent. In the first version, it contains the first particle in the subevent particles. In the second version, the particle with index index-value is returned, where index-value is an integer expression. If its value is negative, the index is counted from the end of the subevent. The order of particles in an event or subevent is not always well-defined, so you may wish to sort the subevent before applying the extract function to it. 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 ref-particle, 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 (unary observables), or two particle momenta (binary observables) or all momenta of the particles (n-ary/subeventary observables) 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 electron-neutrino pair in the event. The last example illustrates why observables are treated like variables, even though they are functions of particles: the eval construct with the particle reference in square brackets after the expression allows to compute derived observables – observables which are functions of new observables – without the need for hard-coding them as new functions. For subeventary observables, e.g. Ht, the momenta of all particles in the subevent are taken to evaluate the observables, e.g. eval Ht/2 [t:T:Z:jet] takes the (half of) the transverse mass of all tops, Zs and jets in the final state. sumThis SINDARIN statement works similar to the eval statement above, with the syntax sum <expr> [<subevt>] It sums the <expr> over all elements of the subevents <subevt>, e.g. sum sqrt(Pt^2 + M^2)/2 [t:T:H:Z] would calculate the transverse mass (square root of the sum of squared transverse momentum and squared mass) of all tops, Higgs and Z bosons in the final state. prodIdentical to sum, but takes the product, not the sum of the expression <expr> evaluated over the full subevent. Syntax: prod <expr> [<subevt>] 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 two-argument case it must be true for all non-overlapping combinations of particles in the two subevents. If one of the arguments is the empty subevent, the result is also true. anyany logical_expr [particles] The any construct is true if the logical expression is true for at least one particle or non-overlapping 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 non-overlapping 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 non-overlapping particle pairs (that pass the test, if any). Predefined observablesThe following real-valued observables are available in SINDARIN for use in eval, all, any, no, and count constructs. The argument is always the subevent or alias enclosed in square brackets.
There are also integer-valued observables:
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 phase-space setup; the vertices used for generating amplitudes and the corresponding Feynman rules are stored in different files within the O’Mega source tree. In a SINDARIN script, a model is a special object of type model. There is always a current model. Initially, this is the SM, so on startup WHIZARD reads the SM.mdl model file and assigns its content to the current model object. (You can change the default model by the --model option on the command line. Also the preloading of a model can be switched off with the --no-model option) Once the model has been loaded, you can define processes for the model, and you have all independent model parameters at your disposal. As noted before, these are intrinsic parameters which need not be declared when you assign them a value, for instance:
Other parameters are derived. They can be used in expressions like any other parameter, they are also intrinsic, but they cannot be modified directly at all. For instance, the electromagnetic coupling ee is a derived parameter. If you change either GF (the Fermi constant), mW (the W mass), or mZ (the Z mass), this parameter will reflect the change, but setting it directly is an error. In other words, the SM is defined within WHIZARD in the GF-mW-mZ scheme. (While this scheme is unusual for loop calculations, it is natural for a tree-level event generator where the Z and W poles have to be at their experimentally determined location3.) The model also defines the particle names and aliases that you can use for defining processes, cuts, or analyses. If you would like to generate a SUSY process instead, for instance, you can assign a different model (cf. Table 10.1) to the current model object:
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 mb 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 mb value has been remembered. If you want both mass values to agree, you have to set them separately in the context of their respective model. Although this might seem cumbersome at first, it is nevertheless a sensible procedure since the parameters defined by the user might anyhow not be defined or available for all chosen models. When using two different models which need an SLHA input file, these have to be provided for both models. Within a given scope, there is only one current model. The current model can be reset permanently as above. It can also be temporarily be reset in a local scope, i.e., the option body of a command or the body of a scan loop. It is thus possible to use several models within the same script. For instance, you may define a SUSY signal process and a pure-SM background process. Each process depends only on the respective model’s parameter set, and a change to a parameter in one of the models affects only the corresponding process. 5.4 ProcessesThe purpose of WHIZARD is the integration and simulation of high-energy physics processes: scatterings and decays. Hence, process objects play the central role in SINDARIN scripts. A SINDARIN script may contain an arbitrary number of process definitions. The initial states need not agree, and the processes may belong to different physics models. 5.4.1 Process definitionA process object is defined in a straightforward notation. The definition syntax is straightforward:
process process-id = incoming-particles Here are typical examples:
Throughout the program, the process will be identified by its process-id, so this is the name of the process object. This identifier is arbitrary, chosen by the user. It follows the rules for variable names, so it consists of alphanumeric characters and underscores, where the first character is not numeric. As a special rule, it must not contain upper-case characters. The reason is that this name is used for identifying the process not just within the script, but also within the Fortran code that the matrix-element generator produces for this process. After the equals sign, there follow the lists of incoming and outgoing particles. The number of incoming particles is either one or two: scattering processes and decay processes. The number of outgoing particles should be two or larger (as 2→ 1 processes are proportional to a δ function they can only be sensibly integrated when using a structure function like a hadron collider PDF or a beamstrahlung spectrum.). There is no hard upper limit; the complexity of processes that WHIZARD can handle depends only on the practical computing limitations (CPU time and memory). Roughly speaking, one can assume that processes up to 2→ 6 particles are safe, 2→ 8 processes are feasible given sufficient time for reaching a stable integration, while more complicated processes are largely unexplored. We emphasize that in the default setup, the matrix element of a physics process is computed exactly in leading-order perturbation theory, i.e., at tree level. There is no restriction of intermediate states, the result always contains the complete set of Feynman graphs that connect the initial with the final state. If the result would actually be expanded in Feynman graphs (which is not done by the O’Mega matrix element generator that WHIZARD uses), the number of graphs can easily reach several thousands, depending on the complexity of the process and on the physics model. More details about the different methods for quantum field-theoretical matrix elements can be found in Chap. 9. In the following, we will discuss particle names, options for processes like restrictions on intermediate states, parallelization, flavor sums and process components for inclusive event samples (process containers). 5.4.2 Particle 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 mb can also be neglected, but this excludes processes where the Higgs couplings of τ or b are relevant. Setting fermion masses to zero enables, furthermore, the possibility to define multi-flavor aliases
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 mass-degenerate particles). If you want to generate inclusive processes with sums over particles of different masses (e.g. summing over W/Z in the final state etc.), confer below the section about process components, Sec. 5.4.4. Assignments of masses, widths and other parameters are actually in effect when a process is integrated, not when it is defined. So, these assignments may come before or after the process definition, with no significant difference. However, since flavor summation requires masses to be zero, the assignments may be put before the alias definition which is used in the process. The muon, tau, and the heavier quarks are actually unstable. However, the width is set to zero because their decay is a macroscopic effect and, except for the muon, affected by hadron physics, so it is not described by WHIZARD. (In the current WHIZARD setup, all decays occur at the production vertex. A future version may describe hadronic physics and/or macroscopic particle propagation, and this restriction may be eventually removed.) 5.4.3 Options for processesThe process definition may contain an optional argument:
process process-id = incoming-particles 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 (Higgs-strahlung) and
W+W−→ H fusion. The restrictions string selects the Higgs-strahlung
graph where the initial electrons combine to a Z boson. Here, the particles
in the process are consecutively numbered, starting with the initial
particles. An alternative for the same selection would be
which is redundant here, however. The restriction keeps the full energy dependence in the intermediate propagator, so the Breit-Wigner shape can be observed in distributions. This breaks gauge invariance, in particular if the intermediate state is off shell, so you should use the feature only if you know the implications. For more details, cf. the Chap. 9 and the O’Mega manual. Other 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 non-Abelian 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 t-channel 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 electron-positron-photon 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 matrix-element 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 t-channel. This is modified by the It is well known that for some processes, e.g., single W production from photon-W fusion, gauge invariance puts constraints on the treatment of the unstable-particle width. By default, O’Mega puts a nonzero width in the s channel only. This correctly represents the resummed Dyson series for the propagator, but it violates QED gauge invariance, although the effect is only visible if the cuts permit the photon to be almost on-shell. An alternative is
which puts zero width in the matrix element, so that gauge cancellations hold, and reinstates the s-channel width in the appropriate places by an overall factor that multiplies the whole matrix element. Note that the fudged width option only applies to charged unstable particles, such as the W boson or top quark. Another possibility is
which puts the width both in the s- and in the t-channel like diagrams. A third option is provided by the running width scheme
which applies the width only for s-channel like diagrams and multiplies it by a factor of p2 / M2. The additional p2-dependent factor mimicks the momentum dependence of the imaginary part of a vacuum polarization for a particle decaying into massles decay products. It is noted that none of the above options preserves gauge invariance. For a gauge preserving approach (at least at tree level), O’Mega provides the complex-mass scheme
However, in this case, one also has to modify the model in usage. For example, the parameter setting for the Standard Model can be changed by,
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 mass-degenerate (in most cases massless) in all particles that are summed over at a certain position. This condition is necessary in order to use the same phase-space parameterization and integration for the flavor-summed process. However, in many applications the user wants to handle inclusive process definitions, e.g. by defining inclusive decays, inclusive SUSY samples at hadron colliders (gluino pairs, squark pairs, gluino-squark associated production), or maybe lepton-inclusive samples where the tau and muon mass should be kept at different values. In WHIZARD from version v2.2.0 on, there is the possibility to define such inclusive process containers. The infrastructure for this feature is realized via so-called process components: processes are allowed to contain several process components. Those components need not be provided by the same matrix element generator, e.g. internal matrix elements, O’Mega matrix elements, external matrix element (e.g. from a one-loop program, OLP) can be mixed. The very same infrastructure can also be used for next-to-leading order (NLO) calculations, containing the born with real emission, possible subtraction terms to make the several components infrared- and collinear finite, as well as the virtual corrections. Here, we want to discuss the use for inclusive particle samples. There are several options, the simplest of which to add up different final states by just using the + operator in SINDARIN, e.g.:
The brackets are not only used for a better grouping of the expressions, they are not mandatory for WHIZARD to interpret the sum correctly. When integrating, WHIZARD tells you that this a process with three different components: | Initializing integration for process multi_comp_1_p1:
| ------------------------------------------------------------------------
| Process [scattering]: 'multi_comp'
| Library name = 'default_lib'
| Process index = 1
| Process components:
| 1: 'multi_comp_i1': e-, e+ => m-, m+ [omega]
| 2: 'multi_comp_i2': e-, e+ => t-, t+ [omega]
| 3: 'multi_comp_i3': e-, e+ => A, A [omega]
| ------------------------------------------------------------------------
A different phase-space setup is used for each different component. The integration for each different component is performed separately, and displayed on screen. At the end, a sum of all components is shown. All files that depend on the components are being attached an _i<n> where <n> is the number of the process component that appears in the list above: the Fortran code for the matrix element, the .phs file for the phase space parameterization, and the grid files for the VAMP Monte-Carlo integration (or any other integration method). However, there will be only one event file for the inclusive process, into which a mixture of events according to the size of the individual process component cross section enter. More options are to specify additive lists of particles. WHIZARD then expands the final states according to tensor product algebra:
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 matrix-element generator O’Mega (and it works similarly if a different matrix element method is chosen) is called to generate matrix element code, the compiler is called to transform this Fortran code into object files, and the linker is called to collect this in a dynamically loadable library. Finally, this library is linked to the program. From version v2.2.0 of WHIZARD this is no longer done by system calls of the OS but steered via process library Makefiles. Hence, the user can execute and manipulate those Makefiles in order to manually intervene in the particular steps, if he/she wants to do so. All this is done automatically when an integrate, unstable, or simulate command is encountered for the first time. You may also force compilation explicitly by the command
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 re-generated and re-compiled even if WHIZARD would think that this is unncessary. The same effect is achieved by calling WHIZARD with a command-line switch,
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 re-running 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 Stand-alone WHIZARD with precompiled processesOnce you have set up a process library, it is straightforward to make a special stand-alone WHIZARD executable which will have this library preloaded on startup. This is a matter of convenience, and it is also useful if you need a statically linked executable for reasons of profiling, batch processing, etc. For this task, there is a variant of the compile command:
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 center-of-momentum 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 single-beam events, where structure functions for single particles indeed do make sense. In the following sections we list the available options for structure functions or spectra inside WHIZARD and explain their usage. More about the physics of the implemented structure functions can be found in Chap. 9. 5.5.2 Asymmetric beams and Crossing 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 proton-proton 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 correspond to a beam of Z bosons with a momentum of 500 GeV. Note, however, that WHIZARD will always do the integration of the particle width in the particle’s rest frame, while the moving beam is then only taken into account for the frame of reference for the simulation. Further options then simply having different beam energies describe a non-vanishing between the two incoming beams. Such concepts are quite common e.g. for linear colliders to improve the beam properties in the collimation region at the beam interaction points. Such crossing angles can be specified in the beam setup, too, using the beams_theta command: beams = e1, E1 It is important that when a crossing angle is being specified, and the collision system consequently never is the center-of-momentum system, the beam momenta have to explicitly set. Besides a planar crossing angle, one is even able to rotate an azimuthal distance: beams = e1, E1 5.5.3 LHAPDFFor incoming hadron beams, the beams statement specifies which structure functions are used. The simplest example is the study of parton-parton scattering processes at a hadron-hadron collider such as LHC or Tevatron. The LHAPDF structure function set is selected by a syntax similar to the process setup, namely the example already shown above:
Note that there are slight differences in using the LHAPDF release series 6 and the older Fortran LHAPDF release series 5, at least concerning the naming conventions for the PDF sets 4. The above beams statement selects a default LHAPDF structure-function set for both proton beams (which is the CT10 central set for LHAPDF 6, and cteq6ll.LHpdf central set for LHAPDF5). The structure function will apply for all quarks, antiquarks, and the gluon as far as supported by the particular LHAPDF set. Choosing a different set is done by adding the filename as a local option to the lhapdf keyword:
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
integer-valued 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 = "<path-to-lhapdf>". Usually, it is not necessary to set this because WHIZARD detects this path via the lhapdf-config script during configuration, but in the case paths have been moved, or special files/special locations are to be used, the user can specify this location explicitly. 5.5.4 Built-in PDFsIn addition to the possibility of linking against LHAPDF, WHIZARD
comes with a couple of built-in 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 [73] for hadron-collider PDF structure functions and their manipulations are correctly linked to WHIZARD, it can be used for advanced calculations and simulations of hadron collider physics. Its main usage inside WHIZARD is for matching schemes between 4-flavor and 5-flavor schemes in b-parton initiated processes at hadron colliders. Note that in versions 2.2.0 and 2.2.1 it only worked together with LHAPDF version 5, while with the LHAPDF version 6 interface from version 2.2.2 on it can be used also with the modern version of PDFs from LHAPDF. Furthermore, from version 2.2.2, the HOPPET b parton matching also works for the builtin PDFs. It depends on the corresponding process and the energy scales involved whether it is a better description to use the 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 higher-order effect for the SLC and LEP I colliders. The soft-collinear and soft photon radiation can indeed be resummed/exponentiated to all orders in perturbation theory [7], while higher orders in hard-collinear photons have to be explicitly calculated order by order [8, 9]. WHIZARD has an intrinsic implementation of the lepton ISR structure function that includes all orders of soft and soft-collinear photons as well as up to the third order in hard-collinear photons. It can be switched on by the following statement:
As the ISR structure function is a single-beam 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 HERA-like setup:
Their are several options for the lepton-collider ISR structure function that are summarized in the following:
The maximal order of the hard-collinear photon emission taken into account by WHIZARD is set by the integer variable isr_order; the default is the maximally available order of three. With the variable isr_alpha, the value of the QED coupling constant αQED used in the ISR structure function can be set. The default is taken from the active physics model. The mass of the radiating lepton (in most cases the electron) is set by isr_mass; again the default is taken from the active physics model. Furthermore, the upper integration border for the ISR structure function which acts roughly as an upper hardness cutoff for the emitted photons, can be set through isr_q_max; if not set, the collider energy (possibly after beamstrahlung, cf. Sec. 5.5.7) √s (or √s) is taken. 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 non-collinear 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 pT distribution to those particles and their interaction products, i.e., all outgoing particles. Cuts that depend on photon pT may be applied to the modified events. For details on the ISR photon handler, cf. Sec. 10.4. The flag ?isr_recoil switches on pT recoil of the emitting lepton against photon radiation during integration; per default it is off. The flag ?isr_keep_energy controls the mode of on-shell projection for the splitting process with pT. 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, pT will become significant, and (i) energy/momentum non-conservation, 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 beam-beam interactions in lepton-collider physics is GuineaPig++ [10, 11, 12]. A direct interface between this tool GuineaPig++ and WHIZARD had existed as an inofficial add-on to the legacy branch WHIZARD1, but is no longer applicable in WHIZARD2. A WHIZARD-internal interface is foreseen for the very near future, most probably within this v2.2 release. Other options are to use parameterizations of the beam spectrum that have been included in the package CIRCE1 [6] which has been interfaced to WHIZARD since version v1.20 and been included in the WHIZARD2 release series. Another option is to generate a beam spectrum externally and then read it in as an ASCII data file, cf. Sec. 5.5.8. More about this can be found in a dedicated section on lepton collider spectra, Sec. 10.3. In this section, we discuss the usage of beamstrahlung spectra by means of the CIRCE1 package. The beamstrahlung spectra are true spectra, so they have to be applied to pairs of beams, and an application to only one beam is meaningless. They are switched on by this beams statement including structure functions:
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 semi-analytic approximation formulae for the spectra, or a Monte-Carlo sampling of the sampling. The second possibility always give a better description of the spectra, and is the default for WHIZARD. It can, however, be switched off by setting the flag ?circe1_generate to false. As the beamstrahlung spectra are sharply peaked at the collider energy, but still having long tails, a mapping of the spectra for an efficient phase-space sampling is almost mandatory. This is the default in WHIZARD, which can be changed by the flag ?circe1_map. Also, the default exponent for the mapping can be changed from its default value 2. with the variable circe1_mapping_slope. It is important to efficiently sample the peak position of the spectrum; the effective ratio of the peak to the whole sampling interval can be set by the parameter circe1_eps. The integer parameter circe1_chat sets the chattiness or verbosity of the CIRCE1 package, i.e. how many messages and warnings from the beamstrahlung generation/sampling will be issued. The actual internal version and revision of the CIRCE1 package are set by the two integer parameters circe1_ver and circe1_rev. The default is in any case always the newest version and revision, while older versions are still kept for backwards compatibility and regression testing. Finally, the geometry and design of the accelerator type is set with the string variable $circe1_acc: it contains the possible options for the old "SBAND" and "XBAND" setups, as well as the "TESLA" and JLC/NLC SLAC design "JLCNLC". The setups for the most important energies of the ILC as they are summarized in the ILC TDR [13, 14, 15, 16] are available as ILC. Beam spectra for the CLIC [18, 19, 20] linear collider are much more demanding to correctly simulate (due to the drive beam concept; only the low-energy modes where the drive beam is off can be simulated with the same setup as the abovementioned machines). Their setup will be supported soon in one of the upcoming WHIZARD versions within the CIRCE2 package. An example of how to generate beamstrahlung spectra with the help of the package CIRCE2 (that is also a part of WHIZARD) is this:
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 pre-generated beam spectrum for such a machine, and make it available for simulation within WHIZARD as an external ASCII data file. Such files basically contain only pairs of energy fractions of the nominal collider energy √s (x values). In WHIZARD they can be used in simulation with the following beams statement:
Note that beam spectra must always be pair spectra, i.e. they are automatically applied to both beam simultaneously. Beam spectra via external files are expected to reside in the current working directory. Alternatively, WHIZARD searches for them in the install directory of WHIZARD in share/beam-sim. There you can find an example file, uniform_spread_2.5%.dat for such a beam spectrum. The only possible parameter that can be set is the flag ?beam_events_warn_eof whose default is true. This triggers the issuing of a warning when the end of file of an external beam spectrum file is reached. In such a case, WHIZARD starts to reuse the same file again from the beginning. If the available data points in the beam events file are not big enough, this could result in an insufficient sampling of the beam spectrum. 5.5.9 Gaussian beam-energy 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 on-shell approximation for the e → eγ collinear splitting to allow the simulation of photon-induced backgrounds in lepton collider physics. The original concept is that of the Weizsäcker-Williams approximation [21, 22, 23]. This is a single-beam structure function that can be applied to both beams, or also to one beam only. Usually, there are some simplifications being made in the derivation. The formula which is implemented here and seems to be the best for the QCD background for low-pT hadrons, corresponds to Eq. (6.17) of Ref. [23]. As this reference already found, this leads to an "overshooting" of accuracy, and especially in the high-x (high-energy) region to wrong results. This formula corresponds to
Here, x is the ratio of the photon energy (called frequency ω in [23] over the original electron (or positron) beam energy E. The energy of the electron (or positron) after the splitting is given by x = 1−x. The simplified version is the one that corresponds to many publications about the EPA during SLC and LEP times, and corresponds to the q2 integration of Eq. (6.16e) in [23], where q2 is the virtuality or momentum transfer of the photon in the EPA:
While Eq. (1) is supposed to be the better choice for simulating hadronic background like low-pT hadrons and should be applied for the low-x region of the EPA, Eq. (2) seems better suited for high-x simulations like the photoproduction of BSM resonances etc. Note that the first term in Eqs. (1) and (2) is the standard Altarelli-Parisi QED splitting function of electron, Pe→ eγ(x) ∝ 1 + (1−x)2, while the last term in both equations is the default power correction. The two parameters Qmax2 and Qmin2 are the integration boundaries of the photon virtuality integration. Usually, they are given by the kinematic limits:
For low-pT hadron simulations, it is not a good idea to take the kinematic limit as an upper limit, but one should cut the simulation off at a hadronic scale like e.g. a multiple of the ρ mass. The user can switch between the two different options using the setting
or
for Eq. (1), while Eq. (2) can be chosen with
Note that a thorough study for high-energy e+e− colliders regarding the suitability of different EPA options is still lacking. For testing purposes also three more variants or simplifications of Eq. (2) are implemented: the first, steered by $epa_mode = log_power uses simply Qmax2 = s. This is also the case for the two other method. But the switch $epa_mode = log_simple uses just epa_mass (cf. below) as Qmin2. The final simplification is to drop the power correction, which can be chosen with $epa_mode = log. This corresponds to the simple formula:
Examples for the application of the EPA in WHIZARD are:
or for a single beam:
The last process allows the reaction of (quasi-) on-shell photons with protons. In the following, we collect the parameters and flags that can be adjusted when using the EPA inside WHIZARD:
The adjustable parameters are partially similar to the parameters in the QED initial-state radiation (ISR), cf. Sec. 5.5.6: the parameter epa_alpha sets the value of the electromagnetic coupling constant, αQED used in the EPA structure function. If not set, this is taken from the value inside the active physics model. The same is true for the mass of the particle that radiates the photon of the hard interaction, which can be reset by the user with the variable epa_mass. There are two dimensionful scale parameters, the minimal momentum transfer to the photon, epa_q_min, which must not be zero, and the upper momentum-transfer 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 soft-photon regulator for the splitting function in x space, epa_x_min, which also has to be explicitly set different from zero. Hence, a minimal viable scenario that will be accepted by WHIZARD looks like this:
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 trade-off between quality of the description on the one hand and numerical stability and speed on the other hand. In the case, the EPA is set after a second structure function like a hadron collider PDF, there is a flavor summation over the quark constituents inside the proton, which are then the radiating fermions for the EPA. Here, the masses of all fermions have to be identical. More about the physics of the equivalent photon approximation can be found in Chap. 9. 5.5.11 Effective W 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 high-energy weak vector-boson fusion and scattering processes at hadron colliders, particularly the Superconducting Super-Collider (SSC). This was at a time when the simulation of 2→ 4 processes war still very challenging and 2→ 6 processes almost impossible, such that this approximation was the only viable solution for the simulation of processes like 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 single-beam structure function that can be applied to one or both beams. We only give an example for both beams here, this is for a 3 TeV CLIC collider:
And this is for LHC or a higher-energy follow-up collider (which also shows the concatenation of the single-beam structure functions, applied to both beams consecutively, cf. Sec. 5.5.14:
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 pT of the electroweak boson, that can be set via eta_pt_max. Indeed, the transversal W/Z structure function is logarithmically divergent in that variable. If it is not set by the user, it is estimated from √s and the splitting kinematics. For the EWA, there is a flag to switch on a recoil for the electroweak boson against the radiating fermion, ?ewa_recoil. Note that this is an experimental feature that is not completely tested. In any case, the non-collinear kinematics violates 4-four momentum conservation, so there are two choices: either to conserve the energy (?ewa_keep_energy = true) or to conserve 3-momentum (?ewa_keep_energy = false). Momentum conservation for the kinematics is the default. This is due to the fact that for energy conservation, there will be a net total momentum in the event including the beam remnants (ISR/EPA/EWA radiated particles) that leeds to unexpected or unphysical features in the energy distributions of the beam remnants recoiling against the rest of the event. More details about the physics can be found in Chap. 9. 5.5.12 Energy scans using structure 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 fall-off 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 lets the scan start at 50 GeV. There are no adjustable parameters for this structure function. How to plot the invariant mass distribution of the final-state muon pair to show the energy scan over the cross section, will be explained in Sec. 5.9. More details can be found in Chap. 9. 5.5.13 Photon collider spectraOne option that has been discussed as an alternative possibility for a high-energy linear lepton collider is to convert the electron and positron beam via Compton backscattering off intense laser beams into photon beams [24, 25, 26]. Naturally, due to the production of the photon beams and the inherent electron spectrum, the photon beams have a characteristic spectrum. The simulation of such spectra is possible within WHIZARD by means of the subpackage CIRCE2, which have been mentioned already in Sec. 5.5.7. It allows to give a much more elaborate description of a linear lepton collider environment than CIRCE1 (which, however, is not in all cases necessary, as the ILC beamspectra for electron/positrons can be perfectly well described with CIRCE1). Here is a typical photon collider setup where we take a photon-initiated process:
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 single-beam structure functions or pair spectra. One important thing is whether there is a phase-space mapping for these structure functions. Also, there are some combinations which do not make sense from the physics point of view, for example using lepton-collider ISR for protons, and then afterwards switching on PDFs. Such combinations will be vetoed by WHIZARD, and you will find an error message like (cf. also Sec. 4.3): ******************************************************************************
******************************************************************************
*** FATAL ERROR: Beam structure: [....] not supported
******************************************************************************
******************************************************************************
Common examples for the concatenation of structure functions are linear collider applications, where beamstrahlung (macroscopic electromagnetic beam-beam interactions) and electron QED initial-state radiation are both switched on:
Another possibility is the simulation of photon-induced 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 high-energy vector boson scattering:
Note that this last case involves a flavor sum over the five active quark (and anti-quark) species u, d, c, s, b in the proton, all of which act as radiators for the electroweak vector bosons in the EWA. This would be an example with three structure functions:
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 statement5:
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 % left-handed electron polarization means that 80 % of the electron beam are polarized, 20 % are unpolarized, i.e. 20 % have half left- and half right-handed polarization each. Hence, 90 % of the electron beam is left-handed, 10 % is right-handed. How does the specification of the polarization work? If there are no entries at all in the polarization constructor, @(), the beam is unpolarized, and the spin density matrix is proportional to the unit/identity matrix. Placing entries into the @() constructor follows the concept of sparse matrices, i.e. the entries that have been specified will be present, while the rest remains zero. Single numbers do specify entries for that particular helicity on the main diagonal of the spin density matrix, e.g. for an electron @(-1) means (100%) left-handed polarization. Different entries are separated by commas: @(1,-1) sets the two diagonal entries at positions (1,1) and (−1,−1) in the density matrix both equal to one. Two remarks are in order already here. First, note that you do not have to worry about the correct normalization of the spin density matrix, WHIZARD is taking care of this automatically. Second, in the screen output for the beam data, only those entries of the spin density matrix that have been specified by the user, will be displayed. If a beams_pol_fraction statement appears, other components will be non-zero, but might not be shown. E.g. ILC-like, 80 % polarization of the electrons, 30 % positron polarization will be specified like this for left-handed electrons and right-handed positrons:
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):
Off-diagonal entries that are equal to one (up to the normalization) of the spin-density matrix can be specified simply by the position, namely: @(m:m′, m″). This would result in a spin density matrix with diagonal entry 1 for the position (m″, m″), and an entry of 1 for the off-diagonal position (m,m′). Furthermore, entries in the density matrix different from 1 with a numerical value <val> can be specified, separated by another colon: @(m:m′:<val>). Here, it does not matter whether m and m′ are different or not. For m = m′ also diagonal spin density matrix entries different from one can be specified. Note that because spin density matrices have to be Hermitian, only the entry (m,m′) has to be set, while the complex conjugate entry at the transposed position (m′,m) is set automatically by WHIZARD. We will give some general density matrices now, and after that a few more definite examples. In the general setups below, we always give the expression for the spin density matrix only for one single beam.
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 multi-channel Monte-Carlo integration. The integration, in turn, requires a phase-space setup, i.e., a collection of phase-space channels, which are mappings of the unit hypercube onto the complete space of multi-particle kinematics. This phase-space information is encoded in the file xxx.phs, where xxx is the process tag. WHIZARD generates the phase-space file on the fly and can reuse it in later integrations. For each phase-space channel, the unit hypercube is binned in each dimension. The bin boundaries are allowed to move during a sequence of iterations, each with a fixed number of sampled phase-space points, so they adapt to the actual phase-space density as far as possible. In addition to this intrinsic adaptation, the relative channel weights are also allowed to vary. All these steps are done automatically when the integrate command is executed. At the end of the iterative adaptation procedure, the program has obtained an estimate for the integral of the matrix element over phase space, together with an error estimate, and a set of integration grids which contains all information on channel weights and bin boundaries. This information is stored in a file xxx.vg, where xxx is the process tag, and is used for event generation by the simulate command. 5.7.1 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 matrix-element function. For each iteration, WHIZARD computes an estimate of the integral and an estimate of the error, based on the binned sums of matrix element values and squares. It also computes an estimate of the rejection efficiency for generating unweighted events, i.e., the ratio of the average sampling function value over the maximum value of this function. After each iteration, both the integration grids (the binnings) and the relative weights of the integration channels can be adapted to minimize the variance estimate of the integral. After each pass of several iterations, WHIZARD computes an average of the iterations within the pass, the corresponding error estimate, and a χ2 value. The integral, error, efficiency and χ2 value computed for the most recent integration pass, together with the most recent integration grid, are used for any subsequent calculation that involves this process, in particular for event generation. In the default setup, during the first pass(es) both grid binnings and channel weights are adapted. In the final (usually second) pass, only binnings are further adapted. Roughly speaking, the final pass is the actual calculation, while the previous pass(es) are used for “warming up” the integration grids, without using the numerical results. Below, in the section about the specification of the iterations, Sec. 5.7.3, we will explain how it is possible to change the behavior of adapting grids and weights. Here is an example of the integration output, which illustrates these properties. The SINDARIN script describes the process e+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 pT 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 √ncalls. The accuracy value has the property that it is independent of ncalls, it describes the quality of adaptation of the current grids. Good-quality grids have a number of order one, the smaller the better. The next column is the estimate for the rejection efficiency in percent. Here, the value should be as high as possible, with 100 % being the possible maximum. In the example, the grids are adapted over ten iterations, after which the accuracy and efficiency have saturated at about 1.0 and 10 %, respectively. The asterisk in the accuracy column marks those iterations where an improvement over the previous iteration is seen. The average over these iterations exhibits an accuracy of 1.22, corresponding to 0.39 % error, and a χ2 value of 1.15, which is just right: apparently, the phase-space for this process and set of cuts is well-behaved. The subsequent five iterations are used for obtaining the final integral, which has an accuracy below one (error 0.3 %), while the efficiency settles at about 10 %. In this example, the final χ2 value happens to be quite small, i.e., the individual results are closer together than the error estimates would suggest. One should nevertheless not scale down the error, but rather scale it up if the χ2 result happens to be much larger than unity: this often indicates sub-optimally adapted grids, which insufficiently map some corner of phase space. One should note that all values are subject to statistical fluctuations, since the number of calls within each iterations is finite. Typically, fluctuations in the efficiency estimate are considerably larger than fluctuations in the error/accuracy estimate. Two subsequent runs of the same script should yield statistically independent results which may differ in all quantities, within the error estimates, since the seed of the random-number generator will differ by default. It is possible to get exactly reproducible results by setting the random-number seed explicitly, e.g.,
at any point in the SINDARIN script. seed is a predefined intrinsic variable. The value can be any 32bit integer. Two runs with different seeds can be safely taken as statistically independent. In the example above, no seed has been set, and the seed has therefore been determined internally by WHIZARD from the system clock. The concluding line with the time estimate applies to a subsequent simulation step with unweighted events, which is not actually requested in the current example. It is based on the timing and efficiency estimate of the most recent iteration. As a default, a graphical output of the integration history will be produced (if both LATEX and MetaPost have been available during configuration). Fig. 5.1 shows how this looks like, and demonstrates how a proper convergence of the integral during the adaptation looks like. The generation of these graphical history files can be switched off using the command ?vis_history = false. 5.7.2 Integration run 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 re-run: all results that have been overwritten have to be recreated. To avoid this, the user may identify a specific run by a string-valued ID, e.g.
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 final-state particles. They are guesses for values that yield good-quality grids and error values in standard situations, where no exceptionally strong peaks or loose cuts are present in the integrand. Actually, the large number of warmup iterations in the previous example indicates some safety margin in that respect. It is possible, and often advisable, to adjust the iteration and call numbers to the particular situation. One may reduce the default numbers to short-cut the integration, if either less accuracy is needed, or CPU time is to be saved. Otherwise, if convergence is bad, the number of iterations or calls might be increased. To set iterations manually, there is the iterations command:
This is a comma-separated list. Each pair of values corresponds to an integration pass. The value before the colon is the number of iterations for this pass, the other number is the number of calls per iteration. While the default number of passes is two (one for warmup, one for the final result), you may specify a single pass
where the relative channel weights will not be adjusted (because this is the final pass). This is appropriate for well-behaved integrands where weight adaptation is not necessary. You can also define more than two passes. That might be useful when reusing a previous grid file with insufficient quality: specify the previous passes as-is, so the previous results will be read in, and then a new pass for further adaptation. In the final pass, the default behavior is to not adapt grids and weights anymore. Otherwise, different iterations would be correlated, and a final reliable error estimate would not be possible. For all but the final passes, the user can decide whether to adapt grids and weights by attaching a string specifier to the number of iterations: "g" does adapt grids, but not weights, "w" the other way round. "gw" or "wg" does adapt both. By the setting "", all adaptations are switched off. An example looks like this:
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 phase-space configuration for the process at hand. The method for the phase-space parameterization is determined by the string variable $phs_method. At the moment there are only two options, "single", for testing purposes, that is mainly used internally, and WHIZARD’s traditional method, "wood". This parameterization is particularly adapted and fine-tuned for electroweak processes and might not be the ideal for for pure jet cross sections. In future versions of WHIZARD, more options for phase-space parameterizations will be made available, e.g. the RAMBO algorithm and its massive cousin, and phase-space parameterizations that take care of the dipole-like emission structure in collinear QCD (or QED) splittings. For the standard method, the phase-space parameterization is laid out in an ASCII file <process-name>_i<comp>.phs. Here, <process-name> is the process name chosen by the user while <comp> is the number of the process component of the corresponding process. This immediately shows that different components of processes are getting different phase space setups. This is necessary for inclusive processes, e.g. the sum of 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 phase-space file and inspect and/or
modify it before proceeding further. To this end, there is the parameter
Things might go wrong with the default phase-space generation, or manual intervention might be necessary to improve later performance. There are a few parameters that control the algorithm of phase-space generation. To understand their meaning, you should realize that phase-space parameterizations are modeled after (dominant) Feynman graphs for the current process. The main phase space setup woodFor the main phase-space parameterization of WHIZARD, which is called "wood", there are many different parameters and flags that allow to tune and customize the phase-space setup for every certain process: The parameter 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.3. 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 right-hand side of the assignment. It may be defined in global scope, so it is applied to all subsequent processes. It may be redefined by another cuts statement. This overrides the first cuts setting: the cuts statement is not cumulative. Multiple cuts should be specified by the logical operators of SINDARIN, for instance
Cuts may also be defined local to an integrate command, i.e., in the options in braces. They will apply only to the processes being integrated, overriding any global cuts. The right-hand side expression in the cuts statement is evaluated at the point where it is used by an integrate command (which could be an implicit one called by simulate). Hence, if the logical expression contains parameters, such as
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
process-independent manner, even if it is not required. The include
command allows for storing a set of cuts in a separate SINDARIN script which
may be read in anywhere. As an example, the system directories contain a file
5.7.6 QCD scale and couplingWHIZARD treats all physical parameters of a model, the coefficients in the Lagrangian, as constants. As a leading-order program, WHIZARD does not make use of running parameters as they are described by renormalization theory. For electroweak interactions where the perturbative expansion is sufficiently well behaved, this is a consistent approach. As far as QCD is concerned, this approach does not yield numerically
reliable results, even on the validity scale of the tree approximation.
In WHIZARD2, it is therefore possible to replace the fixed value of
αs (which is accessible as the intrinsic model variable
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 user-defined scale, the special object scale. This is assigned and used just like the cuts object. The right-hand side is a real-valued expression. Here is an example:
This selects the pT value of the first entry in the list of colored particles sorted by decreasing pT, i.e., the pT of the hardest jet. The scale definition is used not just for running αs (if enabled), but it is also the factorization scale for the LHAPDF structure functions. These two values can be set differently by specifying factorization_scale for the scale at which the PDFs are evaluated. Analogously, there is a variable renormalization_scale that sets the scale value for the running αs. Whenever any of these two values is set, it supersedes the scale value. Just like the cuts expression, the expressions for scale, factorization_scale and also renormalization_scale are evaluated at the point where it is read by an explicit or implicit integrate command. 5.7.7 Reweighting factorIt is possible to reweight the integrand by a user-defined function of the event kinematics. This is done by specifying a weight expression. Syntax and usage is exactly analogous to the scale expression. Example:
We should note that the phase-space setup is not aware of this reweighting, so in complicated cases you should not expect adaptation to achieve as accurate results as for plain cross sections. Needless to say, the default weight is unity. 5.8 EventsAfter the cross section integral of a scattering process is known (or the partial-width integral of a decay process), WHIZARD can generate event samples. There are two limiting cases or modes of event generation:
The choice of parameterizations and the iterative adaptation of the integration grids gradually shift the generation mode from option 2 to option 1, which obviously is preferred since it simulates the actual outcome of an experiment. Unfortunately, this adaptation is perfect only in trivial cases, such that the Monte-Carlo integration yields non-uniform probability still with weighted events. Unweighted events are obtained by rejection, i.e., accepting an event with a probability equal to its own weight divided by the maximal possible weight. Furthermore, the maximal weight is never precisely known, so this probability can only be estimated. The default generation mode of WHIZARD is unweighted. This is controlled by
the parameter 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 machine-dependent binary format, so it is not of
immediate use. Its principal purpose is to serve as a cache: if you re-run
the same script, before starting simulation, it will look for an existing
event file that matches the input. If nothing has changed, it will find the
file previously generated and read in the events, instead of generating them.
Thus you can modify the analysis or any further steps without repeating the
time-consuming task of generating a large event sample. If you change the
number of events to generate, the program will make use of the existing event
sample and generate further events only when it is used up. If necessary, you
can suppress the writing/reading of the raw event file by the parameters
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 built-in capabilities, to produce tables, histograms, or calculate inclusive observables. The basic analysis features of WHIZARD are described below in Sec. 5.9. It can be written to an external file in a standard format that a human or an external program can understand. In Chap. 11, you will find a more thorough discussion of event generation with WHIZARD, which also covers in detail the available event-file formats. Finally, WHIZARD can rescan an existing event sample. The event sample may either be the result of a previous simulate run or, under certain conditions, an external event sample produced by another generator or reconstructed from data.
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 final-state particles
are on-shell and stable, so they appear explicitly in the generated event
files. If events are desired where the decay products of the W bosons
appear, one has to generate another process, e.g., pp→ udūd. In
this case, the intermediate vector bosons, if reconstructed, are off-shell as
dictated by physics, and the process contains all intermediate states that are
possible. In this example, the matrix element contains also ZZ, photon, and
non-resonant intermediate states. (This can be restricted via the
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 computing-intensive. The factorization neglects all off-shell effects and irreducible background diagrams that do not have the decaying particles as an intermediate resonance. While WHIZARD is able to deal with multi-particle processes without factorization, the needed computing resources rapidly increase with the number of external particles. Particularly, it is the phase space integration that becomes the true bottleneck for a high multiplicity of final state particles. In order to use the factorized approach, one has to specify particles as unstable. (Also, the ?allow_decays switch must be true; this is however its default value.) We give an example for a 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 narrow-width approximation. Theoretically, this restriction is necessary since whenever the width plays an important role, the usage of the factorized approach will not be fully justified. In particular, all involved matrix elements must be evaluated on-shell, or otherwise gauge-invariance issues could spoil the calculation. (There are plans for a future WHIZARD version to also include Breit-Wigner or Gaussian distributions when using the factorized approach.) Decays can be concatenated, e.g. for top pair production and decay, e+ 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 anti-top 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 <particle-tag>. 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 branching-fraction
values for the decays of a particle, for instance, NLO-corrected
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 branching-ratio 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 production-decay factorization, there is support for a simplified handling of particle decays. Essentially, there are four options, taking a decay W_ud: W−→ ū d as an example:
Here, the isotropic option completely eliminates spin correlations. The diagonal-decays option eliminates just the off-diagonal entries of the W spin-density matrix. This is equivalent to a measurement of spin before the decay. As a result, spin correlations are still present in the classical sense, while quantum coherence is lost. The definite-helicity option is similar and additional selects only the specified helicity component for the decaying particle, so its decay distribution assumes the shape for an accordingly polarized particle. All options apply in the rest frame of the decaying particle, with the particle’s momentum as the quantization axis. Automatic 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 machine-dependent event format. There are a series of human-readable ASCII event formats that are supported: very verbose formats intended for debugging, formats that have been agreed upon during the Les Houches workshops like LHA and LHEF, or formats that are steered through external packages like HepMC. More details about event formats can be found in Sec. 11.5. 5.9 Analysis and VisualizationSINDARIN natively supports basic methods of data analysis and visualization which are frequently used in high-energy physics studies. Data generated during script execution, in particular simulated event samples, can be analyzed to evaluate further observables, fill histograms, and draw two-dimensional plots. So the user does not have to rely on his/her own external graphical analysis method (like e.g. gnuplot or ROOT etc.), but can use methods that automatically ship with WHIZARD. In many cases, the user, however, clearly will use his/her own analysis machinery, especially experimental collaborations. In the following sections, we first summarize the available data structures, before we consider their graphical display. 5.9.1 ObservablesAnalyses in high-energy physics often involve averages of quantities other than a total cross section. SINDARIN supports this by its observable objects. An observable is a container that collects a single real-valued variable with a statistical distribution. It is declared by a command of the form observable analysis-tag where analysis-tag is an identifier that follows the same rules as a variable name. Once the observable has been declared, it can be filled with values. This is done via the record command: record analysis-tag (value) To make use of this, after values have been filled, we want to perform the actual analysis and display the results. For an observable, these are the mean value and the standard deviation. There is the command write_analysis: write_analysis (analysis-tag) Here is an example:
The result is displayed on screen:
5.9.2 The analysis expressionThe most common application is the computation of event observables – for instance, a forward-backward asymmetry – during simulation. To this end, there is an analysis expression, which behaves very similar to the cuts expression. It is defined either globally analysis = logical-expr or as a local option to the simulate or rescan commands which generate and handle event samples. If this expression is defined, it is not evaluated immediately, but it is evaluated once for each event in the sample. In contrast to the cuts expression, the logical value of the analysis expression is discarded; the expression form has been chosen just by analogy. To make this useful, there is a variant of the record command, namely a record function with exactly the same syntax. As an example, here is a calculation of the forward-backward symmetry in a process ee_mumu with final state µ+µ−:
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 analysis-tag (lower-bound, upper-bound) This creates a histogram data structure for an (unspecified) observable. The entries are organized in bins between the real values lower-bound and upper-bound. The number of bins is given by the value of the intrinsic integer variable n_bins, the default value is 20. The histogram declaration supports an optional argument, so the number of bins can be set locally, for instance histogram pt_distribution (0 GeV, 500 GeV) { n_bins = 50 } Sometimes it is more convenient to set the bin width directly. This can be done in a third argument to the histogram command. histogram pt_distribution (0 GeV, 500 GeV, 10 GeV) If the bin width is specified this way, it overrides the setting of n_bins. The record command or function fills histograms. A single call record (real-expr) puts the value of real-expr into the appropriate bin. If the call is issued during a simulation where unweighted is false, the entry is weighted appropriately. If the value is outside the range specified in the histogram declaration, it is put into one of the special underflow and overflow bins. The write_analysis command prints the histogram contents as a table in
blank-separated fixed columns. The columns are: x (bin midpoint), y (bin
contents), Δ y (error), excess weight, and n (number of entries).
The output also contains comments initiated by a 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 analysis-tag Like observables and histograms, plots are filled by the record command or expression. To this end, it can take two arguments, record (x-expr, y-expr) or up to four: record (x-expr, y-expr, y-error) 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 process-specific information to a file ⟨tag⟩.log, where ⟨tag⟩ is the process name. Furthermore, the write_analysis command dumps analysis data – tables for histograms and plots – to its own set of files, cf. Sec. 5.9. However, there is the occasional need to write data to extra files in a custom format. SINDARIN deals with that in terms of the following commands: 5.10.1 Output 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 ⟨data-objects⟩ according to ⟨format-string-expr⟩ and print the resulting string to standard output if the string variable $out_file is undefined. If $out_file is defined and the file with this name is open for writing, print to this file instead. Print a newline at the end if ?out_advance is true, otherwise don’t finish the line. The ⟨format-string-expr⟩ must evaluate to a string. Formatting follows a subset of the rules for the printf(3) command in the C language. The supported rules are:
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 next-to-leading 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 matrix-element generators that are BLHA-compatible. 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 one-loop matrix elements are the scientific algebra program form and the generator of loop topologies and diagrams, qgraf. These can be accessed via their respective webpages http://www.nikhef.nl/~form/ and http://cfif.ist.utl.pt/~paulo/qgraf.html. Note also that both Java and the Java runtime environment have to be installed in order for Gosam to properly work. Furthermore, libtool needs to be installed. 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 --with-openloops=<openloops_dir>. OpenLoops can be checked out with
Note that WHIZARD only supports OpenLoops version that are at least 2.1.1 or newer. Alternatively, one can use the public beta version of OpenLoops, which can be checked out by the command
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+1-particle tree-level matrix elements in the process directory, so that a complete NLO calculation including real amplitudes can be performed only with OpenLoops. The libraries can then be installed via
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 ppllj, 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 one-loop diagrams can then be used with the SINDARIN option
The functional tests which check the OpenLoops functionality require the libraries ppllj, 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.3.0. Recola can be linked to WHIZARD during configuration by including
In case the Recola library is not in a standard path or a path accessible in the LD_LIBRARY_PATH (or DYLD_LIBRARY_PATH) of the operating system, then the option
can be set, 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 NLO-components to be integrated, i.e. born, real, virtual, dglap, (for hadron collisions) or mismatch (for the soft mismatch in resonance-aware computations) and full. The full option switches on all components and is required if the total NLO result is desired. For example, specifying
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 resonance-aware FKS subtraction,
This mode comes with an additional integration component, the so-called soft mismatch. Note that you can modify the number of iterations in each component with the multipliers:
5.11.3 Fixed-order NLO eventsFixed-order NLO events can also be produced in three different modes: Combined weighted, combined unweighted and separated weighted.
Weighted fixed-order 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 Nphs + 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 events 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. (5) produces either no or one additional emission. These events can then either be analyzed directly or passed on to the parton shower6 for the full simulation. You activate this with
The pTmin of Eq. (5) 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_max. Note that the POWHEG event generation needs an additional grid for efficient event generation that is generated during integration if ?powheg_matching = true is set. Thus, this needs to be set before the integrate statement. Further options that steer the efficiency of this grid are powheg_grid_size_xi, powheg_grid_size_y and powheg_grid_sampling_points. 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 string variable
The finite contributions, which by definition should not contain soft or collinear emissions, will then integrate like an 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
|
