Lilith: a tool for constraining new physics from Higgs measurements

The properties of the observed Higgs boson with mass around 125 GeV can be affected in a variety of ways by new physics beyond the Standard Model (SM). The wealth of experimental results, targeting the different combinations for the production and decay of a Higgs boson, makes it a non-trivial task to assess the compatibility of a non-SM-like Higgs boson with all available results. In this paper we present Lilith, a new public tool for constraining new physics from signal strength measurements performed at the LHC and the Tevatron. Lilith is a Python library that can also be used in C and C++/ROOT programs. The Higgs likelihood is based on experimental results stored in an easily extensible XML database, and is evaluated from the user input, given in XML format in terms of reduced couplings or signal strengths. The results of Lilith can be used to constrain a wide class of new physics scenarios.


Introduction
The discovery of a Higgs boson with properties compatible with those of the SM and mass around 125 GeV at CERN's Large Hadron Collider (LHC) [1,2] was a major breakthrough. Indeed, the Higgs boson was the last elementary particle predicted by the SM remaining to be observed. But, more importantly, the Higgs field has a key role in the SM as it triggers the breaking of the electroweak symmetry and gives masses to the elementary particles. Precision measurements of the properties of the observed boson are of utmost importance to assess its role in the breaking of the electroweak symmetry. They could reveal a more complicated Higgs sector, indicating the presence of more elementary scalars or compositeness of the observed particle, and could also shed light on a large variety of beyond-the-SM (BSM) particles that couple to the Higgs boson. Conversely, precision measurements can be used to rule out new physics scenarios affecting the properties of the Higgs boson.
That the mass of the observed Higgs boson is about 125 GeV is a fortunate coincidence as many decay modes of the SM Higgs boson are accessible with a modest integrated luminosity at the LHC [3]. Hence, complementary information on the properties of the Higgs boson were already obtained from the measurements performed during Run I of the LHC at 7-8 TeV center-of-mass energy [4,5]. A large variety of models of new physics (both effective and explicit ones) can be constrained from the measurements presented in terms of signal strengths. These results were used in a large number of phenomenological studies in the past three years (see Refs.  for a sample of recent studies based on the full data collected at Run I).
However, it is not straightforward to put constraints on new physics from the measured signal strengths. Indeed, a large number of analyses have already been performed by the ATLAS and CMS collaborations. They usually include several event categories, and present signal strength results in different ways. Extracting all necessary information from the figures of the various publications is a tedious and lengthy task. Moreover, as the full statistical models used by the experimental collaborations are not public, a number of assumptions need to be made for constructing a likelihood. The validity of these approximations should be assessed from a comparison with the results provided by ATLAS and CMS.
In order to put constraints on new physics from the LHC Higgs results, many groups have been developing private codes. Moreover, recently a public tool, HiggsSignals [32], became available. HiggsSignals is a FORTRAN code that uses the signal strengths for individual measurements, taking into account the associated efficiencies. In this paper, we present a new public tool, Lilith. 1 Lilith is a library written in Python, that can easily be used in any Python script as well as in C and C++/ROOT codes, and for which we also provide a command-line interface. It follows a different approach than HiggsSignals in that it uses as a primary input results in which the fundamental production and decay modes are unfolded from experimental categories. The experimental results are stored in XML files, making it easy to modify and extend. The user input can be given in terms of reduced couplings or signal strengths, also specified in an XML format.
In Section 2, we present the signal strength framework used to encode deviations from the SM at the LHC, as well as the experimental results that we use as input in Lilith. The parametrization of new physics effects on the observed Higgs boson, as well as derivation of signal strengths, are presented in Section 3. All technical details on how to use Lilith and the XML formats that we use are then given in Section 4. Constraints derived from Lilith are validated in Section 5, and two concrete examples of its capabilities are given in Section 6. Finally, prospects for Run II of the LHC are discussed in Section 7, and conclusions are given in Section 8.

From experimental results to likelihood functions 2.1 Signal strength measurements
Thanks to the excellent operation of the LHC and to the wealth of accessible final states for a 125 GeV SM-like Higgs boson, the properties of the observed Higgs boson have been measured with unforeseeable precision by the ATLAS and CMS collaborations already during Run I of the LHC at 7-8 TeV center-of-mass energy [4,5]. LHC searches are targeting the different combinations for the production and decay modes of a Higgs boson. The SM Higgs boson has five main production mechanisms at a hadron collider: gluon fusion (ggH), vector-boson fusion (VBF), associated production with an electroweak gauge boson (WH and ZH, collectively denoted as VH) and associated production with a pair of top quarks (ttH). 2 Observation of these production modes constrains the couplings of the Higgs to vector bosons (VBF, VH) and to third-generation quarks (ggH, ttH). The main decay modes accessible at the LHC are H → γγ, H → ZZ * → 4 , H → W W * → 2 2ν, H → bb and H → τ τ (with ≡ e, µ). They can provide complementary information on the couplings of the Higgs to vector bosons (from the decay into ZZ * , W W * , and γγ) and to third-generation fermions (from the decay into bb, τ τ , and γγ). Being loop-induced processes, gg → H and H → γγ also have sensitivity to BSM colored particles and BSM electrically-charged particles, respectively.
The results of the Higgs searches at the LHC are given in terms of signal strengths, µ, which scale the number of signal events expected for the SM Higgs, n s . For a given set of selection criteria, the expected number of events is therefore µ · n s + n b , where n b is the expected number of background events, so that µ = 0 corresponds to the no-Higgs scenario and µ = 1 to a SM-like Higgs. Equivalently, signal strengths can be expressed as where A × ε is the product of the acceptance and of the efficiency of the selection criteria. Two assumptions can subsequently be made: first, the signal is a sum of processes that exist for a 125 GeV SM Higgs boson, i.e. σ = X,Y σ(X) B(H → Y ) for the various production modes X ∈ (ggH, VBF, VH, ttH) and decay modes Y ∈ (γγ, ZZ * , W W * , bb, τ τ , . . .). Second, the acceptance times efficiency is identical to the SM one for all processes, that is (A × ε) X,Y = [(A × ε) X,Y ] SM for every X and Y . Under these conditions, signal strengths read where the eff X,Y are "reduced efficiencies", corresponding to the relative contribution of each combination for the production and decay of a Higgs boson to the signal. These can be estimated from the A×ε obtained in a Monte Carlo simulation of individual processes. In the case of an inclusive search targeting a given decay mode Y (i.e. ∀X, (A × ε) X,Y = (A × ε) Y ), eff Y is equal to the ratio of SM cross sections, σ SM X /( X σ SM X ). The signal strength framework used by the ATLAS and CMS collaborations is based on the general form of Eq. (2), hence on the assumption that new physics results only in the scaling of SM Higgs processes. This makes it possible to combine the information from various Higgs searches and assess the compatibility of given scalings of SM production and/or decay processes from a global fit to the Higgs data. This framework is very powerful as it can be used to constrain a wide variety of new physics models (some examples can be found in Ref. [33]). This is the approach that we will follow in Lilith. However, in order to derive constraints on new physics, one first needs to construct a likelihood function from the signal strength information given in the experimental publications. In particular, combining the results from several Higgs searches is non-trivial and deserves scrutiny.

Event categories versus unfolded production and decay modes
The searches for the Higgs boson performed by the ATLAS and CMS collaborations are divided into individual analyses usually focusing on a single decay mode. Within each analysis several event categories are then considered. Among other reasons, these are designed to optimize the sensitivity to the different production mechanisms of the SM Higgs boson (hence, they have different reduced efficiencies eff X,Y ). In order to put constraints on new physics from the results in a given event category, one needs to extract the measurement of the signal strength and the relevant eff X,Y information from the experimental publication. For example, results of the CMS H → γγ analysis [34], in terms of signal strengths for all categories, are shown on the left panel of Fig. 1. With the addition of the reduced efficiencies eff X,γγ , also given in Ref. [34], combinations of σ(X) B(H → γγ) can be constrained.
However, several problems arise when constructing a likelihood. First of all, as can be seen on the left panel of Fig. 1, only two pieces of information are given: the best fit to the data, that will be denoted asμ in the following, and the 68% confidence level (CL) interval or 1σ interval. The full likelihood function category per category is never provided by the experimental collaborations. Assuming that the measurements are approximately Gaussian, it is however possible to reconstruct a simple likelihood, L(µ), from this information. In that case, −2 log L(µ) follows a χ 2 law. From the boundaries of the 68% CL interval, left and right uncertainties at 68% CL, ∆µ − and ∆µ + , with respect to the best fit point can be derived. The likelihood can then be defined as with ∆µ − = ∆µ + in the Gaussian regime. While this is often a valid approximation to the likelihood, it should be pointed out that signal strength measurements are not necessarily Gaussian, depending in particular on the size of the event sample.
Barring this limitation, Eq. (3) can be used to constrain new physics. However, it requires that at least the 68% CL interval and the relevant reduced efficiencies eff X,Y are provided by the experimental collaboration for every individual category. This is not always the case. Categories are sometimes defined without giving the corresponding signal efficiencies (as in, e.g., the CMS ttH analysis [35]), and/or the result is given for a (set of) combined signal strength(s) but not in terms of signal strengths category per category (as in the ATLAS ZZ * and τ τ analyses [36,37] and in the CMS ttH analysis [35]). Such combined µ should in general not be used because they have been obtained under the assumption of SM-like production or decay of the Higgs boson. Whenever the eff X,Y are not given in the experimental publications it is in principe possible to obtain estimates from a reproduction of the selection criteria applied on signal samples generated by Monte Carlo simulation. However, this turns out to be a very difficult or impossible task. Indeed, searches for the Higgs boson typically rely on complex search strategies to optimize the sensitivity, such as multivariate analysis techniques that are impossible to reproduce in practice with the information currently available. Whenever the information on reduced efficiencies is not available, we are left to guesswork, with a natural default choice being that eff X is equal to the ratio of SM cross sections, σ SM X /( X σ SM X ), which would correspond to a fully inclusive search. Constraining new physics from a single LHC Higgs category can already be a non-trivial task and come with some uncertainty because the full information is not provided category per category. However, more severe complications typically arise when using several categories/searches at the same time, as is needed for a global fit to the Higgs data. The simplest solution is to define the full likelihood as the product of individual likelihoods, However, this assumes that all measurements are completely independent. We know that this is not the case as the various individual measurements share common systematic uncertainties. They are divided into two categories: the shared experimental uncertainties, coming from the presence of the same final state objects and from the estimation of the luminosity, and the shared theoretical uncertainties, dominated by the contributions from identical production and/or decay modes to the expected Higgs signal in different categories [38]. The estimation of the experimental uncertainties in ATLAS should be largely independent from the one in CMS, hence these correlations can be treated separately for measurements performed by one collaboration or the other. Conversely, theoretical uncertainties are estimated in the same way in ATLAS and CMS and should be correlated between all measurements.
In the case where all measurements are well within the Gaussian regime, it is possible to take these correlations into account in a simple way, defining our likelihood as where V −1 is the inverse of the n × n covariance matrix, with V ij = cov[μ i ,μ j ] (leading to V ii = σ 2 i ). However, the off-diagonal elements of this matrix are not given by the experimental collaborations and are very difficult to estimate from outside the collaboration. This remarkably simple and compact expression for the likelihood (a n × n matrix) is only valid under the Gaussian approximation; beyond that the expression and the communication of the likelihood become more complicated.
An alternative way for constraining new physics from the experimental results is to consider results in which the fundamental production and decay modes are unfolded from experimental categories. These so-called "signal strengths in the theory plane" are defined as where as before X labels the production mode and Y the decay mode of the Higgs boson. These quantities can be estimated from a fit to the results in several event categories; as the eff X,Y will differ from measurement to measurement, complementary information on the various (X, Y ) couples can be obtained and break the degeneracies. The resulting signal strengths are directly comparable to the predictions in a given new physics model. They have first been used in phenomenological studies in Refs. [39,40].
It has become a common practice of the ATLAS and CMS collaborations to present such results in 2-dimensional likelihood planes for every decay mode. In that case, the five production modes of the SM Higgs boson are usually combined to form just two effective X modes, VBF + VH and ggH + ttH. The likelihood is then shown in the (µ(ggH + ttH, Y ), µ(VBF + VH, Y )) plane. The ATLAS results in this 2-dimensional plane for Y ∈ (γγ, ZZ * , W W * , τ τ ), as given in Ref. [4], are shown on the right panel of Fig. 1 (the corresponding CMS results can be found in Fig. 5 of Ref. [5]). The solid and dashed contours delinate the 68% and 95% CL allowed regions, respectively. As the unfolding of the individual measurements is done by the experimental collaborations themselves, all correlations between systematic uncertainties (both experimental and theoretical) are taken into account for a given decay mode Y , and are encompassed in the correlation between µ(ggH + ttH, Y ) and µ(VBF + VH, Y ). (Other 2-dimensional planes can be relevant, depending on the sensitivity of the searches.) This is a very significant improvement over the naive combination of categories of Eq. (4), in which all measurements are assumed to be independent. Moreover, in this approach no approximation needs to be made because of missing information on the signal efficiencies or signal strengths category per category. For these reasons, we use the results in terms of signal strengths in the theory plane as the primary experimental input in Lilith.
A remark is in order regarding the grouping of the five production modes into just two. First of all, grouping together VBF, WH and ZH is unproblematic for testing the vast majority of the new physics models because custodial symmetry requires that the couplings of the Higgs to W and Z bosons scale in the same way. The combination of the ggH and ttH production modes is more problematic at first sight. While gluon fusion is dominated by the top-quark contribution in the SM, this can be modified drastically if BSM colored particles are present. However, for all decay modes except H → bb (where gluon fusion-initiated production of the Higgs is not accessible) the ttH production mode is currently constrained with much poorer precision than ggH because of its small cross section. Therefore, with the current data it is justified to take µ(ggH + ttH, Y ) = µ(ggH, Y ) for all channels except H → bb, and µ(ggH + ttH, bb) = µ(ttH, bb).

Statistical procedure
We use signal strengths for pure production and decay modes as basic ingredients for the construction of the Higgs likelihood in Lilith. However, the full likelihood in the µ(X, Y ) basis is not accessible as only 1-and 2-dimensional (1D and 2D) results are provided by the experimental collaborations; some of the correlations should therefore be ignored. In Signal strength results from the CMS collaboration: 1D likelihood for VH, H → bb (red curve) [41] (left), and temperature plot in the plane (µ(ggH + ttH, γγ), µ(VBF + VH, γγ)) [34] (right). the currently available 1D and 2D results, the full likelihood is provided in some cases in addition to contours of constant likelihood. This is extremely helpful since the transmission of the result between the collaboration and the reader does not cause any loss of information. Two examples from the CMS collaboration are given in Fig. 2. The 1D likelihood as a function of µ(VH, bb) [41] is shown on the left panel. 3 On the right panel, the full likelihood in the 2D plane (µ(ggH + ttH, γγ), µ(VBF + VH, γγ)) [34] is shown as a "temperature plot". Moreover, likelihood grids have been provided by ATLAS in numerical format in the 2D plane (µ(ggH + ttH, Y ), µ(VBF + VH, Y )) for Y ∈ (γγ, ZZ * , W W * ) [43][44][45].
Any result given in terms of signal strengths can be used in Lilith. Whenever available, we take into account the full likelihood information. The provision of numerical grids for the di-boson final states by the ATLAS collaboration was an important step forward in the communication of the likelihood. Unfortunately, they were derived with previous versions of the analyses, and the same information has not (yet) been given for the corresponding final Run I results [36,46,47]. Moreover, in the CMS H → γγ result shown on the right panel of Fig. 2, the Higgs boson mass has been profiled over instead of being fixed to a given value, making the interpretation of the result very difficult. Limitations of current way of presenting signal strength results, as well as possible improvements, will be discussed in Section 7.
If only contours of constant likelihood (the 68% CL interval in 1D, 68% and 95% CL contours in 2D) are present, assumptions about the shape of the likelihood have to be made in order to reconstruct it in the full plane. The 1D case was already discussed above, and resulted in the likelihood of Eq. (3). In the 2D case, a natural choice is to use a bivariate normal (Gaussian) distribution. For two (combination of) production and decay processes (X, Y )

Analysis
Type Reference ATLAS H → γγ 2D contour [46] H → ZZ * 2D contour [36] H → W W * 2D contour [47] H → τ τ 2D contour [37] VH, H → bb 2D contour [48] ZH, H → invisible full 1D [49] ttH, H → bb 1D interval [50] ttH, H → γγ 1D interval [46] CMS H → γγ, ZZ * , W W * , bb, τ τ 2D contours [5] ttH, H → γγ, τ τ 1D interval [5] ttH, H → leptons 1D interval [35] ZH + VBF, H → invisible full 1D [51] CDF & D0 VH, H → bb 1D interval [52]  and (X , Y ), we obtain the following likelihood: where µ = µ(X, Y ) µ(X , Y ) , and C −1 = a b b c is the inverse of the covariance matrix. Under the bivariate normal approximation, the 68% and 95% CL contours (which are iso-contours of −2 log L) are ellipses. The information on a single contour suffices to reconstruct the likelihood in the full plane: the parameters a, b and c, as well asμ(X, Y ) andμ(X , Y ), can be fitted from points sitting on the 68% CL or 95% CL contours as they have known values of −2 log L (2.30 and 5.99, respectively). In the following, unless stated otherwise, we choose to reconstruct the full likelihood from a fit to the 68% CL contour provided by the experimental collaboration. However, having more than one contour of constant likelihood is very useful for checking the validity of this approximation. This will be presented in Section 5 for the experimental results included in the database of Lilith. Finally, note that generalization of the previous equations is trivial should higher-dimensional signal strength measurements be published by the experimental collaborations.
A database of up-to-date experimental results is shipped with Lilith, along with recommended sets of results to use for computing the likelihood (in the form of list files; all technical details will be given in Section 4). The default set of results, latest.list, includes the latest measurements from the LHC. Its content as of February 2015 is displayed in Table 1.
All considered 2D results are in the plane (µ(ggH, Y ), µ(VBF + VH, Y )) except for Y = γγ in ATLAS, where only VBF is considered instead of VBF+VH, and for Y = bb in CMS, which is given in the plane (µ(ttH, bb), µ(VH, bb)). The CMS 2D results are taken from the combination of Ref. [5], but correspond to the results from Refs. [34,35,[53][54][55]. We also take into account all available searches on production in association with a top-quark pair, as well as searches for invisible decays of the Higgs boson from both ATLAS and CMS. Note the presence of the CDF and D0 combined result for VH, H → bb [52]; only in this channel the precision of the Tevatron result is comparable with the one of the LHC at Run I.
The final Higgs likelihood is the product of the individual (1-or 2-dimensional) likelihoods. Validation of the Higgs likelihood used in Lilith against official LHC results will be presented in Section 5.

Parametrization of new physics
In order to assess the compatibility of a new physics hypothesis with the LHC measurements presented in the previous section, one needs to compute the expected signal strengths µ(X, Y ) (see Eq. (6)) for the relevant production mechanisms X and decay modes Y . This can be achieved in a direct way from σ(X), σ SM (X), B(H → Y ), and B SM (H → Y ), but is often found to be impracticable. Indeed, in order to have well-defined signal strengths (for which µ = 1 corresponds to the SM prediction) one should take the same prescription for computing cross sections and branching fractions in the SM and in the considered new physics scenario [3]. Concretely, one needs to consider the same order in perturbation theory, the same set of parton density functions, etc.
In most new physics scenarios only leading order (LO) computations are available. Thus, all available next-to-LO (NLO) corrections to the SM predictions should be ignored. While this leads to properly defined signal strengths, σ NLO (X)/σ SM NLO (X) will typically differ from σ LO (X)/σ SM LO (X) (and similarly for branching ratios) as soon as one deviates from the SM prediction. This is because the relative contributions of SM particles to the process will be affected by the NLO corrections. For instance, higher-order corrections to the gluon fusion process will change the relative contribution of the top and bottom quark loops. Therefore, considering LO or NLO cross sections will yield different µ(ggH, Y ) if new physics affects the couplings of the Higgs to top and bottom quarks in a different way.
These two problems come from the parametrization of new physics effects from cross sections and branching ratios. As we will see, they can be alleviated if new physics is parametrized instead using reduced couplings.

Scaling factors and reduced couplings
The general signal strength expression given in Eq. (2) can be rewritten as where Γ SM H is the total decay width of the SM Higgs boson, and the cross section (partial width) for each process X (Y ) is scaled with a factor C 2 accounts for the scaling of the total width of the Higgs boson. (We assume that the narrow-width approximation also holds in the new physics scenarios.) Furthermore, we consider the following effective Lagrangian, where C W,Z and C t,b,c,τ are bosonic and fermionic reduced couplings, respectively. Light fermions are not taken into account as their phenomenological impact on the SM Higgs sector is negligeable. In the limit where all reduced couplings go to 1, the SM case is recovered. At leading order in perturbation theory, the scaling factors C X and C Y from Eq. (8) can be directly identified with the reduced couplings C i from Eq. (9) for processes involving just one coupling to the Higgs boson. We obtain where f = b, c, τ and V = W, Z.
For the remaining main processes (ggH and VBF production, decay into gg, γγ and Zγ), there is no direct identification unless the Higgs couplings to all involved SM particles scale in the same way. In the general case, the C X and C Y for these processes will be given by a combination of reduced couplings C i , weighted according to the contribution of the particle i to the process. For the production mechanisms, we have where the σ SM ij are the different contributions to the cross section in the SM. For i = j it corresponds to the cross section from the particle i alone, while i = j comes from the interference between the particles i and j. (We only consider either the term σ SM ij or σ SM ji in the sum, not both, to avoid double counting the interference terms.) Similarly, the reduced couplings for the gg, γγ, and Zγ loop-induced decay modes are computed as where the Γ SM ij are the SM partial widths of the process under consideration. In all cases, all relevant SM contributions have been taken into account. Note that the relative sign of the reduced couplings will affect the interference terms, as they are proportional to C i C j .
At LO, the various σ SM ij and Γ SM ij can be obtained from tree-level amplitudes (for VBF) or from the 1-loop amplitudes (for gg → H and H → gg, γγ, Zγ). 5 It would however be desirable to take into account the NLO corrections to the Higgs cross sections and partial widths as they modify the relations C X,Y (C i ). This can be achieved in a simple way as long as higher-order corrections only rescale the σ SM ij and Γ SM ij that are already existing in Eqs. (11)-(12), and do not induce new couplings to the Higgs boson. This is the case for the QCD corrections, but not for the electroweak corrections. Thus, as will be explained in Section 4.6.2, the QCD corrections for all five processes of Eqs. (11)-(12) will be included in Lilith.
One last remark is in order. The signal strength framework requires that the signal in all searches be a sum of processes that exist for the SM Higgs boson. However, new production or decay modes may exist without spoiling the signal strength interpretation as long as they do not yield sizable contribution in the current Higgs searches. Two particularly interesting cases are Higgs boson decays into undetected particles, or into invisible particles. In the first case, this new decay is simply missed by current searches (as would, e.g., be the case for the decay of the Higgs into light quarks and gluons), while in the second case this new decay mode gives rise to missing energy in the detector. As was shown in Section 2.3, invisible decays of the Higgs boson are constrained by current searches which are taken into account in Lilith. In both cases of undetected and invisible decays, the width of the Higgs boson becomes larger and modifies the signal strength predictions of Eq. (8) as In Lilith, arbitrary invisible and/or undetected decays can be specified, as will be presented in Section 4.6.

CP-violating admixtures
We also consider the case where the observed Higgs boson is a mixture of CP-even and CP-odd states [56,57]. The Higgs coupling to vector bosons has the general form where as above C V measures the departure from the SM: C V = 1 for a pure scalar H + (CP-even) state with SM-like couplings and C V = 0 for a pure pseudoscalar H − (CP-odd) state.
In the fermion sector, we find the general vector and axial-vector structure of the Higgs coupling to fermions. Concretely, we have where in the SM one has Re(C f ) = 1 and Im(C f ) = 0, while a purely CP-odd Higgs would have Re(C f ) = 0 and Im( to a very good approximation [58]. This is what is implement in Lilith. Effects of CP mixing will mainly show up at loop level, in particular in the gg → H and H → γγ rates. A test of the CP properties of the observed Higgs from a global fit to the signal strengths was presented in [9,59]. Following Ref. [9], at leading order the Higgs rates normalized to the SM expectations can be written as with A + 1 [m W ] −8.32 for m H = 125 GeV. For convenience, the contribution from the other quarks has been omitted in the above equations but is taken into account in Lilith. In the case of ttH production, the approximation that we used above for the other fermions does not hold since m t > m H . Instead, the cross section scales as σ(ttH . Following Ref. [60], a factor σ SM (ttH − )/σ SM (ttH + ) ≈ 1/3 is considered in Lilith. Details on how to specify real and imaginary parts for the couplings are given in Section 4.6.2. More precise measurements at Run II of the LHC will ultimately call for an implementation of CP admixture that includes NLO effects in Lilith.

Getting started
Lilith is a library written in Python for constraining model of new physics against the LHC results. The code is distributed under the GNU General Public License v3.0. The latest version of Lilith and of the database of experimental results (as of February 2015, Lilith 1.1 and database version 15.02) as well as all necessary information can be found at http://lpsc.in2p3.fr/projects-th/lilith The archive of Lilith can be unpacked in any directory. It contains a root directory called Lilith-1.1/ where the following directories can be found: • lilith/: the Python package itself. The Lilith application programming interface (API) will be presented in Section 4.2. It also contains the Python/C API that will be presented in Section 4.4.
• data/: contains the database of experimental results in XML format, as well as *.list text files for the recommended sets of results. Details are given in Section 4.5.
• userinput/: where parametrizations of new physics models, in the XML format described in Section 4.6, can be stored. Some basic user input files that include extensive comments are provided with the Lilith distribution.
• examples/: concrete examples on how to use Lilith for constraining new physics. Two of them will be presented in detail in Section 6; an example for using Lilith in C and C++/ROOT programs will be discussed in Section 4.4.
• results/: empty folder where results from Lilith can be stored.
The folder Lilith-1.1/ moreover contains run lilith.py, the command-line interface (CLI) of Lilith that will be presented in Section 4.3, as well as general information, information on the license, and a changelog in the files README, COPYING, and changelog, respectively.
Lilith requires Python 2.6 [61] or more recent, but not the 3.X series. The standard Python scientific libraries, SciPy and NumPy [62], should furthermore be installed. We require SciPy 0.9.0 or more recent, and NumPy 1.6.1 or more recent. Python, SciPy and NumPy are available for the major platforms, including GNU/Linux, Mac OS X, and Microsoft Windows. If two or more updates of the experimental database are provided the same month, from the second release onwards versions will be numbered YY.MM.n (with n starting from 1). The version number of the database can be found in data/version (also accessible via the readdbversion() method of the API, see next section, and printed on the screen when using the CLI). When using the recommended sets of experimental results of Lilith in a publication, the version number of the database must also be cited.

The Lilith API
Lilith provides an API from which all tasks (reading the user and the experimental input, compute the likelihood, print the results in a file, etc.) can be performed, using the methods described below. This is the recommended way of using Lilith. In order to be used in any Python code (or in an interactive session of Python), the package of Lilith, called lilith, first needs to be imported. However, Python needs to know the location of the lilith package. This can be achieved in at least three ways: before importing lilith. Note that the path can also be relative.
The Lilith library can then be imported to the current script by typing import lilith or from lilith import *. In the first case, all classes, methods and attributes defined in the code will be in the namespace lilith, in the second case they will be in the global namespace. We now present all methods and attributes of the API of Lilith. Lilith.readuserinput (userinput) Read the string in XML format given as argument and fill the attribute couplings (if the user input is given in terms of reduced couplings) or user mu and user mu tot (if the user input is given in terms of signal strengths). User input formats are presented in Section 4.6.

Lilith.computecouplings()
Compute from user couplings the following reduced couplings if not already present in user couplings: C VBF , C ggH , C gg , C γγ , and C Zγ .
Lilith.computemufromreducedcouplings() Compute the signal strengths (stored in user mu and user mu tot) from the reduced couplings in user couplings.
Lilith.compute user mu tot() Add up the signals from all Higgs bosons contributing to the signal in user mu; store the result in user mu tot.

Lilith.readexpinput(filepath=default exp list)
Read the experimental input specified in a list file and store the results in exp mu and exp ndf. By default, the list file is data/latest.list. The formats of the experimental results are presented in Section 4.5.

Lilith.readdbversion()
Read the version of the database of experimental results from the file data/version, and store the information in dbversion.
Lilith.compute exp ndf() Compute the number of measurements from exp mu and store the information in exp ndf.
Lilith.computelikelihood(userinput=None, exp filepath=None, userfilepath=None) Evaluate the likelihood function from signal strengths derived from the user input (user mu tot) and the experimental results (exp mu) and store the results in the attribute results. If the arguments userinput and userfilepath are not specified, user mu tot will be assumed to have been filled already. Else, all information will be read from the XML input given in userinput, or from the file located at userfilepath, and user mu tot will be computed. If the exp filepath argument is not specified, the experimental results will be read from the default list file unless exp mu is already filled. Else, experimental results from exp filepath will be read before computing the likelihood.

Lilith.writecouplings(filepath)
Write reduced couplings from the attribute couplings in a file located at filepath in the XML format specified in Section 4.6.2.
Lilith.writesignalstrengths(filepath, tot=False) Write signal strengths from the attribute user mu (if tot=False) or user mu tot (if tot=True) at filepath in the XML format specified in Section 4.6.1.

Lilith.writeresults(filepath, slha=False)
Write the content of the attribute results at the location filepath in the XML format (if slha=False) or the SUSY Les Houches Accord (SLHA)-like format [63] specified in Section 4.7.
Note also that the version of Lilith is stored in the lilith. version variable.
A minimal example of use of the API is as follows: lcal.readuserinputfile( userinput/example_mu.xml ) 5 lcal.computelikelihood() 6 print -2log(likelihood) = , lcal.l The first two lines import the Lilith library into the global namespace and initialize the computations. They are equivalent to The three following lines successively read the experimental input, read the user input from the file userinput/example mu.xml, and compute the likelihood. Alternatively, they could be replaced with a single line, Finally, the value of −2 log L is printed on the screen on the last line.
In the example above, any error (corresponding to an exception in Python) will interrupt the code. This may not be the desired behavior. In particular, if several user inputs are successively given to Lilith (as in the case of a scan of a parameter space), it may be preferable to store the error and move on to the next user input instead of stopping the execution of the code. In Python, the handling of errors can be achieved with try .. except blocks. We provide below a simple example. Here, any error raised by Lilith (of type lilith.LilithError or derived from it) will be catched, in which case the error message will be printed on the screen and the script will continue normally. It is of course also possible to store the error message in a file, or simply replace the last line with the pass statement in order to ignore all errors raised by Lilith and continue the execution of the script. For the definition of all exceptions used in Lilith (that all derive from lilith.LilithError), see lilith/errors.py. This makes it possible to treat each type of error in a different way.

Command-line interface
A command-line interface or CLI is also shipped with Lilith for a more basic usage of the tool. It corresponds to the file run lilith.py located in the directory Lilith-1. where arguments in parentheses are optional.
The first argument, user input file, is the path to the user input file in the XML format described in Section 4.6. New physics can be parametrized in terms of reduced couplings (see Section 3) or directly in terms of signal strengths. Examples are shipped with Lilith in the directory userinput/. The second argument, experimental input file, is the path to the list of experimental results to be used for the construction of the likelihood. If not given, the latest LHC results will be used (data/latest.list; its content is given in Table 1). It is the recommended list of experimental results to be used for performing a global fit. All details about the experimental input will be given in Section 4.5.
If no option is given, basic information as well as the value of the likelihood and the number of measurements is printed on the screen. A number of options are provided to control the information printed on the screen and to print the results of Lilith in output files. They are listed in Table 2.
The option --couplings / -c only works in the reduced couplings mode. In addition to the reduced couplings already present in the input file, it prints scaling factors computed from the input (i.e. C γγ , C Zγ , C gg , C ggH , and C VBF , see Section 3.1). The option --mu / -m prints the complete list of signal strengths in a file in XML format. More specifically, all signal strengths µ(X, Y ) with X ∈ (ggH, VBF, WH, ZH, ttH) and Y ∈ (γγ, ZZ * , W W * , bb, τ τ , cc, Zγ, gg, invisible) are printed. Finally, the option --results / -r prints the value of −2 log L and the number of measurements in a file. If the extension of the filename is .slha (case-insensitive), a file in SLHA-like format is created. Otherwise, an XML file is created, with extra information on the individual contributions to −2 log L from all the experimental results used in the calculation. More details on the structure and content of the output files are given in Section 4.7.
Initialization of Lilith and reading of the database of experimental results is done at each execution of run lilith.py. When computing the Higgs constraints in the context of the scan of a model, this only needs to be done once. Therefore, successive calls to the CLI will be much slower than direct calls to the methods of the API (see previous section), even more  so as internal information is stored to optimize successive computations of the results when using the API. Whenever performance can be an issue, we highly recommend using the API instead of the CLI.

Interface to C and C++/ROOT
Above, we have presented how to use Lilith from the API written in Python and from the CLI. In order to use Lilith within a C, C++, or ROOT code, a first possibility is to call the CLI. However, this will suffer from the performance issues explained at the end of the previous section. Fortunately, Python provides a Python/C API [64] from which each method of the API of Lilith can be called, and each attribute of Lilith can be manipulated, in both C and C++. However, direct use of the functions of the Python/C API can be quite tedious.
For that reason, we have also written a C API for Lilith, consisting in a series of functions using the Python/C API for performing all usual tasks, following closely the methods of the API presented in Section 4.2.
The C API for Lilith is contained in the directory lilith/c-api. Its use requires the libraries and header files needed for Python development. They can be installed from most package managers under the name python-dev or python-devel, depending on the platform. More information on how to install these libraries and header, platform by platform, can be found at [65].
We now present all the functions of the C API: initialize lilith(char* experimental input) Import lilith, instanciate the class Lilith, and read the experimental input file located at experimental input. The function returns the instance object. If experimental input is an empty string(""), the default experimental input file data/latest.list is used.
lilith readuserinput(PyObject* lilithcalc, char* XMLinputstring) Read the user input XML string XMLinputstring and store the information in the object lilithcalc.
lilith readuserinput fromfile(PyObject* lilithcalc, char* XMLinputpath) Read the user input XML file located at XMLinputpath and store the information in the object lilithcalc.

Experimental input
We have seen that the evaluation of the likelihood in Lilith requires the input of a list of experimental results to be considered. It corresponds to a simple text file with a .list extension listing the paths to experimental result files in XML format (each containing a single 1D or 2D signal strength result). Lilith is shipped with the latest LHC Higgs results (plus a Tevatron result), see Table 1 in Section 2.3, in the form of XML files present in subdirectories of data/. Moreover, several lists of experimental results are provided in data/, with latest.list being the default list file. This is the one recommended for a global fit to the LHC+Tevatron Higgs data.
The user can also create his/her own list file in the data/ directory. For instance, in order to put constraints on new physics using only the latest di-boson results from ATLAS [36,46,47], one can create a file list that contains # ATLAS di-boson analyses ATLAS/Run1/HIGG-2013-08_ggH-VBF_gammagamma_n68.xml ATLAS/Run1/HIGG-2013-21_ggH-VVH_ZZ_n68.xml ATLAS/Run1/HIGG-2013-13_ggH-VVH_WW_n68.xml The first line, starting with a #, is a comment and is not read by Lilith. The three following lines indicate the paths to the XML files to be considered. As can be read in the paths, these are published results from the ATLAS collaboration based on Run I data. The conventional, though not mandatory, naming scheme for the files is as follows. The identifier of the analysis (HIGG-2013-XX in this case) comes first, followed by the 2D plane in which results are given; finally, n68 indicates that the likelihood has been reconstructed from the contour at 68% CL under the Gaussian approximation. Note that no consistency check is done by Lilith.
When creating a new list file, the user should make sure that there is no overlapping between experimental results (e.g., that the results of two versions of the same analysis, based on overlapping event sets, are not used at the same time).

XML format
Every single experimental result (1D or 2D) is stored in a different XML file. In this way, modifying and updating the database is an easy process. We now present the format of the experimental results files in Lilith.
The root tag of each file is <expmu>. It has two mandatory attributes, dim and type, that specify the type of signal strength result (1D interval, full 1D, 2D contour, or full 2D, see Section 2.3). Possible values for the attributes are given in Table 3. In addition, the <expmu> tag has two optional attributes: prod and decay. They can be given a value listed in Table 3 if the analysis under consideration is only sensitive to one production mode (e.g., ttH) or to one decay mode (e.g., γγ) of the Higgs boson. In the general case, the prod and decay attributes can be skipped. Indeed, all relevant efficiencies eff X,Y (see Section 2.1) can be specified in <eff> tags.
Before turning to the syntax case by case, we comment on the optional information that can be provided in the experimental files. The following tags can be given: • <experiment>, referring to the experiment that has produced the current result.
• <source>, that contains the name of the analysis, and has attribute type that contains the status of the analysis (published or preliminary).
• <CL>: when the likelihood has been extrapolated from a 2D contour, this can be used to indicate the CL of the contour thas has been used to extract the covariance matrix (usually the 68% or 95% CL contour).
• <mass> contains the Higgs boson mass considered in the results. If not given, m H = 125 GeV is assumed.
Note that, if prod="VH" or prod="VVH" is given as attribute to the <expmu> tag or to an <eff> tag, the relative contributions of WH and ZH (for VH) and of WH, ZH and VBF (for VVH) will be computed internally assuming an inclusive search, i.e., for VVH, where the cross sections are evaluated at the Higgs mass given in the <mass> tag using the LHC Higgs Cross Section Working Group (HXSWG) results for the 8 TeV LHC [3].

full 1D
We consider as an example the H → bb CMS search [41]. The 1D profile likelihood as a function of µ(VH, bb) is provided in Ref. [ The digitized likelihood information is stored in the tag <grid>. The first column is the signal strength (whose nature is determined by the <eff> tags) while the second column is the value of −2 log L.

2D contour
We consider as an example the CMS search H → bb [42]. The 68% and 95% CL contours are provided in the plane (µ(WH, bb), µ(ZH, bb)). As was explained in Section 2.3, we start by fitting the 68% CL contour assuming that the likelihood follows a bivariate normal distribution, and we extract the experimental best-fit point and the inverse of <bestfit> <x>1.123</x> <y>0.997</y> </bestfit> <param> <a>1.393</a> <b>0.190</b> <c>2.217</c> </param> The tag <bestfit> specifies the location of the best-fit point in the (x,y) plane. The tag <param> contains the sub-tags <a>, <b>, and <c>, that parametrize the inverse of the covariance matrix in the (x,y) plane.

full 2D
We consider as an example the H → γγ ATLAS search [66] for which the full likelihood information is given in the plane (µ(ggH + ttH, γγ, VBF + VH, γγ)) [43].  The tag <grid> contains the grid provided by the experimental collaboration. The first and second columns are defined by the axis="x" and axis="y" attributes of the <eff> tag, respectively. The third column is the value of −2 log L.

User model input
The user model input, parametrizing the new physics model under consideration, can be given either in terms of signal strengths µ(X, Y ) directly (defined as in Eq. (6)), or in terms of reduced couplings and scale factors (see Section 3). In the latter case, scale factors that might be missing in the input are computed, and signal strengths are derived from the scale factors.
The user model input has XML syntax and can be provided as a string or in the form of a file (see the methods readuserinput() and readuserinputfile() in Section 4.2). In this section, we present the format that is used in Lilith.

XML format for signal strengths
In the signal strengths mode, the basic inputs are the signal strengths defined as in Eq. (6). An example of XML input file for the signal strengths mode is now presented.
• The <signalstrengths> tag indicates that the user input is given in terms of signal strengths.
-The decay attribute can be gammagamma, Zgamma, WW, ZZ, bb, cc, tautau. As for the prod attribute, multi-particle labels have been defined and are listed in Table 4.
Note that every <mu> tag can be omitted; in such a case the SM value will be assumed. A warning will furthermore be issued in case of missing µ(X, Y ) for Y ∈ (γγ, ZZ * , W W * , bb, τ τ ) (after resolving multi-particle labels).
• Finally, there is the possibility to specify an invisible branching ratio in the <redxsBR> tag. This is defined as As the branching fraction of the SM Higgs boson into invisible particles is very small (B SM (H → 4ν) = 0.11% at m H = 125 GeV [3]) and cannot be probed at the LHC, one usually does not express the results of invisible Higgs searches in terms of signal strengths. Invisible decays of the Higgs boson are currently constrained in association with two jets from VBF, and in association with a Z boson from ZH production, see Table 1.
Note that the signal strengths for several Higgs states contributing to the signal can be defined by specifying an arbitrary number of <signalstrengths> .. </signalstrenths> tags in the input. After reading the input, the signal strengths from each individual state contributing to the signal will be stored in the attribute user mu, and the sum of the signal from the different particules (signal strength per signal strength) will be stored in the attribute user mu tot (for more details, see Section 4.2). We neglect possible interferences between the different states. It can be useful to provide an identifier for each particle. This can be achieved with a part attribute to the <signalstrengths> tag. An example of user input in terms of signal strengths is stored in userinput/example mu.xml for the case of a single Higgs boson contributing to the signal, and in userinput/example mu multiH.xml for the case of two or more particles.

XML format for reduced couplings
New physics can be parametrized in terms of scaling factors that can be identified as (or derived from) reduced couplings, as was presented in Section 3. In this section we present the user input in terms of reduced couplings. Before turning to the format of the user input, that also has XML syntax, we comment on the computation of couplings and of signal strengths. First of all, as we have seen in Section 3.1, predictions for the Higgs boson can be obtained from the reduced couplings C W , C Z , C t , C b , C c , and C τ appearing in Eq. (9). Scaling factors for VBF production and loop-induced processes are function of the C i and can be expressed in as Eq. (11)- (12). In the following, we will consider two possible cases: that these scaling factors are obtained from leading-order calculations (i.e. tree-level results for VBF and one-loop analytical expressions for other processes), or including NLO QCD corrections. The former case will be denoted as LO, the latter one as BEST-QCD. We comment on the computations currently implemented in Lilith:

VBF
The contribution from the W boson, the one from the Z boson, and the interference between them have been obtained from VBFNLO-2.6.3 [67] for Higgs masses in the [123, 128] GeV range with (for BEST-QCD mode) and without (for LO mode) NLO QCD corrections at the LHC 8 TeV, using the MSTW2008 parton distribution functions [68]. The results for σ SM W W (VBF), σ SM ZZ (VBF) and σ SM W Z (VBF) as a function of the Higgs mass were stored in text files shipped with Lilith and read internally when using computereducedcouplings() (see Section 4.2).

ggH
The contributions from the three heaviest quarks (t, b, c) to the SM cross section are taken into account. In the LO mode, we use analytical expressions [58]. In the BEST-QCD mode, those have been generated in the [123, 128] GeV range with HIGLU [69] at the LHC 8 TeV with the MSTW2008 parton distribution functions.

H → gg, γγ, Zγ
The relevant SM partial widths of these processes (taking into account particles listed in Eq. (12)) are obtained from analytical expressions [58] in the LO mode. In the BEST-QCD However, the effective Lagrangian defined in Eq. (9) does not exhaust the possibilities for new physics affecting the properties of the Higgs processes. One particularlity interesting case is that BSM particles enter the loop-induced processes, such as gg → H and H → γγ.
An explicit example will be given in Section 6.2. To account for these cases, we allow direct definition of scaling factors for the four main loop-induced processes (gg → H and H → gg, γγ, Zγ), i.e. direct definition of C ggH and C gg,γγ,Zγ . If some or all of the scaling factors are missing from the input, they will be computed internally using Eq. (11)- (12), i.e. assuming that only SM particles are involved. Finally, note that we use the SM branching ratios provided by the LHC HXSWG [3], at the Higgs mass given in the user input, when computing the signal strengths (see Eq. (8)).
The user input file for the Lilith reduced couplings mode has the following structure.
• The <reducedcouplings> tag is specific to the reduced couplings mode. This is where the reduced couplings are specified. The correspondence between the XML notation and Eq. (9) is given in Table 5. Note the possibility to define common couplings for the up-type fermions, down-type fermions, all fermions, and electroweak gauge bosons.
• The tag <mass> defines the Higgs boson mass at which the likelihood should be computed. The allowed range is [123,128] GeV. This affects the computation of the SM branching ratios and partial cross sections and widths as explained above. If it is not given, a Higgs mass of 125 GeV is assumed.
• Regarding the effective coupling to a pair of gluons, NLO corrections affect gluon fusion (ggH) and the decay into two gluons (H → gg) in a different way. Therefore, scaling factors C ggH and C gg can be specified separately as <C to="gg" for="prod">1.0</C> <C to="gg" for="decay">1.05</C> If for="all" is specified, the same coupling is assigned to the production and decay modes. This is the default behavior if the for attribute is missing.
• CP violation was presented in Section 3.2. In the LO mode, the fermionic couplings C t , C b , C c , C τ can be given a real and an imaginary component. For the top quark for instance, this can be specified as <C to="tt" part="re">0.8</C> <C to="tt" part="im">0.2</C> If part="re|im" is not specified, the coupling is assumed to be purely real. In the BEST-QCD mode, only the real part of the coupling is taken into account.
• The <precision> tag contains either BEST-QCD or LO. If not specified, or wrongly spelled, the BEST-QCD mode is the default mode.
• The <extraBR> tag contains the declaration of the invisible or undetected branching ratios (see Section 3).
As in the case of input in terms of signal strengths, several tags <reducedcouplings> .. </reducedcouplings> can be defined, corresponding to the case where several Higgs states contribute to the observed signal around 125 GeV. These particles can also be given a name with a part attribute to the <reducedcouplings> tag. An example of user input in terms of reduced couplings is stored in userinput/example couplings.xml for the case of a single Higgs boson contributing to the signal, and in userinput/example couplings multiH.xml for the case of two or more Higgs states.

Output
When using the API, all relevant information can be accessed from the public attributes of the class Lilith presented in Section 4.2. The user can manipulate and store this information in any way he/she wants. However, having standardized output formats is important in order to interface Lilith with other programs. We provide three output methods: writecouplings(), writesignalstrengths(), and writeresults() (for details on how to call these methods, see Section 4.2). The first two methods write the reduced coupling information (if available) and the signal strength information, respectively, in a file. The formats specified in Section 4.6.1 and 4.6.2 above are used, such that these output files can also be used as input to a subsequent call to Lilith. In the command-line interface, these two methods can be called with the options -c or --couplings, and -m or --mu, respectively (for more details, see Section 4.3).
The last method, writeresults(), is used to store the results after evaluation of the likelihood. Depending on the second argument, the output will be written in XML format (if slha=False) or in SLHA-like format (if slha=True). By default the output file is in XML format for consistency with what is used otherwise in Lilith; an SLHA-like output was added as it is widely used in BSM phenomenology. This method can be called with the option -r or --results in the command-line interface.
We start with the description of the XML format. Its root tag is <lilithresults>. The results from each analysis are given within an <analysis> tag having two attributes: experiment and source, corresponding to the value of the tags <experiment> and <source> read from the experimental input, if present; otherwise the attribute is given an empty value "". Each <analysis> tag contains a <l> tag, whose value is −2 log L from this experimental result alone, and an <expmu> tag that follows the syntax of the experimental input (see Section 4.5) except that it only provides information on the efficiencies (in <eff> tags).

Validation
Having explained how to use Lilith in the previous section, we now turn to the validation of the likelihood derived from the experimental input shipped with the code. We begin by discussing the validity of the bivariate normal distribution as an approximation to the 2D likelihood functions in the signal strength planes (µ(X, Y ), µ(X , Y )). The use of this approximation is necessary whenever only contours of constant likelihood are provided instead of the full information. Several coupling fits from the ATLAS and CMS collaboration are then reproduced. The results from Lilith are compared to the official ones to assess the validity of the likelihood used in Lilith.

Reconstruction of the experimental likelihoods
In a signal strength plane (µ(X, Y ), µ(X , Y )), an approximation to the likelihood function can be obtained assuming that the measurements follow a bivariate normal distribution, as explained in Section 2.3. Using the 68% CL contour provided by the experimental collaboration, we reconstruct the shape of the likelihood and compare the location of the best-fit point as well as the 68% and 95% CL contours with what is provided by ATLAS or CMS.
Two examples are shown in Fig. 3: the reconstruction of the likelihood for the ATLAS W W * [47] and the CMS γγ [34] final states. 6 In both cases we observe an excellent agreement between the reconstructed likelihood and the official result. The 68% CL regions are perfectly reproduced and the reconstructed best-fit points are very close to the experimental ones. The extrapolation towards the 95% CL regions also shows very good agreement. We find equally good agreements with all other decay modes (with the exception of H → ZZ * ), and we conclude that the Gaussian distribution is a very good approximation to the true distribution.
The largest deviations from the normal approximation are expected to occur for final states with low statistics since the counting of the events, that follows the Poisson distribution, has not yet entered the Gaussian regime. In particular, this is the case for the ZZ * channel. In Fig. 4, we show the comparison between the Lilith reconstructed likelihood in the ZZ * final state and the corresponding ATLAS and CMS ones.
As can be seen, the deviation of the ATLAS likelihood from the bivariate normal approximation can be substantial. In the positive region of the plane (the one that is relevant), the    approximation holds well near the best-fit point. However, going away from it the reconstructed shape fails to reproduce the ATLAS 95% CL contour at large µ(VBF + VH, ZZ * ). Due to non-Gaussian effects, the reconstructed best-fit point is quite distant to the experimental one. For the CMS case, the approximation holds to a better approximation. The reconstructed best-fit point is very close to the experimental one and the shape of the 95% CL contour is very well reproduced although a small shift in the µ(ggH + ttH, ZZ * ) direction is observed. As will be argued in Section 7, provision of the full likelihood information would result in a significant improvement of the Higgs likelihood reconstruction.

Comparison to Higgs coupling fits from ATLAS and CMS
In order to validate the approximate Higgs likelihood used in Lilith, we attempt to reproduce coupling fit results from combination notes from ATLAS [4] and CMS [5]. Note that while the CMS combination [5] makes use of the final Run I results, a number of analyses considered in the combination of the ATLAS results given in Ref. [4] have been updated since then. The final, legacy combination of the Higgs measurements from ATLAS at Run I has not yet been released. For this reason, the lists of recommended experimental results data/CMS-HIG-14-009.list and data/latestCMS.list are identical, as of February 2015, while data/ATLAS-CONF-2014-009.list and latestATLAS.list differ.
First, results for two benchmark scenarios proposed by the LHC HXSWG in [71] are presented. In the first scenario, SM-like tree-level couplings are assumed (i.e., all C i = 1 in Eq. (9)) but two scaling factors are introduced: C γ ≡ C γγ (scaling H → γγ), and C g ≡ C ggH = C gg (scaling ggH production and H → gg). In the second benchmark scenario, two reduced couplings are introduced: C V ≡ C W = C Z , for the coupling of the Higgs boson to a pair massive vector bosons, and C F ≡ C t = C b = C c = C τ , a universal coupling to fermions. In this case, the effective coupling to gluons is simply C F , while C γγ is a function of both C V and C F that was obtained taking into account QCD corrections (see Section 4.6.2).
Let us first discuss the results from ATLAS, obtained using the list of experimental results data/ATLAS-CONF-2014-009.list. Results of the two fits are presented in Fig. 5. In both scenarios, very good agreement is observed between the results from ATLAS and the ones obtained with Lilith. Both the reconstructed best-fit point and contours reproduce very well the ATLAS results. The most significant deviation is a slight deformation of the 95% CL region in the (C γ , C g ) plane. The corresponding results for CMS are shown in Fig. 6. CMS results are well reproduced with Lilith, even for the contour at 99.7% CL. Slight shifts of the best-fit points and minor deformations of the contours are observed. The overall agreement is nevertheless very good.
Let us move on to the 3-parameter fit (C W , C Z , C F ). As in the (C V , C F ) benchmark scenario discussed above a universal coupling to fermions is introduced, but instead of a single coupling to vector boson one defines separately the reduced coupling to W bosons, C W , and to Z bosons,    (left) and CMS [5] (right) data and comparison to the official results. The ATLAS fit considers both signs for the Higgs-fermion-fermion coupling and furthermore defines C The results are given for both signs of C F Z .
using the Higgs measurements alone. The 1-dimensional likelihood profile for C W Z is shown in Fig. 7 for both the ATLAS and CMS combination.
Although the ATLAS result is almost perfectly reproduced, a significant discrepancy is observed in the case of CMS for C W Z > 1. This does not come as a surprise: several experimental results were considered in the (µ(ggH + ttH, Y ), µ(VBF + VH, Y )) plane. The breaking of VBF + VH into the individual production modes VBF, WH and ZH (assumed to be inclusive, see Eq. (17)) becomes relevant for C W = C Z . Moreover, ATLAS results make use of the full numerical likelihood grids that were provided in Refs. [43][44][45] while the bivariate normal approximation is used in the case of CMS. Thus, constraints on models in which C W = C Z should be interpreted with care given the information currently available.
Finally, we present the result of a 3-parameter fit (C γ , C g , B invisible ) in terms of the 1D profile likelihood of B invisible in Fig. 8. A very good agreement is observed in ATLAS, and in CMS for moderate values of B invisible . As explained in Section 3.1, the presence of a branching ratio into invisible particles is constrained by direct searches for invisible decays of the Higgs boson, and also by every Higgs search since it modifies the total Higgs width and therefore scales all signal strengths collectively.

Examples of applications
Having validated the Higgs likelihood in Lilith from results obtained by the ATLAS and CMS collaborations, we now turn to deriving constraints on specific new physics scenarios using the latest LHC results (as of February 2015) present in data/latest.list. The Python routines used to obtain these results are available in the folder examples/python and will be described shortly.

Reduced coupling determination
As a first illustration of the use of Lilith, constraints on the benchmark scenario (C V , C F ) introduced in Section 5.2 are derived. The right panels of Figs. 5 and 6 show results on this scenario in the 2D plane (C V , C F ) using only ATLAS or CMS results. Here we combine the latest ATLAS and CMS results and derive 1D profile likelihood constraints on C V and C F . Results are shown in Fig. 9.
The Python routine used to obtain this result is CVCF 1dprofile.py. It can be executed from the Lilith-1.1/ folder with the command line    given range and number of points. 7 Without the option substract min=True, the "absolute" likelihood −2 log L(C V ) would be returned instead.
The parameter range in which ∆(−2 log L(C V )) < 1 (4) defines the 68% (95%) CL intervals of C V . The constraints on C F are derived in the same way, and all results are plotted and stored in results/CVCF 1dprofile.pdf. In this scenario, the best-fit point is obtained for We also provide example on how to derive constraints and produce figures for a 2D parameter space. The left panel of Fig. 10 presents the 2D contraints obtained from a global fit of the (C γ , C g ) model presented above. The corresponding Python routine is CgammaCg 2d.py. It can be executed from the Lilith-1. where the first, second and third columns contain the values of C γ , C g and −2 log L(C γ , C g ), respectively. The 68%, 95%, 99.7% CL regions in the (C γ , C g ) plane then corresponds to ∆(−2 log L(C γ , C g )) < 2. where xi, yi and Z are list of points defining the grid in the C γ and C g directions, and the corresponding ∆(−2 log L(C γ , C g )) value, respectively. The results are displayed and stored in results/CgammaCg 2d.pdf. For completeness, the 2D constraints on the (C V , C F ) benchmark scenario, using the latest LHC measurements, are also presented in the right panel of Fig. 10. They have been derived in the same way.

Higgs constraints on superpartners of the tau lepton
Supersymmetric scalar partners of the tau leptons, known as staus, can have substantial contribution to the H → γγ decay rate if they are light and have a large mixing [75,76]. Constraints on the parameters controlling this new contribution can therefore be obtained from the Higgs precision measurements. Here, we consider the Minimal Supersymmetric Standard Model (MSSM) and assume that the only deviation from a SM-like Higgs behavior comes from the contribution of staus to the loop-induced process H → γγ. More precisely, it is assumed that the supersymmetric partners of the Higgs boson and of the remaining fermions are decoupled, that the second Higgs doublet is phenomenologically irrelevant and that a Higgs mass of 125 GeV can be obtained for any point of the analysis. In this case, the contribution from staus to the H → γγ decay width is parametrized by the two physical masses mτ 1 and mτ 2 (with mτ 1 < mτ 2 ), the mixing angle θτ and the ratio of vacuum expectation values for the two Higgs doublets, tan β. The corresponding amplitude at LO reads [58,77] Mτ where the sum runs over the two stau mass-eigenstates, A H 0 is a loop form factor and g Hτ iτ i is the Higgs-stau-stau coupling.
Note that the SM amplitude M SM Hγγ appears both in the numerator and denominator of Eq. (20) since SM tree-level couplings are assumed.
The corresponding Python code is stau gammagamma.py. It can be executed by typing the following command line to the shell from the Lilith-1.1/ folder: The routine works as follows. Functions returning C γ according to Eq. (20) and −2 log L(C γ ) are defined. Since tan β and mτ 1 are fixed, a 2-dimensional grid scan is then performed over the two remaining parameters: for each couple (mτ 2 , θτ ), the corresponding ∆(−2 log L) is obtained. The 2-dimensional 68%, 95%, 99.7% CL regions in the plane (mτ 2 , θτ ) are obtained with ∆(−2 log L) < 2.3, 5.99, 11.83, respectively.
Note that direct searches from LEP [78] and vacuum metastability condition [79] impose further constraints on this scenario. Moreover, this simplified SUSY scenario could easily be generalized, e.g. by taking into account H →χ 0 1χ 0 1 . Light staus are especially relevant in the case whereχ 0 1 is light in order to have a viable neutralino dark matter candidate in the MSSM (see, e.g., Ref. [80] and references therein).

Prospects for Run II of the LHC
As we discussed in Section 2, approximations necessarily need to be made when combining signal strength results from several categories or several searches, making it necessary to validate the approach. In Section 5, we have shown that we reproduce well the results of coupling fits from ATLAS and CMS (separately). However, it is clear that the situation will change as more statistics will be collected at Run II of the LHC. Indeed, systematic uncertainties will then dominate over statistical uncertainties in the majority of the channels. Missing correlations between systematic uncertainties (both theoretical and experimental) will thus become a more pressing issue. Moreover, more combinations for production and decay of the Higgs boson, (X, Y ), will be determined with a good precision. This will spoil the simple interpretation we have for a number of the results we currently use, in particular for results given in the plane (µ(ggH + ttH, Y ), µ(VBF + VH, Y )). In this section we recall the main limitations when using the information currently provided by the ATLAS and CMS collaborations for constructing a likelihood. We also discuss new ways of presenting the LHC Higgs results in order to be able to construct a good approximation to the Higgs likelihood at Run II of the LHC. This section is partly based on the note "On the presentation of the LHC Higgs Results" [81] that was put forth by a collaboration of theorists and experimentalists with the aim to maximize the impact of the LHC Higgs results and their utility to the whole high-energy physics community.
First of all, in most cases only contours of constant likelihood (at least the 68% CL interval or contour, sometimes contours at 95% CL) are provided by ATLAS and CMS. This makes it necessary to extrapolate the likelihood assuming, most naturally, a Gaussian shape. When using a given contour to extrapolate the likelihood, the validity of this approximation can be tested from a comparison of the position of the best-fit point and from contours provided by the experimental collaboration. This was done in Section 5.1, where we concluded that the reconstruction is generally very good, although in some cases asymmetrical effects are washed out (see, e.g., Fig. 4). However, in all cases it induces an unnecessary source of error. The ATLAS and CMS collaborations initiated some efforts during Run I to provide the full likelihood information in 1D and 2D planes (see the right panel of Fig. 2 and Refs. [43][44][45]). We strongly hope that this will become standard practice during Run II of the LHC, and that the information will systematically be provided in numerical form.
Another issue is the dependence of the results on the assumed Higgs boson mass m H . Currently, we use results given at a fixed Higgs mass. As not all results are provided at the same m H , a slight inconsistency is introduced in the combination of the different results (the assumed Higgs mass varies within a few hundreds of MeV). Official combination notes allow us to get rid of this inconsistency, as all results are therein given at the same Higgs mass. However, the dependence of the experimental results on the Higgs mass can be very important for the high-resolution channels, that target decay of the Higgs boson into charged leptons and photons (such as H → ZZ * → 4 and H → γγ). Thus, it would be highly desirable to have access to mass-dependent likelihood results.
Current results are presented in 1-or 2-dimensional projections, often corresponding to the combination of production modes (in 2D, typically ggH + ttH and VBF + VH). As we discussed above, this becomes a limitation as measurements get more precise, in which case we would like to investigate deviations in all of the five production modes separately. For such reasons, a total breakdown of the signal strength measurements in terms of the five Higgs production modes (ggH, VBF, WH, ZH, ttH) would be a considerable step forward regarding the interpretation of the LHC Higgs results. We would therefore like to advocate the experimental collaborations to provide the likelihood as a function of the Higgs mass and a full set of production modes, that is to say, in the (m H , µ(ggH, Y ), µ(ttH, Y ), µ(VBF, Y ), µ(ZH, Y ), µ(WH, Y )) (21) parameter space for each final state Y . For some final states, all five production modes are certainly not constrained by the experimental searches and only lower dimensional projections of this space would be relevant.
This would solve most of the limitations currently faced, with the notable exception of correlations between the measurements of different decay modes. For instance, theoretical uncertainties on gluon fusion production affect both the γγ and ZZ * final states. Recently, an interesting proposal was made in this direction in Ref. [82]. Provided experimental collaborations publish likelihoods that are not profiled over a set of theoretical nuisance parameters of interest, but instead given for a fixed scenario, it is possible to build a "recoupled" likelihood incorporating these uncertainties at the later stage. This has the advantage of not being restricted to the Gaussian approximation. It would certainly be of great interest if the information in the 2D plane (µ(ggH + ttH, Y ), µ(VBF + VH, Y )), or even better in the possibly 6D plane discussed above, could be given without profiling over the theoretical uncertainties on the Higgs signal. With the method presented in Ref. [82], one could then fully correlate the theoretical uncertainties between the different channels and experiments, and modify these uncertainties compared to what is done in ATLAS and CMS if desired.

Conclusions
Crucial information on the origin of the electroweak symmetry breaking-and, more generally, on the presence of BSM physics around the electroweak scale-can be obtained from the study of the properties of the Higgs boson with mass around 125 GeV at the LHC. The presentation of the experimental results in terms of signal strengths makes it possible to combine all measurements and perform a global fit to the properties of the Higgs boson. The results of such fits can be used to discriminate between models.
However, using all available information from the ATLAS and CMS collaborations to construct a likelihood is a non-trivial task. Indeed, there is a wealth of experimental searches, from which the necessary information often needs to be extracted and put in numerical form. Care should also be taken in order to include all available correlations between systematic uncertainties. To this aim, we provide a new public tool, Lilith. Lilith is a library written in Python, and for which we provide an API as well as a command-line interface and a basic interface to C and C++/ROOT. The experimental results are read from a database in XML format that is shipped with the code and which is easy to modify and extend. Lilith uses as a primary input results in which the fundamental production and decay modes are unfolded from experimental categories.
New physics can be parametrized in terms of reduced couplings, or signal strengths directly, which are given as input to Lilith in XML format. If needed, scaling factors for the loopinduced processes and VBF production are computed taking into account QCD corrections. CP-violating Higgs couplings can also be given as input to Lilith. The likelihood is evaluated from a set of experimental results and given as output; detailed results can moreover be stored in XML or SLHA-like format. For convenience, Lilith is provided with several applications of the code where constraints on effective or explicit models of new physics are derived. They include extensive comments.
The Higgs likelihood of Lilith obtained from the latest measurements at the LHC has been thoroughly validated against ATLAS and CMS results and can be used to constrain new physics. Future measurements at Run II of the LHC will, however, call for new ways of presenting results in order to derive a good approximation to the Higgs likelihood. In particular, further disentanglement of the different production and decay modes will become necessary. Moreover, correlations between systematic uncertainties, and in particular the treatment of theoretical uncertainties, will become a more pressing issue. The structure of the code is such that Lilith can easily be adapted to handle extended signal strength information.