Higgs combine tutorial based on CMSDAS2014HiggsCombExercise

• • First fits to Higgs search results: use what learned before to perform a first fit to the CMS datacards to produce the signal strength plots like these two
    

• ROOT commands and C++ macro will be embedded in light green box, e.g.

  %CODE{"cpp"}% limit->Draw("2*deltaNLL:r"); %ENDCODE%

• Combinations denoted with the skull (a˜ ) are not thought to be possible to measure mostly because of the lack of distinctive features with which to tag the production of the SMH.
• Many rows and columns are extremely condensed and do not give the full flavour of the signatures explored. For instance since Z bosons are explored in many decay signatures, such as 2e“ or bbì„ and they appear in both the production tag and decay mode. On the other hand, there are analyses such as ttì„H with H→multileptons, where there are contributions from ZZ, WW and ττ which cannot at present be discerned.
• Open stars (a˜†) are taken to mean that no observation has been made in this particular channel and mostly only limits were set.

Now, from the scan you could read the boundaries of the ±1σ interval by looking at where q(μ) crosses the horizontal line at 1. However, since that's a very common task, for models with just one parameter combine can do that for you using the MaxLikelihoodFit module

• Calculate the scale factors for σ·BR of ggH with H→ZZ; VBF with H→WW, VBF with H→γγ; VH with H→bb; ttH with H→bb. Express the results in terms of κ_V, κ_f, κ_H, α, β.

combine c6_comb.root -M MultiDimFit --algo=singles -m 125.7 $RANGES

text2workspace.py comb_hww.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_hww.root --for-fit

HWW="hww0j_8TeV=hwwof_0j_shape_8TeV.txt hww2j_8TeV=hwwof_2j_shape_8TeV.txt"

<--  WE DEFINE HERE THE PATH TO THE GITHUB, USED FOR LINKS LATER 

 Set GITHUB_COMBINE_CMSDAS = https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/tree/CMSDAS-2014-CERN
 
-->

Higgs Combination and Properties: Long Exercise - CMSDAS 2014

Introduction

Overview

The purpose of this long exercise is to understand how some of the key properties of the Higgs boson are measured at CMS: how the couplings are determined from the combination of the results of the searches in different final states (decay modes, production topologies), and how exotic J^CP hypotheses are tested and excluded from the kinematics of the Higgs boson decay products.

For this exercise, we will use the data from the real higgs searches in CMS, and analyze them with the standard CMS statistical toolkit that has been developed for this purpose. In order to keep the amount of computing time down to a manageable level, we will use only the data from the most sensitive channels from each analysis (but including also new channels that were not present at the time of the latest CMS combination), and a version of the statistical tools containing the newest optimized code that is not yet in the production release.

Prerequisites

• Some familiarity with RooFit and basic statistics (e.g. what's a likelihood function) are needed.
• The RooStats short exercise is very strongly recommended

Course Structure Outline

• Introduction
  • The Higgs at CMS: a very brief summary of the higgs physics and of what searches are performed at CMS
  • Datacards: in order to perform a statistical analysis of the results, datacards are used to encode the outcomes of the Higgs searches (i.e. the observed and expected events after the selection, and distribution in whatever variables are used in the analysis to discriminate between signal and background). a description of how these are built and interpreted will be given.
• Higgs Couplings
  • The basics: how to define a simple model to interpret higgs search results fitting some parameters of interest (e.g. the signal strength), how to use basic statistics to define the uncertainties on the extracted parameters in one and two dimensions.
  • First fits to Higgs search results: use what learned before to perform a first fit to the CMS datacards to produce the signal strength plots like these two
    
  • Higgs Coupling models: models to interpret the data for a Higgs with couplings different from the SM ones, and how they are realized in the statistical toolkit
  • Coupling fits from the CMS combination: finally, you will combine all the CMS datacards and fit them to evaluate if the Higgs boson we see is compatible with the SM one or not.
• Higgs J_CP
  • Introduction: the basics of the J_CP studies, and how we can use H→WW→2l2ν to test if the Higgs boson we found is a scalar or a graviton
  • Kinematics: we will use simulated events for a SM Higgs boson and a Graviton to study the kinematics of the events, and understand how the discrimination between the two hypothesis is done
  • Statistics: the basics of the statistical methodology used to separate two hypotheses will be explained
  • Practice: using datacards from the H→WW→2l2ν analysis, we will see if the data favours the SM-like scalar hypothesis or the more exotic graviton one.

Miscellaneous Notes

The color scheme of the Exercise is as follows:

• Shell commands will be embedded in grey box, e.g.

  combine -M MaxLikelihoodFit --justFit comb_hww.txt -m 125.7 --robustFit=1

• Datacards will be embedded in a yellow box, e.g.

  observation 70920

• Python code will be embedded in light blue box, e.g.

  %CODE{"python"}% return self.getHiggsSignalYieldScale(processSource, foundDecay, foundEnergy) %ENDCODE%

• ROOT commands and C++ macro will be embedded in light green box, e.g.

  %CODE{"cpp"}% limit->Draw("2*deltaNLL:r"); %ENDCODE%

• Quizes and important messages will be highlighted in red

Intro I: The Higgs at CMS

Measurements of the Higgs properties in CMS exploit the different production and decay modes expected for the SMH boson.

The SMH boson is produced through gluon-fusion (ggH, red), vector-boson fusion (VBF, blue), associated vector-boson production (VH, purple - where the Higgs is radiated off an energetic V=W,Z boson), and in associated production with a top-quark pair (ttH, green):

As can be seen, the latter three modes have very clear signatures that can be targeted experimentally by, e.g., requiring vector boson candidates in the event or top quark candidates.

One note should be made about the gluon-fusion process: since the SMH boson does not couple to massless particles (like the gluon), the process is loop-induced, with the particles running in the loop being those which couple to the SMH boson and to the gluon. Since the gluon knows about colour charge, the SM particles that can run in the loop are quarks. Since the SMH coupling to quarks is proportional to their mass, the most important contributions to the gluon-fusion cross-section come from the top-loop and bottom-loops, with the contributions form other quarks being negligible. In this way one can immediately see that the gluon-fusion process is sensitive to the top and bottom couplings to the SMH, but could also be affected by BSM particles.

At the LHC, given the very large gluon luminosities, the SMH boson is mainly produced through gluon-fusion and, at much lower rates, the more distinctive processes:

At a mass of m_H ~ 125 GeV, the SMH boson decays into almost every possibly mode except into a pair of top quarks:

Here, the decay into a pair of photons is also loop-mediated as in the case of gluon-fusion production, except that in this case, the photon is sensitive not only to the colour-charge but also to the electric-charge:

While the top and W contributions to the diphoton decay, contributions have also been calculated for other quarks and even the tau lepton. In this way, the diphoton decay branching fraction is sensitive no only to the SMH coupling to those SM particles but as well to BSM particles.

In terms of decay rates, although they seem very small for some decays, in practice and in order to observe the ZZ or WW decays one needs to take into account the branching fractions of the Z and W. After that, the picture is a little different:

Finally, while it is somewhat straightforward to identify the decay signatures and experimentally select pure samples of H→γγ or H→ττ, it is much more complicated to distinguish between production modes. In an nutshell, only leptonic decays of Vs in VH production are extremely pure with the present amount of data. Both VBF and ttì„H selections inevitably include some gluon-fusion events where there was one or more jets radiated from the gluon or quark lines. With more data, stricter cuts can be applied and purer production tags developed.

But for the reasons above we usually refer to "dijet tag" and not "VBF production" when discussing the experimental selections, since such selections will contains some events other than VBF.

In CMS, a large amount of production and decay modes are studied exploring different experimental selections targeting the different production modes. The following table summarises the existing and planned analyses:

There are several notes to this table:

• The H→Zγ decay is a rare, loop-mediated, decay process which is also sensitive to BSM physics as the ggH and H→γγ, with the addition that the Z knows left from right, given the V-A structure of the electroweak interaction.
• In the SM one expects an invisible decay of the SMH from H→ZZ→2ν2ν at the level of 0.1% for m_H=125 GeV. In BSM alternative, the Higgs bosons could decay into massive invisible particles with m_χ<m_H/2 at a substantial rate. From searching those invisible decays one can set stringent limits on this invisible possiblity.
• Combinations denoted with the skull (a˜ ) are not thought to be possible to measure mostly because of the lack of distinctive features with which to tag the production of the SMH.
• Many rows and columns are extremely condensed and do not give the full flavour of the signatures explored. For instance since Z bosons are explored in many decay signatures, such as 2e“ or bbì„ and they appear in both the production tag and decay mode. On the other hand, there are analyses such as ttì„H with H→multileptons, where there are contributions from ZZ, WW and ττ which cannot at present be discerned.
• Open stars (a˜†) are taken to mean that no observation has been made in this particular channel and mostly only limits were set.

But not all ends with the measurements of σ×BR.

In fact there is a much richer and complex structure in terms of the tensor couplings. This is a more difficult arena where the spin (J) and charge-parity (CP) properties of the boson are involved.

The only channel where all information about the boson is available in the detector is the H→ZZ→4e“ where the leptons are electrons or muons (and not neutrinos). This does not mean that other channels cannot contribute to this picture (like H→WW→e“νe“ν or H→γγ), but the fact that all information is readily available renders the analysis much more straightforward in the H→ZZ→4e“ channels.

The following is an illustration from http://arxiv.org/pdf/1208.4018.pdf :

Since the spin-parity tensor structure is much more complicated than the scalar structure of the couplings, much more data will be needed to perform precise measurements. For the moment the strategy is to perform hypotheses tests between the 0⁺ (SMH-like) hypothesis and other, alternative, hypotheses such as pseudoscalar (0^-) bosons (or admixtures), spin-2 bosons, like the graviton, (2⁺), or exotic spin-1 bosons.

As an example, the following figure illustrates the differences in one observable of the ZZ (left), WW (middle), and γγ (right) decays for different J^CP hypotheses (0⁺ red circles, 0^- magenta squares, and two types of spin-2 particle in green and blue):

Intro 2: Datacards

Setting up your computing environment

The Higgs statistical code

The statistical toolkit for the Higgs analyses is provided as a package that is compiled within CMSSW, even if it does not depend on the rest of the CMSSW packages. Tthe code relies critically on the RooFit and RooStats packages, and at the moment it is compatible with the versions packaged in the 6.1.1 and 6.1.2 releases (older releases have bugs, newer releases are backwards-incompatible, not tested and not supported)

The code can be retrieved from the github code hosting side. The code for the development version that we'll be using for this CMSDAS school is visible at https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/tree/CMSDAS-2014-CERN

Log into a SLC5 machine (lxplus5.cern.ch) and do

scramv1 project CMSSW CMSSW_6_1_2      
cd CMSSW_6_1_2/src
cmsenv
git clone https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit.git HiggsAnalysis/CombinedLimit/
cd HiggsAnalysis/CombinedLimit/
git checkout CMSDAS-2014-CERN
scramv1 b -j 5

after the first installation you can instead just do

cd CMSSW_6_1_2/src/HiggsAnalysis/CombinedLimit
cmsenv

The datacards

The datacards contain CMS data, so they cannot be put on github. Normally, they're kept in a protected svn repository at CERN, but for convenience you will be able to get them instead from a tarball on AFS.

cd CMSSW_6_1_2/src/HiggsAnalysis/CombinedLimit
mkdir datacards
cd datacards
tar xvf /afs/cern.ch/user/g/gpetrucc/public/CMSDAS-2014-CERN-cards.tar.gz

A counting experiment: search for invisible higgs decays in the VBF mode

A counting experiment is a search where we just count the number of events that pass a selection, and compare that number with the event yields expected from the signal and background.

For this kind of analysis, all the information needed for the statistical interpretation of the results is encoded in a simple text file An example of search that is done this way is the one for a Higgs boson produced in the VBF mode and decaying to invisible particles (something that does not happen in the standard model but can happen in other scenarios like supersymmetry) You can look at the datacard hinv_vbf_8TeV.txt

Here is a line-by-line commentary of it

# Invisible Higgs analysis for mH=125 GeV

This line is just a comment, describing the datacard. All lines with "#" are comments, and ignored in the processing. Similarly, lines all of minus characters are ignored

imax 1  number of channels
jmax 6  number of backgrounds
kmax 9  number of nuisance parameters (sources of systematical uncertainties)

These lines describe the basics of the datacard: the number of channels, physical processes, and systematical uncertainties Only the first two words (e.g. imax 6) are interpreted, the rest is a comment.

• imax defines the number of final states analyzed (one in this case, but datacards can also contain multiple channels)
• jmax defines the number of independent physical processes whose yields are provided to the code, minus one (i.e. if you have 2 signal processes and 5 background processes, you should put 6)
• kmax defines the number of independent systematical uncertainties (these can also be set to * or -1 which instruct the code to figure out the number from what's in the datacard)

------------
# we have just one channel, in which we observe 0 events
bin vbf
observation 390

These lines describe the number of observed events in each final state. In this case, there's a final state, with label vbf that contains 390 events (the comment above "we observe 0 events" is probably a leftover from an old cut-and-paste...)

------------
bin                          vbf        vbf      vbf     vbf     vbf     vbf     vbf
process                      qqH        zvv      wmu     wel     wtau    qcd     others
process                      0          1        2       3       4       5       6
rate                         207.589    101.879  67.2274 68.2064 54.4137 36.8432 10.3709

These lines describe the number of expected events from each physical process in each final state

The first two rows of each column identify the final state (always "vbf" in this datacard, since it's the only channel considered) and the physics process. Higgs boson signals in the different production modes have special names that combine can recognise: ggH for gluon fusion, qqH for VBF, WH and ZH for associated production with a weak boson (some older datacards that don't separate the two contain VH for the sum of the two, though this is discouraged), and ttH for associated production with a top anti-top pair. Background processes can have any name.

The third line tells the code if a given process is by default to be interpreted as a signal (process number 0 or negative) or as a background (positive process number). For compatibility with other tools that use the same datacard format, in each final state the process numbers should be different for the different process, and increasing from left to right along the columns.

The fourth line contains the expected event yield. In this case, about 200 signal events would be expected (leftmost column) over a background of about 340 events (sum of all the others).

------------
lumi_8TeV             lnN    1.04       -        -       -       -       -       1.04
CMS_vbfhinv_qqh_norm  lnN    1.12981    -        -       -       -       -       -
CMS_vbfhinv_zvv_stat  gmN 12 -          8.4899   -       -       -       -       -
CMS_vbfhinv_zvv_norm  lnN    -          1.25331  -       -       -       -       -
CMS_vbfhinv_wmu_norm  lnN    -          -        1.23748 -       -       -       -
CMS_vbfhinv_wel_norm  lnN    -          -        -       1.29736 -       -       -
CMS_vbfhinv_wtau_norm lnN    -          -        -       -       1.44471 -       -
CMS_vbfhinv_qcd_norm  lnN    -          -        -       -       -       1.84373 -
CMS_vbfhinv_mc_norm   lnN    -          -        -       -       -       -       1.30341

The rest of the datacard describes the systematical uncertainties, and will be explained more in detail below

Systematical uncertainties

The way systematical uncertainties are implemented in the Higgs statistics package is by identifying each independent source of uncertainty, and describing the effect it has on the event yields of the different processes. Each source is identified by a name, and in the statistical model it will be associated with a nuisance parameter.

An indidual source of uncertainty can have an effect on multiple processes, also across different channels, and all these effects will be correlated (e.g. the uncertainty on the production cross section for Higgs from gluon fusion will affect the expected event yield for this process in all datacards). The size of the effect doesn't have to be the same, e.g. a 1% uncertainty on the muon identification efficiency will have a 2% effect on H→WW→2μ2ν but a 4% effect on H→ZZ→4l, and anti-correlations are also possible (e.g. an increase in b-tagging efficiency will simultaneously and coherently increase the signal yield in final states that require tagged b-jets and decrease the signal yield in final states that require no tagged b-jets).

The use of names for each source of uncertainty allow the code to be able to combine multiple datacards recognizing which uncertainties are in common and which are instead separate.

The most common model used for systematical uncertainties is the log-normal, which is identified by the lnN keyword in combine. The distribution is characterized by a parameter κ, and affects the expected event yields in a multiplicative way: a positive deviation of +1σ correspond to a yield scaled by a factor κ compared to the nominal one, while a negative deviation of -1σ correspond to a scaling by a factor 1/κ. For small uncertainties, the log-normal is approximately a gaussian. If δx/x is the relative uncertainty on the yield, κ can be set to 1+δx/x.

For example, the first two uncertaintes of the VBF H→inv datacard are log-normals

bin                          vbf        vbf      vbf     vbf     vbf     vbf     vbf
process                      qqH        zvv      wmu     wel     wtau    qcd     others
process                      0          1        2       3       4       5       6
rate                         207.589    101.879  67.2274 68.2064 54.4137 36.8432 10.3709
------------
lumi_8TeV             lnN    1.04       -        -       -       -       -       1.04
CMS_vbfhinv_qqh_norm  lnN    1.12981    -        -       -       -       -       -

The first one, associated to the normalization of the luminosity in the 8 TeV run ( lumi_8TeV) has a 4% effect on the expected event yield of the signal and of the "others" background (probably this background is estimated directly from simulations, not normalized to data in control regions, and so the expected yield scales with the luminosity). The dash ( -) in the columns of the remaining backgrounds signifies that the uncertainty has no effect on those processes.

The second one, CMS_vbfhinv_qqh_norm is an uncertainty associated to the normalization of only the qqH process, of about 13%.

You can note that the name of the second uncertainty has a prefix CMS_vbfhinv_, while the first one does not: this is used to indicate that the second uncertainty is something pertinent only to the CMS VBF H→inv analysis, while the first uncertainty is common to all CMS or LHC. While the names of the uncertainties don't matter when considering an individual datacard, these naming conventions are very important when combining datacards, to make sure that the effect of an individual source of uncertainty can be tracked across different channels (e.g. other 8 TeV analyses will likely have a lumi_8TeV uncertainty affecting their signal normalization) and that diffent independent sources of uncertainty are not identified in a single one by mistake (which could happen if different datacards used the same name to point to different things).

Another class of systematical uncertainties on the predictions in the signal region are those associated to statistical uncertainties in the control regions used to estimate the backgrounds, or with the limited size of the MC samples used. If the uncertainties are small, Poisson fluctuations become approximately Gaussian, and they be conveniently approximated by log-normal uncertainties with κ = 1 + 1/sqrt(N+1) (where N is the number of events in the control sample or in MC, after the selection). However, if the number of events becomes small, e.g. 10 events or less, but the contribution in the signal region is still significant, then a better statistical treatment can be given by using the Gamma distribution. An example is given below for the Z→νν background, which is estimated using Z→ll events, with a non-negligible statistical uncertainty since only 12 Z→ll events survive the selection.

bin                          vbf        vbf      vbf     vbf     vbf     vbf     vbf
process                      qqH        zvv      wmu     wel     wtau    qcd     others
process                      0          1        2       3       4       5       6
rate                         207.589    101.879  67.2274 68.2064 54.4137 36.8432 10.3709
------------
[...]
CMS_vbfhinv_zvv_stat  gmN 12 -          8.4899   -       -       -       -       -

The gamma distribution, in the form we use it, is characterized by two parameters: N, the number of events in the control region (or in MC), and α the weight by which each event is multiplied to determine the expected yield in the signal region (if the weights are not all the same, one can take α to be the average weight). By definition, the prediciton in the signal region is α·N.

In the datacards, the uncertainty is encoded writing gmN immediately followed by the number of events N, and then putting the factor α in the column corresponding to the process. So, in this case N=12, α=8.5 and the expected yield is about 102 events.

Exercise for the reader

What is the total expected background yield, with its statistical and systematical uncertainties?
For this approximate computation you can treat the log-normals and Gamma distributions as if they were Gaussians (remembering that approximately κ=1+δx/x and δ x/x ~ 1/sqrt(N+1)).

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

A shape analysis using templates: H→ττ

A shape analysis relies not only on the expected event yields for the signals and backgrounds, but also on the distribution of those events in some discriminating variable. This approach is often convenient with respect to a counting experiment, especially when the background cannot be predicted reliably a-priori, since the information from the shape allows a better discrimination between a signal-like and a background-like excess, and provides an in-situ normalization of the background.

The simpler case of shape analysis is the binned one: the observed distribution of the events in data and the expectations from the signal and all backgrounds are provided as histograms, all with the same binning. Mathematically, this is equivalent to just making N counting experiments, one for each bin of the histogram; however, in practice it's much more convenient to provide the the predictions as histograms directly.

For this purpose we will look at the H→ττ search in the τ_μτ_h final state with a VBF di-jet topology. The datacard is htt_mt_5_8TeV.txt, and the full content is:

max    1     number of categories
jmax    9     number of samples minus one
kmax    *     number of nuisance parameters
-------------------------------------------------------------------------------
shapes * * htt_mt.input_8TeV.root $CHANNEL/$PROCESS $CHANNEL/$PROCESS_$SYSTEMATIC
shapes ggH * htt_mt.input_8TeV.root $CHANNEL/$PROCESS$MASS $CHANNEL/$PROCESS$MASS_$SYSTEMATIC
shapes qqH * htt_mt.input_8TeV.root $CHANNEL/$PROCESS$MASS $CHANNEL/$PROCESS$MASS_$SYSTEMATIC
shapes VH * htt_mt.input_8TeV.root $CHANNEL/$PROCESS$MASS $CHANNEL/$PROCESS$MASS_$SYSTEMATIC
-------------------------------------------------------------------------------
bin                                            muTau_vbf
observation                                    199
-------------------------------------------------------------------------------
bin                                            muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf
process                                        -2           -1           0            1            2            3            4            5            6            7
process                                        ggH          qqH          VH           ZTT          QCD          W            ZL           ZJ           TT           VV
rate                                           2.36529      7.04351      0.028046     94.286       39.2563      46.6949      0.828452     3.72816      4.29862      1.48636
-------------------------------------------------------------------------------
lumi_8TeV lnN                                  1.042        1.042        1.042        -            -            -            -            -            -            -
pdf_qqbar lnN                                  -            1.03         1.02         -            -            -            -            -            -            -
pdf_gg lnN                                     1.08         -            -            -            -            -            -            -            -            -
UEPS lnN                                       1.04         1.04         1.04         -            -            -            -            -            -            -
QCDscale_ggH2in lnN                             1.3         -            -            -            -            -            -            -            -            -
QCDscale_qqH lnN                               -            1.04         -            -            -            -            -            -            -            -
QCDscale_VH lnN                                -            -            1.04         -            -            -            -            -            -            -
CMS_eff_m lnN                                  1.02         1.02         1.02         1.020        -            -            1.020        1.020        1.020        1.020
CMS_eff_t_mutau_8TeV lnN                       1.08         1.08         1.08         1.080        -            -            -            -            1.080        1.080
CMS_eff_t_mutau_high_8TeV lnN                  1.015        1.015        1.015        1.015        -            -            -            -            1.015        1.015
CMS_scale_t_mutau_8TeV shape                      1            1            1         1.000        -            -            -            -            -            -
CMS_scale_j_8TeV lnN                            1.1         1.05         1.15         -            -            -            -            -            1.050        1.100
CMS_htt_scale_met_8TeV lnN                     1.05         1.05         1.05         -            -            -            1.050        1.050        1.070        1.070
CMS_htt_zttNorm_8TeV lnN                       -            -            -            1.030        -            -            1.030        1.030        -            -
CMS_htt_ttbarNorm_8TeV lnN                     -            -            -            -            -            -            -            -            1.080        -
CMS_htt_ttbarNorm_mutau_vbf_8TeV lnN           -            -            -            -            -            -            -            -            1.330        -
CMS_htt_DiBosonNorm_8TeV lnN                   -            -            -            -            -            -            -            -            -            1.150
CMS_htt_DiBosonNorm_mutau_vbf_8TeV lnN         -            -            -            -            -            -            -            -            -            1.280
CMS_htt_WNorm_mutau_vbf_8TeV lnN               -            -            -            -            -            1.124        -            -            -            -
CMS_htt_extrap_ztt_mutau_vbf_8TeV lnN          -            -            -            1.100        -            -            -            -            -            -
CMS_htt_QCDSyst_mutau_vbf_8TeV lnN             -            -            -            -            1.190        -            -            -            -            -
CMS_htt_ZLeptonFakeTau_mutau_8TeV lnN          -            -            -            -            -            -            1.300        -            -            -
CMS_htt_ZJetFakeTau_mutau_vbf_8TeV lnN         -            -            -            -            -            -            -            1.400        -            -

Compared to the H→inv datacard, you can see that the format is pretty similar, except for two aspects:

• the lines that start with the shapes keyword at the beginning.
• some systematics have shape as distribution instead of lnN

 *shapes*  _process_  _channel_   _file_    _histogram-name_      _histogram-name-for-systematics_

The first two columns after shapes identify to which process and channel the line refers to; the wildcard * can be used to imply that the rule applies to all processes or to all channels, unless a more specific rule exists for it. The channel names are those used in the bin line in the datacard, and the process names are those in the process line, or data_obs for the observed data.

The third and fourth column are used to point to the shape. The patterns $CHANNEL, $PROCESS, $MASS in the file name or histogram name will be substituted by the name of the channel, the name of the process, and the higgs boson mass hypothesis. The fifth column has a similar use, but is used to find the alternative shapes used to encode systematical uncertainties, and will be described later.

Let's now open the rootfile associated with this datacard in root, and look at the histogram associated to the QCD background. There are no specific shapes rules for this QCD background, so the the first, more generic, rule will apply

shapes * * htt_mt.input_8TeV.root $CHANNEL/$PROCESS $CHANNEL/$PROCESS_$SYSTEMATIC

and for us $CHANNEL is muTau_vbf and $PROCESS is QCD, so we expect the histogram to be in a directory muTau_vbf and be called QCD.

%CODE{"cpp"}% /// root -l root -l htt_mt.input_8TeV.root /// Attaching file htt_mt.input_8TeV.root as _file0... _file0->cd("muTau_vbf"); QCD->Draw("HIST"); %ENDCODE%

Note that we're not drawing the error bars on the histograms, since combine is not using them so it would be misleading to show them.

We can also check the normalization of the histogram, to see that it matches the expected yield in the datacard as it should.

%CODE{"cpp"}% cout << QCD->Integral() << endl; /// ----> 39.2563 %ENDCODE%

Note that underflow and overflow bins are not used by combine.

Similarly we can look at how the 199 events observed in data are distributed: the first rule again will give muTau_vbf/data_obs

%CODE{"cpp"}% cout << data_obs->Integral() << endl; /// ----> 199 data_obs->Draw(); %ENDCODE%

We can also at how the signal looks like for different Higgs boson mass hypotheses; we will use the useful Form("pattern", args) method of ROOT to create the proper string with the mass written inside.

%CODE{"cpp"}% qqH110->Draw("HIST"); for (int m = 110; m <= 140; m += 10) { TH1F h = (TH1F) gDirectory->Get(Form("qqH%d",m)); h->SetLineColor(206+4*(m-110)/10); h->SetLineWidth(2); h->Draw("HIST SAME"); } %ENDCODE%

and in a similar way we can see how do the mean and rms of the signal vary with the mass

%CODE{"cpp"}% for (int m = 110; m <= 140; m += 5) { TH1F h = (TH1F) gDirectory->Get(Form("qqH%d",m)); printf("for mH=%3d, signal mean = %5.1f, rms = %5.1f\n", m, h->GetMean(), h->GetRMS()); } %ENDCODE%

for mH=110, signal mean = 100.5, rms =  15.1
for mH=115, signal mean = 104.6, rms =  16.8
for mH=120, signal mean = 108.8, rms =  17.5
for mH=125, signal mean = 113.2, rms =  17.6
for mH=130, signal mean = 117.3, rms =  18.0
for mH=135, signal mean = 121.3, rms =  19.1
for mH=140, signal mean = 125.5, rms =  21.1

Shape uncertainties using templates

Still using the htt_mt_5_8TeV.txt datacard, we can now consider a systematical uncertainty that affects not just the normalization but also the shape of the expected distribution for a process. The uncertainty on the tau energy scale is a good example of this. The relevant lines of the datacard are

shapes * * htt_mt.input_8TeV.root $CHANNEL/$PROCESS $CHANNEL/$PROCESS_$SYSTEMATIC
[...]
-------------------------------------------------------------------------------
bin                                            muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf    muTau_vbf
process                                        -2           -1           0            1            2            3            4            5            6            7
process                                        ggH          qqH          VH           ZTT          QCD          W            ZL           ZJ           TT           VV
rate                                           2.36529      7.04351      0.028046     94.286       39.2563      46.6949      0.828452     3.72816      4.29862      1.48636
-------------------------------------------------------------------------------
[...]
CMS_scale_t_mutau_8TeV shape                      1            1            1         1.000        -            -            -            -            -            -

This tells us that there is a shape uncertainty called CMS_scale_t_mutau_8TeV which affects the signals and the ZTT background (this background is determined embedding simulated taus into Z→μμ events, so any data/mc discrepancy on the energy response for taus will affect the prediction for it just like it does for the simulated taus in signal MC)

For shape analyses using templates, the effect of a shape-affecting systematical uncertainty is encoded in the datacard by providing two additional templates corresponding to ±Nσ shifts of the corresponding parameter. For each process affected by this uncertainty, the column in the datacard will contain a value 1/N. In this case, as often, ±1σ shapes have been used, and so the columns of the datacard contain a value of 1.

The associated shapes are determined using the last column of the shapes line, which for this datacard reads $CHANNEL/$PROCESS_$SYSTEMATIC, where $SYSTEMATIC is taken to be the name of the uncertainty in the card, and the postixes Up and Down are added to get the up and down templates. Let's see how much does the mean of the expected ZTT distribution change for the ±1σ shifts, and also look at the shapes in ROOT. We can also see that the tau energy scale does not affect only the shape of the background, but also the normalization (probably because if taus are reconstructed as more energetic, they have a higher probability of passing the p_T cuts or other kinematic cuts of the analsysis)

%CODE{"cpp"}% /// root -l htt_mt.input_8TeV.root /// Attaching file htt_mt.input_8TeV.root as _file0... _file0->cd("muTau_vbf"); TH1F nominal = (TH1F) gDirectory->Get("ZTT"); TH1F tau_up = (TH1F) gDirectory->Get("ZTT_CMS_scale_t_mutau_8TeVUp"); TH1F tau_dn = (TH1F) gDirectory->Get("ZTT_CMS_scale_t_mutau_8TeVDown");

printf("Nominal shape: mean = %6.2f (norm: %7.3f events)\n", nominal->GetMean(), nominal->Integral()); // --> Nominal shape: mean = 88.03 (norm: 94.286 events) printf("Up shape: mean = %6.2f (norm: %7.3f events)\n", tau_up->GetMean(), tau_up->Integral()); // --> Up shape: mean = 88.86 (norm: 97.172 events) printf("Down shape: mean = %6.2f (norm: %7.3f events)\n", tau_dn->GetMean(), tau_dn->Integral()); // --> Down shape: mean = 87.20 (norm: 91.995 events)

nominal->SetLineWidth(2); nominal->SetLineColor(kBlack); tau_up->SetLineWidth(2); tau_up->SetLineColor(kBlue); tau_dn->SetLineWidth(2); tau_dn->SetLineColor(kRed); nominal->Draw("HIST"); tau_up->Draw("HIST SAME"); tau_dn->Draw("HIST SAME"); nominal->GetYaxis()->SetRangeUser(0,45); %ENDCODE%

A parametric shape analysis: H→γγ

In some cases, it can be convenient to describe the expected signal and background shapes in terms of analytical functions rather than templates; a typical example are the searches where the signal is apparent as a narrow peak over a smooth continuum background. In this context, uncertainties affecting the shapes of the signal and backgrounds can be implemented naturally as uncertainties on the parameters of those analytical functions. It is also possible to adapt an agnostic approach in which the parameters of the background model are left freely floating in the fit to the data, i.e. only requiring the background to be well described by a smooth function.

Technically, this is implemented by means of the RooFit package, that allows writing generic probability density functions, and saving them into ROOT files. The pdfs can be either taken from RooFit's standard library of functions (e.g. Gaussians, polynomials, ...) or hand-coded in C++, and combined together to form even more complex shapes.

A prototypical case for this kind of analysis is H→γγ analysis. For this excercise, we will use a datacard that contains only the 8 TeV data for four event categories: cat0 and cat1 are categories of untagged events containing good quality diphotons, the former of the two with a higer purity but lower event yield obtained preferentially selecting high p_T diphotons; cat4 and cat5 are categories of di-jet events with different levels of tightness (and so also different level of signal contamination from gluon fusion).

The datacard is the following:

imax 4 number of bins
jmax 5 number of processes minus 1
kmax * number of nuisance parameters
----------------------------------------------------------------------------------------------------------------------------------
shapes WH        cat0      hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_wh_cat0
shapes ZH        cat0      hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_zh_cat0
shapes bkg_mass  cat0      hgg.inputbkgdata_8TeV_MVA.root cms_hgg_workspace:pdf_data_pol_model_8TeV_cat0
shapes data_obs  cat0      hgg.inputbkgdata_8TeV_MVA.root cms_hgg_workspace:roohist_data_mass_cat0
shapes ggH       cat0      hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_ggh_cat0
shapes qqH       cat0      hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_vbf_cat0
shapes ttH       cat0      hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_tth_cat0
[... same as above for cat1, cat4, cat5 ...]
----------------------------------------------------------------------------------------------------------------------------------
bin          cat0         cat1         cat4         cat5       
observation  -1.0         -1.0         -1.0         -1.0       
----------------------------------------------------------------------------------------------------------------------------------
bin                                      cat0         cat0         cat0         cat0         cat0         cat0         cat1         cat1         cat1         cat1         cat1         cat1         cat4         cat4         cat4         cat4         cat4         cat4         cat5         cat5         cat5         cat5         cat5         cat5       
process                                  ZH           qqH          WH           ttH          ggH          bkg_mass     ZH           qqH          WH           ttH          ggH          bkg_mass     ZH           qqH          WH           ttH          ggH          bkg_mass     ZH           qqH          WH           ttH          ggH          bkg_mass   
process                                  -4           -3           -2           -1           0            1            -4           -3           -2           -1           0            1            -4           -3           -2           -1           0            1            -4           -3           -2           -1           0            1          
rate                                     6867.0000    19620.0000   12753.0000   19620.0000   19620.0000   1.0000       7259.4000    19620.0000   12360.6000   19620.0000   19620.0000   1.0000       7063.2000    19620.0000   12556.8000   19620.0000   19620.0000   1.0000       3924.0000    19620.0000   15696.0000   19620.0000   19620.0000   1.0000     
----------------------------------------------------------------------------------------------------------------------------------
CMS_eff_j               lnN              0.999125     0.964688     0.999125     0.998262     0.996483     -            0.999616     0.980982     0.999616     0.99934      0.999012     -            1.02         1.02         1.02         1.02         1.02         -            1.02         1.02         1.02         1.02         1.02         -          
CMS_hgg_JECmigration    lnN              -            -            -            -            -            -            -            -            -            -            -            -            0.853986     0.995971     0.853986     0.846677     0.927283     -            1.025        1.005        1.025        1.025        1.025        -          
CMS_hgg_UEPSmigration   lnN              -            -            -            -            -            -            -            -            -            -            -            -            0.737174     0.991941     0.737174     0.724019     0.86911      -            1.045        1.01         1.045        1.045        1.045        -          
CMS_hgg_eff_MET         lnN              -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -          
CMS_hgg_eff_e           lnN              -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -          
CMS_hgg_eff_m           lnN              -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -            -          
CMS_hgg_eff_trig        lnN              1.01         1.01         1.01         1.01         1.01         -            1.01         1.01         1.01         1.01         1.01         -            1.01         1.01         1.01         1.01         1.01         -            1.01         1.01         1.01         1.01         1.01         -          
CMS_hgg_n_id            lnN              1.034/0.958  1.039/0.949  1.034/0.958  1.053/0.915  1.035/0.958  -            1.042/0.936  1.038/0.948  1.042/0.936  1.053/0.909  1.038/0.954  -            1.016/0.972  1.022/0.963  1.016/0.972  1.017/0.967  1.019/0.964  -            1.024/0.959  1.030/0.948  1.024/0.959  1.014/0.975  1.023/0.962  -          
CMS_hgg_n_pdf_1         lnN              -            0.998/0.996  -            -            1.002/0.998  -            -            0.992/0.999  -            -            1.001/1.000  -            -            0.999/0.998  -            -            1.003/0.999  -            -            0.996/1.000  -            -            1.002/0.999  -          
CMS_hgg_n_pdf_10        lnN              -            1.028/1.004  -            -            0.998/0.998  -            -            1.051/0.997  -            -            1.004/0.999  -            -            1.020/1.000  -            -            1.000/0.998  -            -            1.010/0.999  -            -            0.999/1.000  -          
[... and more rows like this ...]
CMS_hgg_n_sc_gf         lnN              -            -            -            -            0.842/1.123  -            -            -            -            -            0.976/1.031  -            -            -            -            -            0.858/1.095  -            -            -            -            -            0.880/1.083  -          
CMS_hgg_n_sc_vbf        lnN              -            0.993/1.010  -            -            -            -            -            0.994/0.994  -            -            -            -            -            0.997/1.001  -            -            -            -            -            0.999/0.996  -            -            -            -          
CMS_hgg_n_sigmae        lnN              0.950/1.099  0.943/1.114  0.950/1.099  0.956/1.089  0.944/1.112  -            0.954/1.093  0.948/1.104  0.954/1.093  0.971/1.057  0.919/1.165  -            0.996/1.006  0.992/1.016  0.996/1.006  0.995/1.010  0.994/1.012  -            0.991/1.019  0.985/1.031  0.991/1.019  0.995/1.010  0.988/1.026  -          
CMS_id_eff_eb           lnN              1.01999      1.020047     1.01999      1.02002      1.020022     -            1.018535     1.019332     1.018535     1.018642     1.019654     -            1.017762     1.017952     1.017762     1.018824     1.01812      -            1.016723     1.016872     1.016723     1.01884      1.017172     -          
CMS_id_eff_ee           lnN              1.000284     1.000137     1.000284     1.000206     1.0002       -            1.004035     1.001979     1.004035     1.003759     1.001148     -            1.006058     1.005556     1.006058     1.003308     1.005127     -            1.008752     1.008351     1.008752     1.003254     1.007586     -          
JEC                     lnN              0.995187     0.938204     0.995187     0.990442     0.980656     -            0.997891     0.966719     0.997891     0.996368     0.994567     -            1.11         1.035        1.11         1.11         1.11         -            1.11         1.035        1.11         1.11         1.11         -          
QCDscale_VH             lnN              0.982/1.021  -            0.982/1.021  -            -            -            0.982/1.021  -            0.982/1.021  -            -            -            0.982/1.021  -            0.982/1.021  -            -            -            0.982/1.021  -            0.982/1.021  -            -            -          
QCDscale_ggH            lnN              -            -            -            -            0.918/1.076  -            -            -            -            -            0.918/1.076  -            -            -            -            -            0.918/1.076  -            -            -            -            -            0.918/1.076  -          
QCDscale_qqH            lnN              -            0.992/1.003  -            -            -            -            -            0.992/1.003  -            -            -            -            -            0.992/1.003  -            -            -            -            -            0.992/1.003  -            -            -            -          
QCDscale_ttH            lnN              -            -            -            0.906/1.041  -            -            -            -            -            0.906/1.041  -            -            -            -            -            0.906/1.041  -            -            -            -            -            0.906/1.041  -            -          
UEPS                    lnN              0.988624     0.858751     0.988624     0.977408     0.954278     -            0.995014     0.923929     0.995014     0.991416     0.987158     -            1.26         1.08         1.26         1.26         1.26         -            1.26         1.08         1.26         1.26         1.26         -          
lumi_8TeV               lnN              1.044        1.044        1.044        1.044        1.044        -            1.044        1.044        1.044        1.044        1.044        -            1.044        1.044        1.044        1.044        1.044        -            1.044        1.044        1.044        1.044        1.044        -          
pdf_gg                  lnN              -            -            -            0.920/1.080  0.930/1.076  -            -            -            -            0.920/1.080  0.930/1.076  -            -            -            -            0.920/1.080  0.930/1.076  -            -            -            -            0.920/1.080  0.930/1.076  -          
pdf_qqbar               lnN              0.958/1.042  0.972/1.026  0.958/1.042  -            -            -            0.958/1.042  0.972/1.026  0.958/1.042  -            -            -            0.958/1.042  0.972/1.026  0.958/1.042  -            -            -            0.958/1.042  0.972/1.026  0.958/1.042  -            -            -          
vtxEff                  lnN              0.991/1.025  0.989/1.030  0.991/1.025  0.993/1.020  0.989/1.030  -            0.990/1.026  0.990/1.027  0.990/1.026  0.994/1.016  0.984/1.042  -            1.000/0.999  0.997/1.007  1.000/0.999  1.000/1.003  0.998/1.006  -            0.999/1.004  0.998/1.006  0.999/1.004  1.000/1.000  0.997/1.007  -          
CMS_hgg_nuissancedeltamcat4  param  0.0 0.001458
CMS_hgg_nuissancedeltafracright_8TeV  param  1.0 0.002000
CMS_hgg_nuissancedeltamcat1  param  0.0 0.001470
CMS_hgg_nuissancedeltamcat0  param  0.0 0.001530
CMS_hgg_nuissancedeltasmearcat4  param  0.0 0.001122
CMS_hgg_nuissancedeltasmearcat1  param  0.0 0.001167
CMS_hgg_nuissancedeltasmearcat0  param  0.0 0.001230
CMS_hgg_globalscale  param  0.0 0.004717

The first difference compared to the template datacard is in the shapes line; let's take for example these two lines

shapes ggH       cat0      hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_ggh_cat0
shapes data_obs  cat0      hgg.inputbkgdata_8TeV_MVA.root cms_hgg_workspace:roohist_data_mass_cat0

In the datacard using templates, the column after the file name would have been the name of the histogram. Here, instead, we found two names, separated by a colon ( :): the first part identifies the name of the RooWorkspace containing the pdf, and the second part the name of the RooAbsPdf inside it (or, for the observed data, the RooAbsData). Let us remark that parameters of pdfs must be set to be constants, otherwise the combine tool will yield weird results.

Let's inspect this workspace, starting with the data

%CODE{"cpp"}% TFile *fDat = TFile::Open("hgg.inputbkgdata_8TeV_MVA.root"); RooAbsData *data = cms_hgg_workspace->data("roohist_data_mass_cat0"); data->Print("") // --> RooDataHist::roohist_data_mass_cat0[CMS_hgg_mass] = 160 bins (1449 weights)

// so, we have a binned dataset, whose variable is called CMS_hgg_mass: RooRealVar *mass = cms_hgg_workspace->var("CMS_hgg_mass"); mass->Print(""); // RooRealVar::CMS_hgg_mass = 140 L(100 - 180)

// we can make a plot of the dataset with the following RooPlot *plot = mass->frame(); data->plotOn(plot); plot->Draw(); %ENDCODE%

Now let's look also at the signal

%CODE{"cpp"}% TFile *sig = TFile::Open("hgg.inputsig_8TeV_MVA.root"); RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0"); ggH->Print(""); // --> RooAddPdf::hggpdfrel_ggh_cat0[ hist_func_frac_g0_ggh_cat0 * hgg_gaus_g0_ggh_cat0 + hist_func_frac_g1_ggh_cat0 * hgg_gaus_g1_ggh_cat0 + const_func_frac_g2_ggh_cat0 * hgg_gaus_g2_ggh_cat0 ] = 0.604097 // this appears to be a linear combination of multiple gaussian pdfs (hgg_gaus_g0_ggh_cat0, hgg_gaus_g1_ggh_cat0, ...) with coefficents hist_func_frac_g0_ggh_cat0, hist_func_frac_g1_ggh_cat0

// let's get the list of the parameters that describe it RooArgSet *params = ggH->getParameters(*data); params->Print(""); // --> RooArgSet::parameters = (CMS_hgg_globalscale,CMS_hgg_nuissancedeltamcat0,CMS_hgg_nuissancedeltasmearcat0,MH)

// MH is a special parameter, which combine and text2workspace set to the Higgs mass hypothesis. // now since we're just looking into the input workspace, we can set it by hand wsig_8TeV->var("MH")->setVal(125.7);

// Now we can make a plot of the pdf, and show also the contributions from the different gaussians from which it is composed RooPlot *plot = wsig_8TeV->var("CMS_hgg_mass")->frame(); ggH->plotOn(plot); ggH->plotOn(plot, RooFit::Components("hgg_gaus_g0_ggh_cat0"), RooFit::LineColor(kRed)); ggH->plotOn(plot, RooFit::Components("hgg_gaus_g1_ggh_cat0"), RooFit::LineColor(209)); ggH->plotOn(plot, RooFit::Components("hgg_gaus_g2_ggh_cat0"), RooFit::LineColor(222)); plot->Draw(); %ENDCODE%

Parametric signal normalization

There is also another feature of the H→γγ datacard that should catch your eye quicky: the event yields are remarkably strange: the the signal yield is 19620.0000 events for ggH, qqH and ttH, the background yield is 1.0, and observed event yield is -1.0.

Let's start with the simple case: an event yield of -1 just instructs text2workspace and combine to take the yield from the corresponding dataset in the input rootfile, avoiding the need of writing it in the text datacard, but also making the datacard less human-readable. Incidentally, this feature can be used also for datacards that use root histograms like the H→ττ example above.

Now, the signal and background yields: 19k signal events from each production mode with just one background event would be really nice to have, but of course it can't be true. The way the H→γγ works is by relying on an additional feature of text2workspace: in addition to providing generic shapes for the signals and backgrounds, it is also possible to provide generic functions that describe the expected signal and background yields. This additional functions are multiplied by the number in the rate column of the datacard to obtain the final expected event yield. This feature allows e.g. to use the very same datacard to describe all possible different Higgs boson mass hypotheses by just parameterizing properly the expected yield as function of the MH variable.

At present, this feature is not indicated by any special line the datacard: simply, whenever a RooAbsPdf is loaded, text2workspace will search also for a RooAbsReal object with the same name but a _norm postfix, and if present it will use it to scale the event yield.

Armed with this new piece of knowledge, we can now determine what is the expected signal yield for ggH in category 0:

%CODE{"cpp"}% TFile *sig = TFile::Open("hgg.inputsig_8TeV_MVA.root"); wsig_8TeV->var("MH")->setVal(125.7); RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0"); RooAbsReal *ggH_norm = wsig_8TeV->function("hggpdfrel_ggh_cat0_norm"); cout << ggH_norm->getVal()*19620.0000 << endl; // --> 12.3902 %ENDCODE%

This approach can also be used to make a background freely floating, by associating to them as normalization term a floating RooRealVar. However, this can be achieved in a more transparent way by putting instead a normalization uncertainty on that background using a flat pdf: e.g. to leave the background floating between 50% and 200% of its input prediction, a lnU systematic can be used with κ=2 ( lnU has a syntax like lnN, but produces a uniform pdf between 1/κ and κ rather than a log-normal; see SWGuideHiggsAnalysisCombinedLimit)

Shape uncertainties using parameters

The part of the H→γγ datacard related to the systematics starts with many lines of log-normals that should already be familiar, except possibly for the notation with two numbers separated by a slash (e.g. 0.950/1.099). This notation is used for asymmetrical uncertainties: a log-normal with 0.950/1.099 means that at -1σ the yield is scaled down by a factor 0.95, while at +1σ the yield is scaled up by a factor 1.099.

The last part of the datacard contains some lines that use a different syntax, e.g.

CMS_hgg_globalscale  param  0.0 0.004717

These lines encode uncertainties on the parameters of the signal and background pdfs. The example line quoted here informs text2workspace that the parameter CMS_hgg_globalscale is to be assigned a Gaussian uncertainty of ±0.004717 around its mean value of zero (0.0).

The effect can be visualized from RooFit, e.g.

%CODE{"cpp"}% TFile *sig = TFile::Open("hgg.inputsig_8TeV_MVA.root"); wsig_8TeV->var("MH")->setVal(125.7); RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0");

// prepare the canvas RooPlot *plot = wsig_8TeV->var("CMS_hgg_mass")->frame();

// plot nominal pdf ggH->plotOn(plot, RooFit::LineColor(kBlack));

// plot minus 3 sigma pdf wsig_8TeV->var("CMS_hgg_globalscale")->setVal(-3*0.004717); ggH->plotOn(plot, RooFit::LineColor(kBlue));

// plot plus 3 sigma pdf wsig_8TeV->var("CMS_hgg_globalscale")->setVal(+3*0.004717); ggH->plotOn(plot, RooFit::LineColor(kRed)); plot->Draw(); %ENDCODE%

Exercises for the reader

• What is the number of observed events in each category?

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

• Which channel has the highest purity of VBF+VH signal compared to ggH signal?

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

Unbinned, multi-dimensional analysis: H→ZZ

The datacard is identical in form to that of H→γγ, i.e. it provides shapes as RooAbsPdf from a RooWorkspace, and implements shape uncertainties through parameters (for some a fixed range is given to the parameter, e.g. CMS_zz4l_bkgMELA is constrained to be in [-3,3]; this can be useful if for too large deviations of the parameter the shapes start to misbehave)

imax 1
jmax 7
kmax *
------------
shapes * * hzz4l_4muS_8TeV_0.input.root w:$PROCESS
------------
bin a1_0
observation 13
------------
## mass window [105.7,140.7]
bin a1_0 a1_0 a1_0 a1_0 a1_0 a1_0 a1_0 a1_0
process ggH qqH WH ZH ttH bkg2d_qqzz bkg2d_ggzz bkg2d_zjets
process -4 -3 -2 -1 0 1 2 3
rate 1.0000 1.0000 1.0000 1.0000 1.0000 7.2035 0.1360 0.7119
------------
lumi_8TeV lnN 1.044 1.044 1.044 1.044 1.044 1.044 1.044 -
pdf_gg lnN 1.0720 - - - 1.0780 - 1.0710 -
pdf_qqbar lnN - 1.0270 1.0350 1.0349 - 1.0342 - -
pdf_hzz4l_accept lnN 1.02 1.02 1.02 1.02 1.02 - - -
QCDscale_ggH lnN 1.0750 - - - - - - -
QCDscale_qqH lnN - 1.0020 - - - - - -
QCDscale_VH lnN - - 1.0040 1.0155 - - - -
QCDscale_ttH lnN - - - - 1.0655 - - -
QCDscale_ggVV lnN - - - - - - 1.2435 -
QCDscale_VV lnN - - - - - 1.0285 - -
BRhiggs_hzz4l lnN 1.02 1.02 1.02 1.02 1.02 - - -
CMS_eff_m lnN 1.043 1.043 1.043 1.043 1.043 1.043 1.043 -
CMS_hzz4mu_Zjets lnN - - - - - - - 0.6/1.4
CMS_zz4l_bkgMELA param 0  1  [-3,3]
QCDscale_ggH2in lnN - - - - - - - -
QCDscale_qqH lnN - - - - - - - -
CMS_zz4l_ggH_PToverM_sys param 0  1  [-3,3]
CMS_zz4l_qqH_PToverM_sys param 0  1  [-3,3]
CMS_zz4l_ttH_PToverM_sys param 0  1  [-3,3]
CMS_zz4l_WH_PToverM_sys param 0  1  [-3,3]
CMS_zz4l_ZH_PToverM_sys param 0  1  [-3,3]
CMS_zz4l_qqZZ_PToverM_sys param 0  1  [-3,3]
CMS_zz4l_ggZZ_PToverM_sys param 0  1  [-3,3]
CMS_zz4l_ZX_PToverM_sys param 0  1  [-3,3]
CMS_zz4l_mean_m_sig param 0.0 1.0
## CMS_zz4l_mean_m_sig = 0.001
CMS_zz4l_sigma_m_sig param 0.0 0.2
CMS_zz4l_n_sig_1_8 param 0.0 0.01
interf_ggH param 0 1 [-1,1]

The interesting difference however is the observed dataset:

%CODE{"cpp"}% gSystem->Load("libHiggsAnalysisCombinedLimit.so") gFile = TFile::Open("hzz4l_4muS_8TeV_0.input.root"); RooAbsData *data = w->data("data_obs"); data->Print("") // --> RooDataSet::data_obs[CMS_zz4l_mass,melaLD,CMS_zz4l_PToverM] = 13 entries %ENDCODE%

In the Hγγ the dataset was a RooDataHist with one entry per bin, each entry having a weight equal to the number of events in that bin. Here instead it's a RooDataSet, with one entry per event.

Let's explore its contents, both building some binned 1D projections and by making a 2D scatter plot.

%CODE{"cpp"}% w->var("CMS_zz4l_mass")->Print(""); // --> RooRealVar::CMS_zz4l_mass = 114.157 L(104 - 142) B(200) RooPlot *mass = w->var("CMS_zz4l_mass")->frame(RooFit::Bins(19)); data->plotOn(mass); mass->Draw();

RooPlot *mela = w->var("melaLD")->frame(RooFit::Bins(10)); data->plotOn(mass); mass->Draw();

TGraph points(data->numEntries()); for (int i = 0; i < data->numEntries(); ++i) { const RooArgSet *ev = data->get(i); points->SetPoint(i, ev->getRealValue("CMS_zz4l_mass"), ev->getRealValue("melaLD")); } points->SetMarkerStyle(20); points->GetXaxis()->SetTitle("m_{4l} (GeV)") points->GetYaxis()->SetTitle("K_{D}") .x tdrstyle.cc points->Draw("AP"); %ENDCODE%

1D (mass)	1D (MELA KD)	2D

Combining datacards

In order to analyze the outcomes of multiple searches together, it is necessary to combine their datacards into a single one.

This is done by the combineCards.py script, whose syntax is simple:

combineCards.py Name1=card1.txt Name2=card2.txt .... > card.txt

If the input datacards had just one bin each, then the output channels will be called Name1, Name2, and so on. Otherwise, a prefix Name1_ ... Name2_ will be added to the bin labels in each datacard. The supplied bin names Name1, Name2, etc. must themselves conform to valid C++/python identifier syntax.

The package containing the datacards also has a simple script mkcomb.sh which will construct different combinations of input datacards.

For example, for the H→WW search, we will combine the 0-jets opposite-flavour channel and the di-jet opposite-flavour channels, which have respectively the best sensitivity to gluon fusion and vbf higgs production. The combineCards.py $HWW > comb_hww.txt command will be expanded in

combineCards.py hww0j_8TeV=hwwof_0j_shape_8TeV.txt hww2j_8TeV=hwwof_2j_shape_8TeV.txt > comb_hww.txt

The combined datacard looks like the following

Combination of hww0j_8TeV=hwwof_0j_shape_8TeV.txt  hww2j_8TeV=hwwof_2j_shape_8TeV.txt
imax 2 number of bins
jmax 23 number of processes minus 1
kmax 89 number of nuisance parameters
----------------------------------------------------------------------------------------------------------------------------------
shapes *           hww0j_8TeV  hwwof_0j.input_8TeV.root histo_$PROCESS histo_$PROCESS_$SYSTEMATIC
shapes data_obs    hww0j_8TeV  hwwof_0j.input_8TeV.root histo_Data
shapes *           hww2j_8TeV  hww-19.47fb.mH125.of_2j_shape_mll.root histo_$PROCESS histo_$PROCESS_$SYSTEMATIC
shapes data_obs    hww2j_8TeV  hww-19.47fb.mH125.of_2j_shape_mll.root histo_Data
----------------------------------------------------------------------------------------------------------------------------------
bin          hww0j_8TeV  hww2j_8TeV
observation  5729.0      24.0
----------------------------------------------------------------------------------------------------------------------------------
bin                                                          hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww0j_8TeV  hww2j_8TeV  hww2j_8TeV  hww2j_8TeV  hww2j_8TeV  hww2j_8TeV  hww2j_8TeV  hww2j_8TeV  h
process                                                      ZH          qqH         WH          ggH         qqWW        WjetsE      Wgamma      Top         WH_SM       WjetsM      ggWW        ZH_SM       VV          Wg3l        qqH_SM      Ztt         Zjets       ggH_SM      qqH         ggH         WW          WWewk       Vg          WJet        Top         g
process                                                      -3          -2          -1          0           1           2           3           4           5           6           7           8           9           10          11          12          13          14          -2          0           15          16          17          18          4           7
rate                                                         1.7481      3.0668      5.8955      237.1930    3969.6300   282.7930    115.5780    498.7010    0.0000      331.7970    210.6270    0.0000      132.6090    167.8020    0.0000      45.9960     0.0000      0.0000      5.0290      1.4400      4.0230      1.9080      0.0389      2.4540      14.0300     0
----------------------------------------------------------------------------------------------------------------------------------
[...]
CMS_8TeV_hww_Top_2j_of_2j_extr   lnN                         -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           1.188       -
[...]
CMS_hwwof_0j_MVAWHStatBounding_8TeV   shape                  -           -           1.0         -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -           -
[....]
lumi_8TeV               lnN                                  1.044       1.044       1.044       1.044       1.0         -           1.044       -           1.044       -           1.0         1.044       1.044       1.044       1.044       1.044       -           1.044       1.044       1.044       1.044       1.044       1.044       -           -           1

You can see that the bin and observation rows now contain two columns (one for each channel), and the common lumi_8TeV uncertainty row affects the columns of of processes from the different cards, while e.g. CMS_8TeV_hww_Top_2j_of_2j_extr or CMS_hwwof_0j_MVAWHStatBounding_8TeV have no effect on the columns of the other channel.

In the last command you see that an option -S is passed to the combineCards.py program. This option is needed when combining together counting experiments (such as VBF H→inv) and shape analyses, and tells the code to reintepret the counting experiment as if it were a shape analysis.

Higgs Couplings

The basics

From Datacards to Likelihoods

The datacards from each search contain the information about the observed events and the expected yields for the various physical processes, backgrounds and signals. For a SM Higgs of fixed mass, there is no free parameter in the theory, and so nothing to fit.

To search for deviations from the SM Higgs predictions, we define a model that can predict different expected yields for the signal, and fit the parameters of this model to the data.

We will start with the simplest possible model, with only one parameter μ: the predictions for the signal in a given final state in this model will be μ times the yield expected for a SM Higgs (this parameter can be interpreted as a scale factor for all the production cross sections of the signal, so we sometimes write μ = σ/σ_SM)

In the case of a simple counting experiment with no systematical uncertainties, we can write the likelihood function of the observed data as function of the parameter of interest in terms of the observed events ( N), expected signal events for a SM Higgs ( S), and expected background events ( B)
L(\,\mathrm{data}\,|\, \mu\cdot s + b\,) = \mathrm{Poisson}(\,N\,|\,\mu\cdot S + B\,)

If there are systematical uncertainties, they are modelled as additional nuisance parameters θ that affect the expected signal and background, and the likelihood contains some additional constraint term for θ that models our a-priori knowledge of the parameter (e.g. in the case of the normalization of the integrated luminosity for the VBF H→inv datacard given before, the signal and background yields will be proportional to the luminosity and the constranint will encode the fact we knew the luminosity with a relative 4% uncertainty). The likelihood can be written now as
L(\,\mathrm{data}\,|\, \mu\cdot s(\theta) + b(\theta)\,) = \mathrm{Poisson}(\,N\,|\,\mu\cdot S(\theta) + B(\theta)\,)\ \cdot\ L(\theta)

When we include multiple channels, or when we do a binned shape analysis, the likelihood becomes the product over the channels or bins
L(\,\mathrm{data}\,|\, \mu\cdot s(\theta) + b(\theta)\,) = \prod_{i} \huge[ \mathrm{Poisson}(\,N_i\,|\,\mu\cdot S_i(\theta) + B_i(\theta)\,)\huge]\ \cdot\ L(\theta)

For unbinned analyses, there is an additional term that includes a product on all the events of the channel involving the expected probability density functions for the distribution of the signal and background events (f_S and f_B) weighted with the expected events for each
L(\,\mathrm{data}\,|\, \mu\cdot s + b\,) = \mathrm{Poisson}(\,N\,|\,\mu\cdot S + B\,)\ \cdot\ \prod_e \left[ \frac{\mu\cdot S}{\mu\cdot S + B}\,f_S(e) + \frac{B}{\mu\cdot S + B}\,f_B(e) \right]

Likelihoods in combine and RooStats

The first step of the statistical analysis is to convert the text datacard into a likelihood function, in the format used by RooStats. For some simple cases, combine does this automatically when given a text file, but it is possible to do it also manually beforehand - and it's necessary to do it manually whenever we want to use a different model than the one with just μ.

The conversion from datacard to likelihood is done with the text2workspace.py program, also contained in the HiggsAnalysis/CombinedLimit package. The syntax is

text2workspace.py datacard.txt -m mass  [ -P physicsModel [ --PO physicsModelOptions .... ] ]  [other options]   -o output.root

The physicsModel is the reference to the python class that describes what are the parameters of the model (e.g. μ) and how those parameters affect the yields of the signal and background (e.g. that the signal is multiplied by μ while the background not)

The most basic Higgs model, defining a single parameter of interest μ is below. As python and roofit don't like greek letters, μ is called r in the code. The code defines two functions, one that declares the parameters, and one that defines how the expected Higgs boson signal yields for each production mode, decay mode and energy depend on those parameters.

• doParametersOfInterest, which defines the parameters
• getHiggsSignalYieldScale), which defines how the expected Higgs boson signal yields for each production mode, decay mode and energy depend on those parameters. The implementation of the latter is trivial: all signal yields scale as r, so the function always returns r.

def doParametersOfInterest(self): """Create POI and other parameters, and define the POI set.""" # --- Signal Strength as only POI --- self.modelBuilder.doVar("r[1,0,20]"); self.modelBuilder.doSet("POI","r") # --- Higgs Mass as other parameter ---- if self.options.mass = 0: if self.modelBuilder.out.var("MH"): self.modelBuilder.out.var("MH").removeRange() self.modelBuilder.out.var("MH").setVal(self.options.mass) else: self.modelBuilder.doVar("MH[%g]" % self.options.mass);

def getHiggsSignalYieldScale(self,production,decay, energy): return "r" %ENDCODE%

This class StrictSMLikeHiggsModel is defined in the HiggsAnalysis.CombinedLimit.PhysicsModel file, with name strictSMLikeHiggs (see the lines at the bottom of the file), so the argument to pass to the -P option will be HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs

Let's use this model to convert the H→WW + 2 jets datacard into a likelihood

text2workspace.py hwwof_2j_shape_8TeV.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs  -o hww2j.root

the beginning of the datacard was

----------------------------------------------------------------------------------------------------
bin         of_2j
observation 24
shapes  *           * hww-19.47fb.mH125.of_2j_shape_mll.root     histo_$PROCESS histo_$PROCESS_$SYSTEMATIC
shapes  data_obs    * hww-19.47fb.mH125.of_2j_shape_mll.root     histo_Data
----------------------------------------------------------------------------------------------------
bin                                                       of_2j     of_2j     of_2j     of_2j     of_2j     of_2j     of_2j     of_2j     of_2j     of_2j     of_2j
process                                                   ggH       qqH       ggWW      VgS       Vg        WJet      Top       WW        WWewk     VV        DYTT
process                                                   0         -1        1         2         3         4         5         6         7         8         9
                                         rate          1.44     5.029    0.4295    0.0195    0.0389     2.454     14.03     4.023     1.908    0.4833     1.879
----------------------------------------------------------------------------------------------------

Now we can open the hww2j.root file in ROOT and see its contents. Start a root prompt with the command root.exe -b -l -n and

%CODE{"cpp"}% // Load the libraries needed to interpret the likelihood gSystem->Load("libHiggsAnalysisCombinedLimit.so"); // Open the file gFile = TFile::Open("hww2j.root"); // Retrieve the workspace from the file RooWorkspace *w = (RooWorkspace *) gFile->Get("w"); // Retrieve the ModelConfig from the workspace RooStats::ModelConfig *mc = (RooStats::ModelConfig *) w->genobj("ModelConfig");

// Now we get the list of parameters of the workspace mc->GetParametersOfInterest()->Print("V"); // the output will be // 1) RooRealVar:: r = 1

// Now we get the expected signal yield for ggH and DYTT cout << w->function("n_exp_final_binof_2j_proc_ggH")->getVal() << endl; cout << w->function("n_exp_final_binof_2j_proc_DYTT")->getVal() << endl; // the output will be 1.44, and 1.879, as expected

// Now we change the value of r to be 2.0... we'd expect the ggH yield to grow to 2.88 and the DYTT to remain the same w->var("r")->setVal(2.0); cout << w->function("n_exp_final_binof_2j_proc_ggH")->getVal() << endl; // ----> 2.88 cout << w->function("n_exp_final_binof_2j_proc_DYTT")->getVal() << endl; // ----> 1.879

// We can also get the number of observed events, that should match the "observation" line in the datacard cout << w->data("data_obs")->sumEntries() << endl; // ----> 24 %ENDCODE%

One-dimensional fits and profile likelihood

Now that we have written the likelihood as function of a parameter μ we can proceed further to see how well different values of μ describe the observed data.

The first step is to find the value of μ that describes the data best, i.e. that maximises the likelihood. This value of μ will be denoted as \hat{\mu}. To ask combine to compute this value, we can use the MultiDimFit method, without specifying any option. We will now use the datacard with all H→WW channels, and tell combine to call the output "HWW_FIT"

text2workspace.py comb_hww.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs  -o comb_hww.root
combine -M MultiDimFit comb_hww.root -m 125.7 -n HWW_FIT

The output of the program will be

>>> including systematics
>>> method used to compute upper limit is MultiDimFit
>>> random number generator seed is 123456
Computing limit starting from observation

 --- MultiDimFit ---
best fit parameter values:
   r :    +0.599
Done in 0.01 min (cpu), 0.02 min (real)

The output file from combine will contain a tree with a single entry, corresponding to this value.

%CODE{"cpp"}% // root.exe -b -l higgsCombineHWW_FIT.MultiDimFit.mH125.7.root // Attaching file higgsCombineHWW_FIT.MultiDimFit.mH125.7.root as _file0... root [1] limit->Scan("r"); //************************ //* Row * r * //************************ //* 0 * 0.5986483 * //************************ %ENDCODE%

Now, to compare different values of μ, we will use the negative log of the profile likelihood function λ
q(\mu) = - 2 \ln \lambda(\mu) = - 2 \ln \frac{L(\,\mathrm{data}\,|\, \mu\cdot s(\hat{\theta}_\mu) + b(\hat{\theta}_\mu)\,)}{L(\,\mathrm{data}\,|\, \hat{\mu}\cdot s(\hat{\theta}) + b(\hat{\theta})\,)}

• The denominator of the ratio is the likelihood, maximised as function of μ and also of the nuisance parameters θ
• The numerator of rhe ratio is the likelihood function for a fixed value of μ; the nuisance parameters θ are set to the values that maximise the likelihood for that fixed μ

We can now use combine to compute q(μ) for different values of μ. Using the same input file as before, let's scan 100 values of mu along a line (i.e. a 1-dimensional grid), in the range between 0 and 3. We will tell combine to call the output file HWW_SCAN_MU

combine -M MultiDimFit comb_hww.root -m 125.7 -n HWW_SCAN_MU --algo=grid --points=100 --setPhysicsModelParameterRanges r=0,3

This command will be slower than the previous one, since the code now has to maximise the likelihood 101 times: once first to get the best fit μ then 100 times one for each point along the grid. Some of those fits will not converge at the first attempt, so the putput could contain some lines like Minimization did NOT converge, Edm is above max or Covar was made pos def but this is normally ok.

We can inspect the output from ROOT: the output tree in the file higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root will contain one entry for each of the 101 fits,

%CODE{"cpp"}% // root -l higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root

limit->GetEntries() // (const Long64_t)101

limit->Scan("r:2*deltaNLL") //************************************ //* Row * r * 2*deltaNL * //************************************ //* 0 * 0.5986473 * 0 * //* 1 * 0.0149999 * 5.9730005 * //* 2 * 0.0450000 * 5.3444776 * //* 3 * 0.0750000 * 4.7531828 * //* 4 * 0.1049999 * 4.1995768 * //* 5 * 0.1350000 * 3.6845009 *

limit->Draw("2*deltaNLL:r"); %ENDCODE%

The quality of the default plot looks awful, so we will retrieve that information as a TGraph object and make it a bit more pretty; also, we will zoom into the region closer to the minimum which is more interesting.

%CODE{"cpp"}% // Draw limit->Draw("2*deltaNLL:r");

// Make a copy of the graph TGraph *graph = gROOT->FindObject("Graph")->Clone();

// sort the points for increasing values of r graph->Sort();

// set the style graph->SetLineWidth(2); graph->SetLineColor(214); graph->SetMarkerStyle(7); graph->SetTitle("Likelihood scan as function of #mu;#mu;-2 #Delta ln L"); graph->GetXaxis()->SetRangeUser(0,2); graph->GetYaxis()->SetRangeUser(0,8); c1->SetGridy();

// draw as a smooth line plus the points graph->Draw("APC"); %ENDCODE%

Now, from the scan you could read the boundaries of the ±1σ interval by looking at where q(μ) crosses the horizontal line at 1. However, since that's a very common task, for models with just one parameter combine can do that for you using the MaxLikelihoodFit module

combine -M MaxLikelihoodFit comb_hww.root -m 125.7 -n HWW_FIT_MU --justFit --robustFit=1

The two options on the command line instruct combine to do just a fit for μ ( --justFit), using a slow but robust algorithm ( --robustFit=1; without that, for a datacard with many nuisance parameters the calculation of the uncertainties will often fail with error messages like Error in : Minuit2Minimizer::MINOS failed due to invalid function minimum or Error in : Minos error calculation failed for all parameters). In the output, you can see various values of r being tried, until it finds the right ones

>>> including systematics
>>> method used to compute upper limit is MaxLikelihoodFit
>>> random number generator seed is 123456
Computing limit starting from observation
Searching for crossing at nll = -18788 in the interval 0.598648, 20
      r      lvl-here  lvl-there   stepping
2.538783    -11.83519  +0.49447    1.940135
1.180689    -1.63413  -11.83519    0.582041
0.773260    +0.23757  -1.63413    0.174612
0.947873    -0.38172  +0.23757    0.174612
0.825644    +0.07707  -0.38172    0.052384
0.878028    -0.13010  +0.07707    0.052384
0.841359    +0.02136  -0.13010    0.015715
0.857074    -0.04154  +0.02136    0.015715
0.846074    +0.00325  -0.04154    0.004715
0.850788    -0.01552  +0.00325    0.004715
0.847488    -0.00232  -0.01552    0.001414
0.846498    +0.00158  -0.00232    0.000424
0.846922    -0.00008  +0.00158    0.000424
Searching for crossing at nll = -18788 in the interval 0.598648, 0
      r      lvl-here  lvl-there   stepping
0.538783    +0.45835  +0.49447    -0.059865
0.478919    +0.35676  +0.45835    -0.059865
0.419054    +0.19666  +0.35676    -0.059865
0.359189    -0.04147  +0.19666    -0.059865
0.401094    +0.13250  -0.04147    -0.017959
0.383135    +0.06353  +0.13250    -0.017959
0.365175    -0.01304  +0.06353    -0.017959
0.377747    +0.04198  -0.01304    -0.005388
0.372359    +0.01921  +0.04198    -0.005388
0.366971    -0.00477  +0.01921    -0.005388
0.370743    +0.01215  -0.00477    -0.001616
0.369127    +0.00497  +0.01215    -0.001616
0.367510    -0.00232  +0.00497    -0.001616
0.368642    +0.00280  -0.00232    -0.000485
0.368157    +0.00061  +0.00280    -0.000485
0.367672    -0.00158  +0.00061    -0.000485
0.368011    -0.00005  -0.00158    -0.000145

 --- MaxLikelihoodFit ---
Best fit r: 0.598648  -0.230684/+0.248194  (68% CL)
nll S+B -> -8.65939  nll B -> -1
Done in 0.11 min (cpu), 0.12 min (real)

So, from this combination we have measured μ =0.60 +0.25/–0.23. Since to keep things simpler we didn't use all the datacards of the H→WW analysis but only a few ones, the result is different compared to the one to appear in the legacy H→WW paper and the error bars are a bit larger, but we are in the right ballpark.

Exercise for the reader

Do a fit of the signal strength separately for the 0-jet and 2-jet datacards of H→WW, and compute the weighted average using the textbook formulas for averaging numbers accounting for their uncertainties (for this approximate computation, you can deal with the asymmetric uncertainties by just taking the one that goes in the direction of the average result). How different it is from what obtained from the combined datacard?

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

text2workspace.py -m 125.7  -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs hwwof_0j_shape_8TeV.txt -o hww0j.root
text2workspace.py -m 125.7  -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs hwwof_2j_shape_8TeV.txt -o hww2j.root

<--/twistyPlugin-->

Expectations: the Asimov dataset

Determining what is the expected outcome of a fit if the data were to match the SM Higgs expectations can be useful for various purposes, e.g. which to quantify the sensitivity of a given analysis in a way that is independent of the presence of fluctuations in the observed data (or when the data is still blinded, or not yet collected). The behaviour ofa fit with a complex model on the Asimov is also easier to understand with respect to the fit to the observed data, especially if the data departs significantly from the expectations.

Normally, the procedure used to determine these expectations is to generate a large sample of possible experimental outcomes, all under the SM Higgs hypothesis, and use those to determine the median outcome of the analysis. This however would require substantial computational resources. Fortunately, a much cheaper but good approximation exists, through what is called the "Asimov dataset"

In a parallel reality imagined by the science fiction writer I. Asimov, politics was run in a peculiar way: instead of mobilizing millions of people to cast their vote to deliberate on something, an algorithm was used to select an individual "average" person, and then this person was asked to take the decision on that matter.

In more practical terms, an Asimov dataset can be defined as the one that matches perfectly the median expectations from the SM Higgs (i.e. the event yields in each bin are exactly as expected). The idea is simple and intuitive, and apart from a few subtleties like the need of having a non-integer number of observed events in each bin and of dealing efficiently with multi-dimensional models, which we won't discuss in this exercise.

When running on a datacard, combine can be instructed to use an Asimov dataset instead of the observed one simply by providing the options -t -1 --expectSignal=1 (the last of which instructs the code to create an Asimov for signal+background and not for background-only).

We can for example produce the expected likelihood scan vs μ for H→WW with

combine -M MultiDimFit comb_hww.root -m 125.7 -n HWW_SCAN_MU_ASIMOV --algo=grid --points=100 --setPhysicsModelParameterRanges r=0,3 -t -1 --expectSignal=1

A comparison of the two likelihood scans, observed and Asimov, is below: the observed data is shown as a thick line, the Asimov as a light dotted line. The shapes are similar, but the Asimov is centered at μ=1 (by construction, since that's the SM expectation), while the observed data is not.

Show code to make it Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

Two-dimensional fits: μ_ggH,ttH vs μ_VBF,VH

Now we will look at a more elaborate model, which has two parameter:

• μ_ggH,ttH scaling the Higgs boson production cross sections in modes relying on the fermion couplings (gluon fusion, which is mediated mostly by top quarks; ttH, in which the Higgs is radiated by a top quark)
• μ_VBF,VH scaling the cross section in modes related to the coupling to electroweak gauge bosons (VBF, WH, and ZH, all of which involve the vertex with a Higgs boson and two gauge bosons).

RvRfXSHiggs

%CODE{"python"}% class RvRfXSHiggs(SMLikeHiggsModel): "Float ggH and ttH together and VH and qqH together" def __init__(self): SMLikeHiggsModel.__init__(self) # not using 'super(x,self).__init__' since I don't understand it self.floatMass = False

def setPhysicsOptions(self,physOptions): ## [ ... code for floating the Higgs mass, irrelevant for us now .... ]

def doParametersOfInterest(self): """Create POI out of signal strength and MH""" # --- Signal Strength as only POI --- self.modelBuilder.doVar("RV[1,-5,15]") self.modelBuilder.doVar("RF[1,-4,8]") if self.floatMass: # [ ... not interesting ... ] else: # [ ... not interesting ... ] self.modelBuilder.doSet("POI",'RV,RF')

def getHiggsSignalYieldScale(self,production,decay, energy): if production in ['ggH', 'ttH']: return 'RF' if production in ['qqH', 'WH', 'ZH', 'VH']: return 'RV' raise RuntimeError, "Unknown production mode '%s'" % production %ENDCODE%

Analogously to the one-dimensional case, we can define a likelihood function that depends on μ_ggH,ttH ( RF in the code) and μ_VBF,VH ( RV in the code), search for the maximum, and define a q(μ_ggH,ttH,μ_VBF,VH) from the profile likelihood
q(\mu_\mathrm{ggH,ttH}, \mu_\mathrm{VBF,VH}) = - 2 \ln \lambda(\mu) = - 2 \ln \frac{L(\,\mathrm{data}\,|\, \mu_\mathrm{ggH,ttH}\cdot s_\mathrm{ggH,ttH} + \mu_\mathrm{VBF,VH}\cdot s_\mathrm{VBF,VH} + b\,)}{L(\,\mathrm{data}\,|\, \hat{\mu}_\mathrm{ggH,ttH}\cdot s_\mathrm{ggH,ttH} + \hat{\mu}_\mathrm{VBF,VH}\cdot s_\mathrm{VBF,VH} + b\,)}

Let's apply this analysis to the H→γγ search, to find first the best fit point and then do a two-dimensional scan of q.

text2workspace.py comb_hgg.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.PhysicsModel:rVrFXSHiggs  -o rVrF_hgg.root
combine -M MultiDimFit rVrF_hgg.root -m 125.7 -n HGG_RVRF_FIT

>>> including systematics
>>> method used to compute upper limit is MultiDimFit
>>> random number generator seed is 123456
ModelConfig 'ModelConfig' defines more than one parameter of interest. This is not supported in some statistical methods.
Computing limit starting from observation

 --- MultiDimFit ---
best fit parameter values:
   RV :    +0.259
   RF :    +1.229
Done in 0.01 min (cpu), 0.01 min (real)

Now, for the 2D scan we'll do 400 points (20x20), using the ranges [-1,3] and [-2,5] for the two parameters. This requires 400+1 fits, so it will take a few minutes.

combine -M MultiDimFit rVrF_hgg.root -m 125.7 -n HGG_RVRF_SCAN_2D --algo=grid --points=400 --setPhysicsModelParameterRanges RF=-1,3:RV=-2,5

>>> including systematics
>>> method used to compute upper limit is MultiDimFit
>>> random number generator seed is 123456
ModelConfig 'ModelConfig' defines more than one parameter of interest. This is not supported in some statistical methods.
Set Range of Parameter RF To : (-1,3)
Set Range of Parameter RV To : (-2,5)
Computing limit starting from observation
Point 0/400, (i,j) = (0,0), RV = -1.825000, RF = -0.900000
Point 2/400, (i,j) = (0,2), RV = -1.825000, RF = -0.500000
Point 4/400, (i,j) = (0,4), RV = -1.825000, RF = -0.100000
Point 6/400, (i,j) = (0,6), RV = -1.825000, RF = 0.300000
Point 8/400, (i,j) = (0,8), RV = -1.825000, RF = 0.700000
Point 10/400, (i,j) = (0,10), RV = -1.825000, RF = 1.100000
Point 12/400, (i,j) = (0,12), RV = -1.825000, RF = 1.500000
[...]
Point 388/400, (i,j) = (19,8), RV = 4.825000, RF = 0.700000
Point 390/400, (i,j) = (19,10), RV = 4.825000, RF = 1.100000
Point 392/400, (i,j) = (19,12), RV = 4.825000, RF = 1.500000
Point 394/400, (i,j) = (19,14), RV = 4.825000, RF = 1.900000
Point 396/400, (i,j) = (19,16), RV = 4.825000, RF = 2.300000
Point 398/400, (i,j) = (19,18), RV = 4.825000, RF = 2.700000
Done in 1.52 min (cpu), 1.53 min (real)

A simple plot from the output of the bidimensional scan can be done with

%CODE{"cpp"}% // plot the likelihood in the 20x20 grid, selecting the interesting region of small deltaNLL limit->Draw("2*deltaNLL:RV:RF>>h2(20,-1,3,20,-2,5)","deltaNLL > 0 && deltaNLL < 11","PROF COLZ");

// make the plot less horrible gStyle->SetOptStat(0); h2->SetContour(200); h2->SetTitle("Profile Likelihood in 2D;#mu_{ggH,ttH};#mu_{VBF,VH}");

/// Now we add the best fit point limit->Draw("RV:RF","deltaNLL == 0","P SAME"); TGraph *bestFit = gROOT->FindObject("Graph"); bestFit->SetMarkerStyle(34); bestFit->SetMarkerSize(3); c1->Update(); %ENDCODE%

Contours for different confidence levels can be constructed starting from the 2d likelihood scan, using the thresholds for a chisquare with two degrees of freedom (pdf review). In particular, the 1σ region is given by q ≤ 2.30, the 95% one by q ≤ 5.99 and the 3σ one by q ≤ 11.83; the C++ function ROOT::Math::chisquared_quantile(cl,ndf) can be used to compute the thresholds for each CL.

ROOT provides a way to extract a smooth contour out of a two-dimensional histogram, and this functionality is packaged in a script test/plotting/contours2D.cxx. The top level function, that you will call directly, is

%CODE{"cpp"}% void contour2D(TString xvar, int xbins, float xmin, float xmax, TString yvar, int ybins, float ymin, float ymax, float smx=1.0, float smy=1.0, TFile *fOut=0, TString name="contour2D") ; %ENDCODE%

• gFile should point to the TFile containing the limit TTree
• xvar should be the variable to use on the X axis, with xbins bins in the [xmin, xmax] range
• yvar should be the variable to use on the Y axis, with ybins bins in the [ymin, ymax] range
• (smx, smy) are the coordinates at which to put a diamond representing the SM expectation
• if fOut is not null, then the output objects will be saved to fOut (see example later)

We can now use this macro interactively to make a nicer plot of the previous 2D scan. For convenience, copy the contours2D.cxx script in your working directory, and then

%CODE{"cpp"}% // root -l higgsCombineHGG_RVRF_SCAN_2D.MultiDimFit.mH125.7.root contours2D.cxx // Attaching file higgsCombineHGG_RVRF_SCAN_2D.MultiDimFit.mH125.7.root as _file0... // Processing contours2D.cxx...

// before we had done // limit->Draw("2*deltaNLL:RV:RF>>h2(20,-1,3,20,-2,5)","deltaNLL > 0 && deltaNLL < 11","PROF COLZ"); // now we do contour2D("RF",20,-1.,3., "RV",20,-2.,5., 1.,1.); %ENDCODE%

Fitting only a subset of the parameters of a model

Sometimes, one is interested in determining confidence regions for one or two parameters in a model that has more parameters. There are two conceptually different ways in which this can be done: the other parameters can be either kept fixed at some nominal value (e.g. the SM prediction), or they can be profiled (i.e. left freely floating in the fit)

Still using the H→γγ analysis with the μ_ggH,ttH vs μ_VBF,VH model as an example, we can compare the two approaches by making 1d likelihood scans as function of μ_ggH,ttH keeping μ_VBF,VH fixed or profiled: we use the option -P to select the parameter to fit, and the option --floatOtherPOI to configure what to do with the others (μ_VBF,VH, in our case)

combine -M MultiDimFit rVrF_hgg.root -m 125.7 --algo=grid --points=50 --setPhysicsModelParameterRanges RF=-1,3 -P RF --floatOtherPOI=0 -n HGG_RVRF_SCAN_1D_RF_FIX     
combine -M MultiDimFit rVrF_hgg.root -m 125.7 --algo=grid --points=50 --setPhysicsModelParameterRanges RF=-1,3 -P RF --floatOtherPOI=1 -n HGG_RVRF_SCAN_1D_RF_PROF 

%CODE{"cpp"}% // root -l higgsCombineHGG_RVRF_SCAN_1D_RF_* /// Attaching file higgsCombineHGG_RVRF_SCAN_1D_RF_FIX.MultiDimFit.mH125.7.root as _file0... /// Attaching file higgsCombineHGG_RVRF_SCAN_1D_RF_PROF.MultiDimFit.mH125.7.root as _file1... _file0->cd(); limit->Draw("2*deltaNLL:RF"); TGraph *graphFix = gROOT->FindObject("Graph")->Clone(); graphFix->Sort();

_file1->cd(); limit->Draw("2*deltaNLL:RF"); TGraph *graphProf = gROOT->FindObject("Graph")->Clone(); graphProf->Sort();

graphProf->SetLineWidth(3); graphProf->Draw("AC");

graphFix->SetLineWidth(2); graphFix->SetLineStyle(9); graphFix->SetLineColor(206); graphFix->Draw("C SAME"); %ENDCODE%

You can see that the two give different results, both in terms of the preferred value and in terms of the uncertainties. The interpretation of the two results is also different from the physics point of view:

• the scan with μ_VBF,VH profiled implies that you assume that both parameters can deviate from the SM prediciton, but you are interested only in the information about μ_ggH,ttH
• the scan with μ_VBF,VH fixed instead implies that you believe that deviations from the SM are possible only in μ_ggH,ttH, not in μ_VBF,VH (e.g. you could take the stance that the coupling of the Higgs to vector bosons has to be necessarily like in the SM since it's needed for EW symmetry breaking and unitariety, while the coupling to the fermions is just an ad-hoc addition to the theory to explain also fermion masses with the same scalar field, not enforced by anything)

First fits to Higgs search results

• Separately for each decay mode, WW, ZZ, γγ, bb, and ττ, do a fit of the signal strength. Similarly, do a fit of the signal strength for the combination of the ttH datacards.
  Suggestion: each student can run one or two of them, and then you can put together the results.

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

bash mkcomb.sh

combineCards.py $TTH > comb_ttH.txt

text2workspace.py  -m 125.7  -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs datacard.txt   -o datacard.root
combine -M MaxLikelihoodFit --justFit --robustFit=1 -m 125.7 datacard.root

combine -M MaxLikelihoodFit --justFit --robustFit=1 -m 125.7 datacard.txt

<--/twistyPlugin-->

• Separately for each decay mode, WW, ZZ, γγ, bb, and ττ, do bidimensional scan of the likelihood as function of μ_ggH,ttH vs μ_VBF,VH
  Again, each student could do one of these and then you can put the results together.
  Some practical suggestions are:
  • You can use the grid3x3 algorithm instead of the grid one in MultiDimFit: this algorithm effectively replaces each point of the grid with a 3x3 set of points, but if the central point of the 3x3 set is far away from the interesting region it will compute the remaining 8 points with a slightly approximate formula; so a 20x20 grid becomes 60x60, but the time to fill it is much less than what you would need for a real 60x60 grid).
  • You can pass the option --for-fit to text2workspace.py; this will create only the signal+background pdf and not the background-only pdf, and so save some time

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

for X in h{bb,tt,gg,ww,zz}; do text2workspace.py comb_${X}.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.PhysicsModel:rVrFXSHiggs  -o rVrF_${X}.root --for-fit ; done
for X in h{bb,tt,gg,ww,zz}; do combine -M MultiDimFit rVrF_${X}.root -m 125.7 -n ${X}_RVRF_SCAN_2D --algo=grid3x3 --points=400 --setPhysicsModelParameterRanges RF=-1,3:RV=-2,5; done

To see a single plot, interactively in ROOT,

%CODE{"cpp"}% /// root -l tdrstyle.cc contours2D.cxx higgsCombinehww_RVRF_SCAN_2D.MultiDimFit.mH125.7.root // Processing tdrstyle.cc... // Processing contours2D.cxx... // Attaching file higgsCombinehww_RVRF_SCAN_2D.MultiDimFit.mH125.7.root as _file0... contour2D("RF",60,-1,3, "RV",60,-2,5, 1.0,1.0); %ENDCODE%

H → WW	H → ZZ	H → γγ	H → bb	H → ττ

A macro to make all the plots including the summary one is below:

%CODE{"cpp"}% // run with root.exe -b -l tdrstyle.cc contours2D.cxx plotall_rvrf.cxx -q void plotall_rvrf() { TCanvas *c1 = new TCanvas("c1","c1"); TFile *fOut = TFile::Open("rvrf_all.root", "RECREATE"); const int decays = 5; const char *decay[decays] = { "hww", "hzz", "hgg", "hbb", "htt" }; const int color[decays] = { 4 , 2 , 209 , 67 , 221 }; for (int i = 0; i < decays; ++i) { gFile = TFile::Open(Form("higgsCombine%s_RVRF_SCAN_2D.MultiDimFit.mH125.7.root", decay[i])); contour2D("RF",60,-1,3,"RV",60,-2,5,1.0,1.0,fOut,Form("rvrf_%s",decay[i])); c1->Print(Form("output_rvrf_%s.png",decay[i])); gFile->Close(); } fOut->cd(); TH2F *frame = new TH2F("frame",";#mu_{ggH,ttH};#mu_{VBF,VH}",60,-1.,3.,60,-2.,5.); frame->Draw("AXIS"); for (int i = 0; i < decays; ++i) { TGraph *best = fOut->Get(Form("rvrf_%s_best",decay[i])); best->SetMarkerColor(color[i]); best->SetMarkerStyle(34); // apparently ROOT forgets this best->SetMarkerSize(2.0); // so we do this again best->Draw("P SAME"); TList *c68 = fOut->Get(Form("rvrf_%s_c68",decay[i])); styleMultiGraph(c68, /*color=*/color[i], /*width=*/3, /*style=*/1); c68->Draw("L SAME"); } TMarker m; m.SetMarkerSize(3.0); m.SetMarkerColor(97); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); m.SetMarkerSize(1.8); m.SetMarkerColor(89); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); c1->Print("output_rvrf_all.png"); fOut->Close(); } %ENDCODE%

<--/twistyPlugin-->

Higgs Coupling models

The previous physics models were considering only some common signal strength modifiers aff Under the hypothesis that the boson observed is a narrow scalar state with the same structure of couplings as the SM Higgs boson, the Higgs production and decay modes decouple and we can parameterize the expected yields in terms of production cross sections and decay widths:
\sigma\cdot\mathrm{BR}(X \to \mathrm{H}\to Y) = \sigma_{X}\cdot\frac{\Gamma_{Y}}{\Gamma_\mathrm{tot}}
where the total width is given by the sum of the partial widths in all SM decay modes, plus possibly some BSM decays
\Gamma_\mathrm{tot} = \sum_{X} \Gamma_{X}\ +\ \Gamma_\mathrm{BSM}

To search for possible deviations in the data from the expectations for the SM Higgs boson, we introduce some coupling modifiers κ_X, and express the cross sections and partial widths as function of these κ_X (where κ_X=1 corresponds to the SM Higgs). We will consider different models, each characterized by a different set of parameters and by different assumptions.

The κ_V, κ_f benchmark model.

In this model, we introduce just two parameters κ_V, κ_f that scale the couplings between the Higgs boson and the electroweak gauge bosons and between the Higgs boson and the fermions.

The loop-induced effective couplings of the Higgs boson to gluons and to photons are computed in terms of κ_V, κ_f by assuming that only the SM particles participate to the loop.

• for the gluon coupling, at the leading order the process is induced only by loops containing quarks, so the coupling will scale as κ_f
• for the photon coupling, the W and quark loops both contribute, and interfere destructively in the SM. The effective coupling scales as |α κ_V + β κ_f|, with α ~ 1.28 and β ~ -0.28 (for a SM Higgs boson mass of 125.7 GeV)

In this model we assume that there are no BSM decays, so
\Gamma_\mathrm{tot} = \Gamma_\mathrm{ff} + \Gamma_\mathrm{gg} + \Gamma_\mathrm{VV} = |\kappa_\mathrm{f}|^2(\Gamma_\mathrm{ff}^\mathrm{SM} + \Gamma_\mathrm{gg}^\mathrm{SM}) + |\kappa_\mathrm{V}|^2\Gamma_\mathrm{VV}^\mathrm{SM}
where we have neglected the γγ and Zγ decays that contribute only to the few per-mille level. We can then dividing by the SM width, so that we can express a scale factor for the total width κ_H in terms of SM BRs
\kappa_\mathrm{H}^2 := \frac{\Gamma_\mathrm{tot}}{\Gamma_\mathrm{tot}^\mathrm{SM}} = \ |\kappa_\mathrm{f}|^2 (\mathrm{BR}_\mathrm{ff}^\mathrm{SM} + \mathrm{BR}_\mathrm{gg}^\mathrm{SM}) + |\kappa_\mathrm{V}|^2\mathrm{BR}_\mathrm{VV}^\mathrm{SM}
Since in the SM for a 125.7 GeV Higgs we have BR(H→VV) = 0.254, to a good approximation
\kappa_\mathrm{H}^2 = \frac14 |\kappa_\mathrm{V}|^2 + \frac34 |\kappa_\mathrm{f}|^2

Exercises for the reader

• Calculate the scale factors for σ·BR of ggH with H→ZZ; VBF with H→WW, VBF with H→γγ; VH with H→bb; ttH with H→bb. Express the results in terms of κ_V, κ_f, κ_H, α, β.

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

• For each of the processes before, calculate the first derivative of the scale factor vs κ_f² and κ_V² around the SM point (κ_f = κ_V = 1), using the numbers provided earlier for α, β and κ_H. What is the channel whose event yield depends more strongly on κ_V? and on κ_f.

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

Implementation in combine

The PhysicsModel that implements the κ_V and κ_f model is CvCfHiggs (when the model was first concieved, κ_V was called C_V and κ_f was called C_F). The relevant python code is below (for clarity, the parts of the code which allow the possibility of leaving the Higgs mass floating, and of setting the ranges for CV and CF, as well have been hidden). A few more comments have also been added here to make it more easy to understand.

%CODE{"python"}% class CvCfHiggs(SMLikeHiggsModel): def __init__(self): SMLikeHiggsModel.__init__(self) # not using 'super(x,self).__init__' since I don't understand it self.floatMass = False self.cVRange = ['0','2'] self.cFRange = ['-2','2']

def doParametersOfInterest(self): """Create POI out of signal strength and MH""" # --- Signal Strength as only POI --- self.modelBuilder.doVar("CV[1,%s,%s]" % (self.cVRange[0], self.cVRange[1])) self.modelBuilder.doVar("CF[1,%s,%s]" % (self.cFRange[0], self.cFRange[1]))

if self.floatMass: #[...] else: #[...] self.modelBuilder.doSet("POI",'CV,CF')

## Create an object "SMHiggsBuilder" that can compute the SM predictions for the Higgs ## (e.g. it encodes SM BRs and the the alpha and beta parameters for the diphoton decay) self.SMH = SMHiggsBuilder(self.modelBuilder) self.setup()

def setup(self): ## Define the four possible modes according to which a decay can scale: the two simple ones (HVV and Hff, which also works for the gluon decay) ## and the two complex ones that require to know what particles contribute in the SM loops (diphoton, Zgamma) self.decayScaling = 'hgg':'hgg', 'hzg':'hzg', 'hww':'hvv', 'hzz':'hvv', 'hbb':'hff', 'htt':'hff', }

## Define which coupling scales which production cross section (always trivial) self.productionScaling = { 'ggH':'CF', 'ttH':'CF', 'qqH':'CV', 'WH':'CV', 'ZH':'CV', 'VH':'CV', }

## Use tables provided by the theory community to compute the scaling for the diphoton ## and Zgamma decay as function of the various couplings self.SMH.makeScaling('hgg', Cb='CF', Ctop='CF', CW='CV', Ctau='CF') self.SMH.makeScaling('hzg', Cb='CF', Ctop='CF', CW='CV', Ctau='CF')

## partial widths, normalized to the SM one, for decays scaling with F, V and total for d in [ "htt", "hbb", "hcc", "hww", "hzz", "hgluglu", "htoptop", "hgg", "hzg", "hmm", "hss" ]: self.SMH.makeBR(d) self.modelBuilder.factory_('expr::CvCf_Gscal_sumf("@0*@0 * (@1+@2+@3+@4+@5+@6+@7)", CF, SM_BR_hbb, SM_BR_htt, SM_BR_hcc, SM_BR_htoptop, SM_BR_hgluglu, SM_BR_hmm, SM_BR_hss)') self.modelBuilder.factory_('expr::CvCf_Gscal_sumv("@0*@0 * (@1+@2)", CV, SM_BR_hww, SM_BR_hzz)') self.modelBuilder.factory_('expr::CvCf_Gscal_gg("@0 * @1", Scaling_hgg, SM_BR_hgg)') self.modelBuilder.factory_('expr::CvCf_Gscal_Zg("@0 * @1", Scaling_hzg, SM_BR_hzg)')

## Total width self.modelBuilder.factory_('sum::CvCf_Gscal_tot(CvCf _Gscal_sumf, CvCf _Gscal_sumv, CvCf _Gscal_gg, CvCf _Gscal_Zg)')

## BRs, normalized to the SM ones: they scale as (coupling/coupling_SM)^2 / (totWidth/totWidthSM) self.modelBuilder.factory_('expr::CvCf_BRscal_hgg("@0/@1", Scaling_hgg, CvCf _Gscal_tot)') self.modelBuilder.factory_('expr::CvCf_BRscal_hzg("@0/@1", Scaling_hzg, CvCf _Gscal_tot)') self.modelBuilder.factory_('expr::CvCf_BRscal_hff("@0*@0/@1", CF, CvCf _Gscal_tot)') self.modelBuilder.factory_('expr::CvCf_BRscal_hvv("@0*@0/@1", CV, CvCf _Gscal_tot)')

## Print out everything self.modelBuilder.out.Print()

def getHiggsSignalYieldScale(self,production,decay,energy): ## create a function that provides the scale factor for a given (production,decay) pair, unless it exists already

name = "CvCf_XSBRscal_%s_%s" % (production,decay) if self.modelBuilder.out.function(name): return name

XSscal = self.productionScaling[production] BRscal = self.decayScaling[decay] self.modelBuilder.factory_('expr::%s("@0*@0 * @1", %s, CvCf _BRscal_%s)' % (name, XSscal, BRscal)) return name

## This lines below are instead in HiggsCouplings.py from HiggsAnalysis.CombinedLimit.HiggsBenchmarkModels.VectorsAndFermionsModels import CvCfHiggs, CvCfXgHiggs, CfXgHiggs cVcF = CvCfHiggs()

%ENDCODE%

Working example: H→WW

Let's now investigate the constraints in the (κ_V,κ_f) provided by the H→WW analysis separately for 0 jets and 2 jets, and for the combination of the two

We can create the three workspaces:

text2workspace.py hwwof_0j_shape_8TeV.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_hww0j.root --for-fits
text2workspace.py hwwof_2j_shape_8TeV.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_hww2j.root --for-fits
text2workspace.py comb_hww.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_hww.root --for-fits

and then we run the scans on the Asimov dataset, since they are easier to interpret

for X in hww hww0j hww2j; do 
   ( combine -M MultiDimFit cVcF_${X}.root -m 125.7 -n ${X}_CVCF_SCAN_2D_ASIMOV --algo=grid3x3 --points=225 --setPhysicsModelParameterRanges CF=0,2:CV=0,2 -t -1 --expectSignal=1 > log.cvf$X 2>&1 & );
done

The resulting plot is the following:

Show code to make it Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

The shape of the H→WW + 2 jets can be explained considering that the expected signal is mostly from VBF. As computed in the previous questions from the reader, the expected scale factor for VBF H→WW is
\mu(\mathrm{VBF H}\to\mathrm{WW}) = \frac{\kappa_\mathrm{V}^2 \cdot \kappa_\mathrm{V}^2}{\frac14 \kappa_\mathrm{V}^2 + \frac34\kappa_\mathrm{f}^2 }
as long as κ_V ≤ 2κ_f, the κ_V at the denominator can be neglected, and so the measurement of μ is a measurement of κ_V²/κ_f: the contour will then have a parabolic shape κ_f ~ κ_V².

The shape of the H→WW + 0 jets is also understandable, considering that the signal is mostly gluon fusion: if κ_f is large enough, we can again neglect the κ_V at denominator, but for this channel now the κ_f at denominator cancels with the one in the numerator (from gluon fusion) and what is left is a measure of κ_V alone, giving an allowed region shaped as a vertical band at constant κ_V. However if κ_f is instead small with respect to κ_V, as it happens in the bottom-right part of the plot, it is now κ_V to dominate the denominator; this κ_V will cancel the one at the numerator (from the Γ_WW) leaving just a measure of κ_f (i.e. a horizontal band).

C6: a general six-parameter model for the couplings of a SM-like Higgs boson

The κ_V, κ_f model is quite restrictive, assuming e.g. the same coupling strength modifier for both quarks and leptons, and for both up-type and down-type fermions, assumptions which are not true in many BSM scenarios (e.g. in supersymmetry the couplings to the top and bottom quarks have different modifiers). We will consider therefore a more generic model which relies on a much smaller set of assumptions, namely:

• that the modifier to the W and Z couplings is the same (κ_V), as common also in most BSM scenarios.
• the coupling modifier to the charm quark is the same as the one to the top quark: the Higgs coupling to the charm is not measured directly at LHC, but some constraint on it has to be set as it indirectly affects the prediction through the total width. Relating it to the one of the top is motivated by the fact that most BSM scenarios alter the two in the same way since they're both up-type quarks. Similar assumptions are also made to identify the coupling modifier to the strange quark with the one of the bottom quark, and the muon to the tau, though these have even less importance since the BR of these decays is very small (for the muon we could rely also on the direct measurement of H→μμ, but we're not including that channel in this exercise)
• there are no BSM decays of the Higgs: the total width is taken as the sum of the partial widths of the SM-like decays (we will relax this assumption later in the final project of this part of the exercise)

In this model we therefore have six parameters:

Exercise for the reader

• Express the total width as function of the coupling modifiers, for a Higgs boson mass of 125.7 GeV. You can take the BRs from CERN Yellow Report 3, or from the data/lhc-hxswg/sm/br/BR2bosons.txt and data/lhc-hxswg/sm/br/BR2fermions.txt text files in the combination package.

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

<-- 

 Consider these three processes: ggH Ha†’WW, VBF Ha†’Î3Î3, VH Ha†’bb, and 
 
-->

Implementation in combine

The PhysicsModel that implements the c6 model is C6. The relevant python code is below (for clarity, the parts of the code which allow the possibility of leaving the Higgs mass floating, as well have been hidden, and a few more comments have also been added here to make it more easy to understand (and some misleading obsolete comments removed...).

%CODE{"python"}% class C6(SMLikeHiggsModel): def __init__(self): SMLikeHiggsModel.__init__(self) # not using 'super(x,self).__init__' since I don't understand it self.floatMass = False ## This model has the option to have a separate coupling for the HZgamma vertex, we won't use it in this exercise self.doHZg = False def setPhysicsOptions(self,physOptions): for po in physOptions: if po.startswith("higgsMassRange="): # [...] if po = 'doHZg': self.doHZg = True def doParametersOfInterest(self): """Create POI out of signal strength and MH""" self.modelBuilder.doVar("kV[1,0.0,2.0]") self.modelBuilder.doVar("ktau[1,0.0,2.0]") self.modelBuilder.doVar("ktop[1,0.0,2.0]") self.modelBuilder.doVar("kbottom[1,0.0,2.0]") self.modelBuilder.doVar("kgluon[1,0.0,2.0]") self.modelBuilder.doVar("kgamma[1,0.0,2.0]") pois = 'kV,ktau,ktop,kbottom,kgluon,kgamma' if self.doHZg: self.modelBuilder.doVar("kZgamma[1,0.0,30.0]") pois + ",kZgamma" if self.floatMass: # [...] else: # [...] self.modelBuilder.doSet("POI",pois)

## Create the SMHiggsBuilder, to get the SM BRs self.SMH = SMHiggsBuilder(self.modelBuilder) self.setup()

def setup(self):

# SM BRs, needed to compute the modified total width as funciton of the partial widths for d in [ "htt", "hbb", "hcc", "hww", "hzz", "hgluglu", "htoptop", "hgg", "hzg", "hmm", "hss" ]: self.SMH.makeBR(d)

## Partial decay witdhs, normalized to the SM one self.modelBuilder.factory_('expr::c6_Gscal_Vectors("@0*@0 * (@1+@2)", kV, SM_BR_hzz, SM_BR_hww)') self.modelBuilder.factory_('expr::c6_Gscal_tau("@0*@0 * (@1+@2)", ktau, SM_BR_htt, SM_BR_hmm)') self.modelBuilder.factory_('expr::c6_Gscal_top("@0*@0 * (@1+@2)", ktop, SM_BR_htoptop, SM_BR_hcc)') self.modelBuilder.factory_('expr::c6_Gscal_bottom("@0*@0 * (@1+@2)", kbottom, SM_BR_hbb, SM_BR_hss)') self.modelBuilder.factory_('expr::c6_Gscal_gluon("@0*@0 * @1", kgluon, SM_BR_hgluglu)') if not self.doHZg: self.modelBuilder.factory_('expr::c6_Gscal_gamma("@0*@0 * (@1+@2)", kgamma, SM_BR_hgg, SM_BR_hzg)') else: self.modelBuilder.factory_('expr::c6_Gscal_gamma("@0*@0 *@1+@2*@2*@3", kgamma, SM_BR_hgg, kZgamma, SM_BR_hzg)')

## Modifier for the total width, from the sum of the partial widths self.modelBuilder.factory_('sum::c6_Gscal_tot(c6_Gscal_Vectors, c6_Gscal_tau, c6_Gscal_top, c6_Gscal_bottom, c6_Gscal_gluon, c6_Gscal_gamma)')

## BRs, normalized to the SM ones: they scale as (coupling/coupling_SM)^2 / (total/total_SM) self.modelBuilder.factory_('expr::c6_BRscal_hvv("@0*@0/@1", kV, c6_Gscal_tot)') self.modelBuilder.factory_('expr::c6_BRscal_htt("@0*@0/@1", ktau, c6_Gscal_tot)') self.modelBuilder.factory_('expr::c6_BRscal_hbb("@0*@0/@1", kbottom, c6_Gscal_tot)') self.modelBuilder.factory_('expr::c6_BRscal_hgg("@0*@0/@1", kgamma, c6_Gscal_tot)') if self.doHZg: self.modelBuilder.factory_('expr::c6_BRscal_hzg("@0*@0/@1", kZgamma, c6_Gscal_tot)')

def getHiggsSignalYieldScale(self,production,decay,energy): name = "c6_XSBRscal_%s_%s" % (production,decay) print '[LOFullParametrization::C6]' print name, production, decay, energy if self.modelBuilder.out.function(name) = None: XSscal = "kgluon" if production in ["WH","ZH","VH","qqH"]: XSscal = "kV" if production = "ttH": XSscal = "ktop" BRscal = "hgg" if decay in ["hbb", "htt"]: BRscal = decay if decay in ["hww", "hzz"]: BRscal = "hvv" if self.doHZg and decay == "hzg": BRscal = "hzg" self.modelBuilder.factory_('expr::%s("@0*@0 * @1", %s, c6_BRscal_%s)' % (name, XSscal, BRscal)) return name

## This lines below are instead in HiggsCouplings.py from HiggsAnalysis.CombinedLimit.LOFullParametrization import C5, C6, C7, PartialWidthsModel c6 = C6() %ENDCODE%

Working example 1: κ_γ

We will take the full combination but excluding the H→inv search, and extract the κ_γ modifier, profiling the others. However, to profile them correctly we need to select ranges that are large enough so that the parameter doesn't end on the boundary. We can look at the past CMS Higgs combination, HIG-13-005, to see what was the sensitivity to the other parameters, e.g from the summary plot

In addition to the observed outcome, we will compute also the Asimov one

text2workspace.py comb.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.HiggsCouplings:c6 -o c6_comb.root --for-fit

RANGES="--setPhysicsModelParameterRanges kbottom=0,5:ktop=0,8:ktau=0,4:kgluon=0,4:kgamma=0,3"

combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma
/// note: 40 points take about 4 minutes (on lxplus)

combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV -t -1 --expectSignal=1
/// note: 40 points take about 20 minutes (on lxplus)

Speeding up things: running scans in parallel

When the model becomes complex, evaluating all the points of the scan in a single job can be inefficient. For this purpose, MultiDimFit allows one to fraction the task by specifying a range of points to run: different ranges can be run in parallel, and at the end the root files can be merged with hadd to get the final one. For example, for the Asimov scan of κ_γ we could do the 40 points in 5 jobs, doing

combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.0 -t -1 --expectSignal=1 --firstPoint 0 --lastPoint 7
combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.1 -t -1 --expectSignal=1 --firstPoint 8 --lastPoint 15
combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.2 -t -1 --expectSignal=1 --firstPoint 16 --lastPoint 23
combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.3 -t -1 --expectSignal=1 --firstPoint 24 --lastPoint 31
combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.4 -t -1 --expectSignal=1 --firstPoint 32 --lastPoint 39
hadd higgsCombineC6_comb_SCAN_1D_kgamma_ASIMOV.MultiDimFit.mH125.7.root higgsCombineC6_comb_SCAN_1D_kgamma_ASIMOV.*.MultiDimFit.mH125.7.root

where the first 5 commands can be executed in parallel.

For convenience, there is even already a script that takes care of running all the jobs in parallel splitting the points, suitable for running on a multi-core machine: test/parallelScan.py. For the previous case, we could for example do

cp -s ${CMSSW_BASE}/src/HiggsAnalysis/CombinedLimit/test/parallelScan.py .
parallelScan.py -j 5 c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV -t -1 --expectSignal=1
hadd higgsCombineC6_comb_SCAN_1D_kgamma_ASIMOV.MultiDimFit.mH125.7.root higgsCombineC6_comb_SCAN_1D_kgamma_ASIMOV.*.MultiDimFit.mH125.7.root

The option -j N is used to specify the number of jobs; if not present, it will use one job per CPU of the machine, or half that number if it detects that is being run on a lxplus or cmslpc node (one should be respectful of other users logged in the same machine)

Working example 2: κ_g vs κ_b degeneracy

The channels with the best sensitivity to the signal strength are mostly driven by the gluon fusion production cross section (e.g. H→WW and H→ZZ, H→γγboosted H→ττ). Conversely, the H→bb is the dominant decay mode, so that the total width is mostly sensitive to κ_b. The combination of these two features results in a strong correlation between the uncertainties of the two measurements of κ_g and κ_b: a simultaneous increase of both couplings would produce an increase in the gluon fusion production cross section but also a decrease of the branching ratio in the other modes (WW, ZZ, γγ and ττ), leaving the expected σ·BR unchanged, and so provides still a good fit to the data.

We can see this feature with a 2D likelihood scan:

./parallelScan.py c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgluon -P kbottom --floatOtherPOI=1 $RANGES --points=400 -n C6_comb_SCAN_2D_kgluon_kbottom

The contour in these two parameters is very elongated, exemplifying this correlation.

Coupling fits from the CMS combination

• Make 2D likelihood scans in the κ_V,κ_f plane for each of the five decay modes and for the full combination (but not including the H→invisible searches)

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

text2workspace.py comb.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_comb.root --for-fit ; 
for X in h{bb,tt,gg,ww,zz}; do text2workspace.py comb_${X}.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_${X}.root --for-fit ; done
for X in h{bb,tt,gg,ww,zz,comb}; do combine -M MultiDimFit cVcF_${X}.root -m 125.7 -n ${X}_CVCF_SCAN_2D --algo=grid3x3 --points=225 --setPhysicsModelParameterRanges CV=0,2:CF=0,2; done

<--/twistyPlugin-->

• Determine the best fit values for all the 6 couplings in the C6 model
  Note: if you want to get also the uncertainties on each parameter individually by running MultiDimFit with --algo=singles (this is the same as running MaxLikeilhoodFit but applies also to models with more parameters)

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

combine c6_comb.root -M MultiDimFit -m 125.7 $RANGES

>>> including systematics
>>> method used to compute upper limit is MultiDimFit
>>> random number generator seed is 123456
ModelConfig 'ModelConfig' defines more than one parameter of interest. This is not supported in some statistical methods.
Set Range of Parameter kbottom To : (0,5)
Set Range of Parameter ktop To : (0,8)
Set Range of Parameter ktau To : (0,4)
Set Range of Parameter kgluon To : (0,4)
Set Range of Parameter kgamma To : (0,3)
Computing limit starting from observation

 --- MultiDimFit ---
best fit parameter values:
        kV :    +1.044
      ktau :    +1.274
      ktop :    +2.053
   kbottom :    +1.401
    kgluon :    +1.039
    kgamma :    +1.143
Done in 0.11 min (cpu), 0.11 min (real)

combine c6_comb.root -M MultiDimFit --algo=singles -m 125.7 $RANGES

>>> including systematics
>>> method used to compute upper limit is MultiDimFit
>>> random number generator seed is 123456
ModelConfig 'ModelConfig' defines more than one parameter of interest. This is not supported in some statistical methods.
Set Range of Parameter kbottom To : (0,5)
Set Range of Parameter ktop To : (0,8)
Set Range of Parameter ktau To : (0,4)
Set Range of Parameter kgluon To : (0,4)
Set Range of Parameter kgamma To : (0,3)
Computing limit starting from observation
Running Minos for POI
Real time 0:00:11, CP time 11.700

Running Minos for POI
Real time 0:00:11, CP time 11.700

Running Minos for POI
Real time 0:00:11, CP time 11.700

Running Minos for POI
Real time 0:00:11, CP time 11.700

Running Minos for POI
Real time 0:00:11, CP time 11.700

Running Minos for POI
Real time 0:00:11, CP time 11.700


 --- MultiDimFit ---
best fit parameter values and profile-likelihood uncertainties:
        kV :    +1.044   -0.234/+0.204 (68%)
      ktau :    +1.274   -0.497/+0.504 (68%)
      ktop :    +2.053   -0.897/+1.702 (68%)
   kbottom :    +1.401   -0.674/+1.227 (68%)
    kgluon :    +1.039   -0.324/+0.811 (68%)
    kgamma :    +1.143   -0.339/+0.388 (68%)
Done in 1.55 min (cpu), 1.64 min (real)

<--/twistyPlugin-->

Going beyond: introducing H→invis decays into C6

The C6 model assumes that there are no BSM decays of the Higgs. We will now take one step forward, and allow the possibility of Higgs decays into invisible particles, relying on the searches in the VBF mode to constrain this.

We will introduce a new parameter BR_Inv, representing the BR of the Higgs boson decays to invisible particles.

• How do the BRs depend on the κs and BR_Inv ?

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

• Define a C6Inv PhysicsModel in combine copying from C6 but introducing the H→inv decay

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

• Using this C6Inv model and taking only the VBF H→inv datacard, do a 2D likelihood scan of κ_V vs BR_Inv, keeping all the other parameters fixed to the SM values.
• Same as above, but for the combination of the H→WW datacards.
• Same as above, but for the combination of the VBF H→inv and H→WW datacards.

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

text2workspace.py hinv_vbf_8TeV.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_vbfinv.root --for-fit
combine -m 125.7 -M MultiDimFit c6i_vbfinv.root -n C6I_vbfinv_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=0 --points=400 --algo=grid3x3 $RANGES

text2workspace.py comb_hww.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_hww.root --for-fit
./parallelScan.py -M MultiDimFit c6i_hww.root -n C6I_hww_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=0 --points=100 --algo=grid3x3 $RANGES
hadd -f higgsCombineC6I_hww_SCAN_2D_kV_BRInv.MultiDimFit.mH125.7.root higgsCombineC6I_hww_SCAN_2D_kV_BRInv.*.MultiDimFit.mH125.7.root

HWW="hww0j_8TeV=hwwof_0j_shape_8TeV.txt hww2j_8TeV=hwwof_2j_shape_8TeV.txt"
HIN="hinv_8TeV=hinv_vbf_8TeV.txt"
combineCards.py -S $HWW $HIN > comb_test2.txt
text2workspace.py comb_test2.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_test2.root --for-fit
./parallelScan.py  -m 125.7  -M MultiDimFit c6i_test2.root -n C6I_test2_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=0 --points=100 --algo=grid3x3 $RANGES
hadd -f higgsCombineC6I_test2_SCAN_2D_kV_BRInv.MultiDimFit.mH125.7.root higgsCombineC6I_test2_SCAN_2D_kV_BRInv.*.MultiDimFit.mH125.7.root

VBF H→inv	H→WW	Combined

<--/twistyPlugin-->

• Using this C6Inv model and taking the combination of all datacards do a 2D likelihood scan of κ_V vs BR_Inv, profiling all other nuisance parameters; do this twice, including and not including the VBF H→Inv search in the combination.

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

text2workspace.py comb.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_comb.root --for-fit; 
text2workspace.py combx.txt -m 125.7  -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_combx.root --for-fit; 
./parallelScan.py  -m 125.7  -M MultiDimFit c6i_comb.root  -n C6I_comb_SCAN_2D_kV_BRInv  -P kV -P BRInv --floatOtherPOI=1 --points=400 --algo=grid $RANGES --hadd
./parallelScan.py  -m 125.7  -M MultiDimFit c6i_combx.root -n C6I_combx_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=1 --points=400 --algo=grid $RANGES ---hadd

<--/twistyPlugin-->

• Using this C6Inv model and taking the combination of all datacards, including the VBF H→Inv search, get the fitted values of all the 7 parameters, and their 1D uncertainties. In addition to the observed ones, compute also the expected uncertainties under the SM Higgs hypothesis.
  If you have time, do also some plots of 2D likeilhood scans that you think could be interesting

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

Best fit and one-dimensional uncertainties (observed): takes O(1h)

combine  -m 125.7  -M MultiDimFit c6i_combx.root -n C6I_comb_SINGLES --algo=singles $RANGES --robustFit=1

 --- MultiDimFit ---
best fit parameter values and profile-likelihood uncertainties:
        kV :    +1.146   -0.207/+0.199 (68%)
      ktau :    +1.399   -0.466/+0.428 (68%)
      ktop :    +2.000   -0.683/+0.000 (68%)
   kbottom :    +1.428   -0.537/+0.572 (68%)
    kgluon :    +1.075   -0.265/+0.407 (68%)
    kgamma :    +1.231   -0.306/+0.345 (68%)
     BRInv :    +0.177   -0.177/+0.157 (68%)
Done in 54.14 min (cpu), 55.02 min (real)

One-dimensional expected uncertainties: one parameter at a time in parallel, since it takes O(1h) per parameter.

for D in k{V,gamma,gluon,bottom,top,tau} BRInv; do 
     ( combine  -m 125.7  -M MultiDimFit c6i_combx.root -n C6I_comb_SINGLES_${D}_Asimov --algo=singles $RANGES --robustFit=1 -t -1 -P $D --floatOtherPOI=1 --expectSignal=1 2>&1  | tee log.C6I_Asimov_$D > /dev/null & ); 
done
tail -n 2 log.C6I_Asimov_*

==> log.C6I_Asimov_BRInv <==
   BRInv :    +0.039   -0.039/+0.190 (68%)
Done in 58.08 min (cpu), 65.68 min (real)

==> log.C6I_Asimov_kV <==
   kV :    +1.020   -0.201/+0.196 (68%)
Done in 62.85 min (cpu), 70.62 min (real)

==> log.C6I_Asimov_kbottom <==
   kbottom :    +1.020   -0.433/+0.585 (68%)
Done in 66.84 min (cpu), 74.38 min (real)

==> log.C6I_Asimov_kgamma <==
   kgamma :    +1.021   -0.225/+0.255 (68%)
Done in 61.08 min (cpu), 68.59 min (real)

==> log.C6I_Asimov_kgluon <==
   kgluon :    +1.017   -0.241/+0.457 (68%)
Done in 74.38 min (cpu), 81.99 min (real)

==> log.C6I_Asimov_ktau <==
   ktau :    +1.022   -0.366/+0.323 (68%)
Done in 74.65 min (cpu), 82.36 min (real)

==> log.C6I_Asimov_ktop <==
   ktop :    +1.019   -1.019/+0.750 (68%)
Done in 57.58 min (cpu), 65.44 min (real)

One-dimensional likelihood scans (observed and expected)

for D in k{V,gamma,gluon,bottom,top,tau} BRInv; do 
    ./parallelScan.py -m 125.7 -M MultiDimFit c6i_combx.root -n C6I_combx_SCAN_1D_$D          -P ${D} --floatOtherPOI=1 $RANGES --algo=grid --points 50 --hadd;
    ./parallelScan.py -m 125.7 -M MultiDimFit c6i_combx.root -n C6I_combx_SCAN_1D_${D}_ASIMOV -P ${D} --floatOtherPOI=1 $RANGES --algo=grid --points 30 -t -1 --expectSignal=1 --hadd; 
done

Results:

κV	κg	κγ	BRInv
κt	κb	κτ

Plotting macro Hide

<--/twistyPlugin twikiMakeVisibleInline-->

%CODE{"cpp"}% void plot1ds_c6i() { TCanvas c1 = new TCanvas("c1","c1"); gStyle->SetOptStat(0); TFile *fOut = TFile::Open("plot1d_c6i.root", "RECREATE"); const int params = 7; const char *param[params] = { "kV", "kbottom", "ktop", "ktau", "kgluon", "kgamma", "BRInv" }; const float min[params] = { 0, 0, 0, 0, 0, 0, 0, }; const float max[params] = { 2, 5, 8, 4, 4, 3, 1, }; const char *label[params] = { "#kappa_{V}", "#kappa_{b}", "#kappa_{t}", "#kappa_{#tau}", "#kappa_{g}", "#kappa_{#gamma}","BR_{Inv}" }; for (int i = 0; i < params; ++i) { cout << param[i] << endl; gFile = TFile::Open(Form("higgsCombineC6I_combx_SCAN_1D_%s.MultiDimFit.mH125.7.root", param[i])); TTree *limit = (TTree) gFile->Get("limit"); limit->Draw(Form("2*deltaNLL:%s",param[i]),""); TGraph gO = gROOT->FindObject("Graph")->Clone(); gO->Sort(); gFile->Close(); gFile = TFile::Open(Form("higgsCombineC6I_combx_SCAN_1D_%s_ASIMOV.MultiDimFit.mH125.7.root", param[i])); limit = (TTree) gFile->Get("limit"); limit->Draw(Form("2*deltaNLL:%s",param[i]),""); TGraph *gA = gROOT->FindObject("Graph")->Clone(); gA->Sort(); gFile->Close();

gO->SetName(Form("c6i_1d_%s_obs", param[i])); gA->SetName(Form("c6i_1d_%s_asimov",param[i])); fOut->WriteTObject(gO); fOut->WriteTObject(gA);

TH1F *frame = new TH1F("frame",Form(";%s;- 2 #Delta ln L",label[i]), 10, min[i], max[i]); frame->Draw(); frame->GetYaxis()->SetRangeUser(0,7); gO->SetLineWidth(3); gO->SetLineColor(1); gA->SetLineWidth(2); gA->SetLineColor(1); gA->SetLineStyle(2); gO->Draw("C SAME"); gA->Draw("C SAME");

c1->Print(Form("output_C6I_1D_%s.png", param[i])); delete frame; } fOut->Close(); } %ENDCODE%

<--/twistyPlugin-->

<--/twistyPlugin-->

Higgs J^CP

Introduction

The narrow resonance with an invariant mass of approximately 125 GeV that has been discovered by ATLAS and CMS, are mainly via diboson decay modes γγ, WW and ZZ. The determination of the quantum numbers of this Higgs-like particle, such as its spin and parity, is currently one of the most important tasks in high energy physics and thus is important for future studies.

Kinematics

The scattering amplitudes of Higgs or any exotic boson decays to the WW final state have been performed in several literatures. It was demonstrated that the WW chaneel has a good sensitivity to distinguish between the standard model Higgs boson hypothesis and a spin-2 resonance, which couples to the dibosons through minimal couplings, referred to as 2+. Some sensitivity in this final state can also distinguish between the min0+ and the pseudoscalar 0− boson hypotheses. The typical experimental observables are m(ll), m_T(ll,E_T^miss), and Δφ(ll).

The tdefinition of transverse mass m_T(ll,E_T^miss) used in this analysis is
m_T(ll,E_T^miss) := sqrt( 2 * p_T(ll) * E_T^miss * (1 - cos Δ#phi;(ll, E_T^miss) )

m(ll)	Δφ(ll)	mT(ll,ETmiss)

SM Higgs	Graviton

Exercise for the reader

Two datasets of simulated X→WW→2l2ν events have been produced, with X being a SM Higgs and spin 2 graviton, and a very basic event selection has already been applied to retain only the events with at least two high-pt leptons, and no hadronic jets.

What you can do is to run on those simulated events and make plots of the kinematic variables like the ones above (or even others, e.g. ΔR(ll)).

You can do this exercise either starting from simple ROOT trees, or from PAT files (for the latter, use CMSSW_5_3_14), using any of the analysis tools you should have learned from the pre-exercises (root, fwlite, python, ...).

• The files for the trees are HWW_Tree_SMHiggs.root and HWW_Tree_Graviton.root.
  Each file contains a TTree called tree/probe_tree which has Float_t branches for the 4-vectors of the two leading leptons ( lep1_pt, lep1_eta, lep1_phi, lep1_mass and the same for lep2) and the missing energy ( met_pt, met_phi).
• The PAT files are HWW_Tree_SMHiggs.root and HWW_Tree_Graviton.root.
  Each file contains the standard PAT collections of cleatPatMuons ( vector<pat::Muon>), cleanPatElectrons ( vector<pat::Electron>), patMETs ( vector<pat::MET>).
  The leptons have the selection already applied.

• /SMHiggsToWWTo2L2Nu_M-125_8TeV-jhu-pythia6/Summer12_DR53X-PU_S10_START53_V7A-v1/AODSIM
• /Graviton2PMToWWTo2L2Nu_M-125_8TeV-jhu-pythia6/Summer12_DR53X-PU_S10_START53_V7A-v1/AODSIM

Answers (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

Statistics

The statistics of hypotheses testing is based on the likelihood ratio between the hypotheses being tested.

Likelihood ratio

Consider the null (N) and alternate (A) hypotheses which correspond in this case to the SMH-like 0⁺+background (N) and the alternative J^CP+background (A). The likelihood ratio is defined as q = -2 ln( L(x|A) / L(x|N) ), where x represents a set of events. This set of events, x, can be:

• (the) observed data: leading to the observed value qobs.
• pseudo-data (toys): sampled from the expectation from a given hypothesis, like N or A.

You can find the example of such distributions below for the ZZ+WW combination performed in HIG-13-005:

Measures of separation and the CLs criterium

Given the distributions above, one can then characterise the expectations and observation in a number of ways:

• Separation, defined as 1-A_tail, where A_tail is the tail area calculated at the value of q' such that P( q < q' | N ) = P( q > q' | A ). The smaller A_tail, the more separated the two distributions are.
• The probability for one hypothesis to fluctuate beyond the median of the other hypothesis: P( q > qexp(N) | A ) and P( q < qexp(A) | N ), where qexp(X) stands for the median expected value of q under the hypothesis X as estimated from the distribution P(q|X).
• The CLs value.

• CLs(obs) = P( q > qobs | A ) / P( q > qobs | N ) for the observation.
• CLs(exp) = P( q > qexp(N) | A ) / 0.5, for the expectation. (By definition P( q > qexp(N) | N ) = 0.5.)

Pre-fit and post-fit expectations

A final complication relates to the normalisation of the processes entering the N and A hypotheses, especially the signal normalisation. We consider the "pre-fit" and "post-fit" cases:

• pre-fit: what you expect without looking at the data.
• post-fit what you expect after matching the signal yield to the best-fit signal yield in data.

the ZZ+WW combination performed in HIG-13-005

It can be seen that the when going from left to right, the distributions P(q|A) and P(q|N) become closer, less discernible, and the separation is smaller. This is due to the fact that the WW channel has a deficit of signal.

Finally, in the following table you can find the results of the hypothesis test performed for HIG-13-005, including all the quantities discussed above:

Exercise for the reader

What does it mean to show -0.9σ in the Table above?
Hint: how do you describe the case in which P( q < qobs | N ) < 0.5?

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

<--/twistyPlugin-->

Practice

For the spin-2 (2⁺_m in this case) vs spin-0 (0⁺) hypotheses test, the full H→WW→2a„“2ν analysis includes 4 channels corresponding to {0-jet,1-jet}⊗{7 TeV, 8 TeV}. For simplicity we will perform the exercise on the 0-jet 8 TeV data card.

head jcp_hwwof_0j_8TeV.txt

[...]
Observation 5747
shapes *   *   jcp_hwwof_0j.input_8TeV.root  histo_$PROCESS histo_$PROCESS_$SYSTEMATIC
shapes data_obs * jcp_hwwof_0j.input_8TeV.root  histo_Data
bin j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev
process qqbarH_ALT ggH_ALT ZH WH qqH ggH qqWW ggWW VV Top Zjets WjetsE Wgamma Wg3l Ztt WjetsM qqWW2j
process -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11
rate  129.869  146.427   1.728   5.833   3.033  234.884  3973.542  210.835  131.912  499.856   0.000  279.465  115.142  166.923  46.430  327.212   0.288

You will note that besides the usual SM signal contributions (ggH, qqH, ZH, and WH) there are two new processes (qqbarH_ALT and ggH_ALT) corresponding to the production of 2⁺_m via quark-antiquark annihilation and gluon-fusion.

Let's look at the shape of the templates:

root -l jcp_hwwof_0j.input_8TeV.root

%CODE{"cpp"}% //Compare main signal and background histo_ggH->GetXaxis()->SetRangeUser(-1,1) histo_ggH->SetLineColor(kBlack) histo_ggH->DrawNormalized("H") histo_qqWW->SetLineColor(kGray) histo_qqWW->DrawNormalized("Hsame") %ENDCODE%

You will notice multiple peaks. These correspond to an unrolled 2D histogram encoding the information of mT and mll as explained in https://twiki.cern.ch/twiki/bin/view/CMSPublic/Hig12042TWiki Basically, inside each peak mll is increasing, while different peaks correspond to larger mT values.

You can compare the plots you obtained with the following:

Let's now compare the different signal model shapes:

root -l jcp_hwwof_0j.input_8TeV.root

%CODE{"cpp"}% //Define a region where it is easy to see the differences between the alternate signals histo_ggH->;GetXaxis()->SetRangeUser(-1,0) histo_ggH->SetLineColor(kBlack) histo_ggH->DrawNormalized("H") histo_ggH_ALT->SetLineColor(kBlue) histo_ggH_ALT->DrawNormalized("Hsame") histo_qqbarH_ALT->SetLineColor(kRed) histo_qqbarH_ALT->DrawNormalized("Hsame") %ENDCODE%

You can see that the shapes are different, with different hypotheses having more or less signal in the different "peaks" (which correspond to different mT values), as shown in:

Building a nested model

The hypothesis test is performed by constructing a nested model which has two relevant parameters:

• fqq, the fraction qqbarH_ALT / (qqbarH_ALT + ggH_ALT)
• x, the fraction 0⁺ / (0⁺ + 2⁺_m)

You can immediately see that this general (nested) signal model encodes several interesting limiting cases:

The PhysicsModel that encodes the signal model above is the twoHypothesisHiggs.

Let's create the binary workspace. Suppose we blindly run:

text2workspace.py jcp_hwwof_0j_8TeV.txt -P HiggsAnalysis.CombinedLimit.HiggsJPC:twoHypothesisHiggs -m 125.7 --PO verbose -o jcp_hww.root

MH (not there before) will be assumed to be 125.7
Process  ggH_ALT  will get norm  x
Process  ZH  will get norm  not_x
Process  qqH  will get norm  not_x
Process  qqbarH_ALT  will get norm  x
Process  WH  will get norm  not_x
Process  ggH  will get norm  not_x
Process  ggH_ALT  will get norm  x
Process  ZH  will get norm  not_x
Process  qqH  will get norm  not_x
Process  qqbarH_ALT  will get norm  x
Process  WH  will get norm  not_x
Process  ggH  will get norm  not_x

See how the ggH_ALT and qqbarH_ALT are both being scaled with x? That is not what we want, since we would like each of the processes to be differentiated via fqq.

To do this properly we need to tell the model to also include fqq:

text2workspace.py jcp_hwwof_0j_8TeV.txt -P HiggsAnalysis.CombinedLimit.HiggsJPC:twoHypothesisHiggs -m 125.7 --PO verbose --PO fqqIncluded -o jcp_hww.root

Will consider fqq = fraction of qqH in Alt signal (signal strength will be left floating)
MH (not there before) will be assumed to be 125.7
Process  ggH_ALT  will get scaled by  r_times_x_times_not_fqq
Process  ZH  will get scaled by  r_times_not_x
Process  qqH  will get scaled by  r_times_not_x
Process  qqbarH_ALT  will get scaled by  r_times_x_times_fqq
Process  WH  will get scaled by  r_times_not_x
Process  ggH  will get scaled by  r_times_not_x
Process  ggH_ALT  will get scaled by  r_times_x_times_not_fqq
Process  ZH  will get scaled by  r_times_not_x
Process  qqH  will get scaled by  r_times_not_x
Process  qqbarH_ALT  will get scaled by  r_times_x_times_fqq
Process  WH  will get scaled by  r_times_not_x
Process  ggH  will get scaled by  r_times_not_x

In this case you see how ggH_ALT is scaled with r_times_x_times_not_fqq, i.e. r*x*(1-fqq), and qqbarH_ALT is scaled with r_times_x_times_fqq, i.e. r*x*fqq.

Running toys

Let's dive immediately into throwing toys for the three cases of interest:

• pre-fit expected.
• post-fit expected.
• post-fit observed.

GEN_OPTS="-m 125.7 jcp_hww.root --seed 8192"
WHAT_OPTS="--testStat=TEV --singlePoint 1 --saveHybridResult"
HOW_OPTS="-M HybridNew -T 500 -i 2 --fork 6 --clsAcc 0 --fullBToys"
combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --generateExt=1 --generateNuis=0  --expectedFromGrid 0.5 -n jcp_hww_pre-fit_exp  &> jcp_hww_pre-fit_exp.log
combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --frequentist                     --expectedFromGrid 0.5 -n jcp_hww_post-fit_exp &> jcp_hww_post-fit_exp.log
combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --frequentist                     -n jcp_hww_post-fit_obs &> jcp_hww_post-fit_obs.log
grep -H "^CLs" jcp_hww_{pre,post}*.log

While the commands are running (they'll take a few minutes), let's decode the command lines (taken from "combine --help").

There is a first set of common parameters which relate to what we are doing:

• --testStat=TEV "Test statistics: LEP, TEV, LHC (previously known as Atlas), Profile.": The TEV test-statistic definition matches the likelihood ratio that we want to use.
• --singlePoint 1 "Just compute CLs for the given value of r": Since we are doing an hypothesis test, it is enough to compute CLs and not do any scan/fit of signal strength.
• --saveHybridResult "Save result in the output file": Save the test-statistic values in a RooStats::HypoTestResult object that we can later use for plotting.

• -M HybridNew "Method to extract upper limit. Supported methods are: [...]": This method is used to throw toys. The name refers to the hybrid frequentist-bayesian treatment of toys. We'll configure it so as to be frequentist.
• -T 500 "Number of Toy MC extractions to compute CLs+b, CLb and CLs": 500 toys per iteration.
• -i 2 "Number of times to throw 'toysH' toys to compute the p-values (for --singlePoint if clsAcc is set to zero disabling adaptive generation)": 2 iterations of 500 toys, so 1000 toys in total.
• --fork=6 "Fork to N processes before running the toys (set to 0 for debugging)": Profit from multicore processors.
• --clsAcc 0 "Absolute accuracy on CLs to reach to terminate the scan": Turn off the CLs accuracy termination condition, since we are throwing a fixed number of toys.
• --fulBToys "Run as many B toys as S ones (default is to run 1/4 of b-only toys)": combine only generally knows about the background and the signal+background hypotheses. This ensures that we get as many null hypothesis as alternative hypothesis toys.

• --expectedFromGrid 0.5 "Use the grid to compute the expected limit for this quantile": Instead of the likelihood ratio in data (qobs) this will use the median from toys (qexp) when calculating CLs.
• --frequentist "Shortcut to switch to Frequentist mode (--generateNuisances=0 --generateExternalMeasurements=1 --fitNuisances=1)": This is a shortcut, so let's look at the individual options.
  • --generateNuisances=0 "Generate nuisance parameters for each toy": Nuisances stay fixed in the frequentist approach. In the hybrid approach, they are sampled randomly.
  • --generateExternalMeasurements=1 "Generate external measurements for each toy, taken from the GlobalObservables of the ModelConfig": Observables are sampled randomly as part of the frequentist procedure.
  • --fitNuisances=1 "Fit the nuisances, when not generating them.": This is the bit that makes these toys the post-fit ones. For pre-fit, this must be turned off.

jcp_hww_pre-fit_exp.log:CLs = 0.0582329 +/- 0.0108136
jcp_hww_post-fit_exp.log:CLs = 0.206827 +/- 0.0203793
jcp_hww_post-fit_obs.log:CLs = 0.089172 +/- 0.0171273

You can see two things:

• the post-fit expected CLs is 3 to 4 times larger than the pre-fit expected.
• the post-fit observed CLs is 2 times smaller than the post-fit expected.

The difference between the post-fit expected and the pre-fit expected is related to fact that the data does not accommodate a large signal.

We can check this also by performing fits for the signal strength for the two hypotheses: we create a workspace that has the signal strength as an additional parameter of interest, and then run fits for the signal strength r at a fixed value of the mixing x corresponding to the two hypotheses

text2workspace.py jcp_hwwof_0j_8TeV.txt -P HiggsAnalysis.CombinedLimit.HiggsJPC:twoHypothesisHiggs -m 125.7 --PO verbose --PO fqqIncluded -o jcp_hww_2poi.root --PO muAsPOI
combine jcp_hww_2poi.root -M MultiDimFit -m 125.7 -P r --floatOtherPOI=0 --setPhysicsModelParameters x=0 -n MU_SCALAR

 --- MultiDimFit ---
best fit parameter values:
   r :    +0.673

combine jcp_hww_2poi.root -M MultiDimFit -m 125.7 -P r --floatOtherPOI=0 --setPhysicsModelParameters x=1 -n MU_SPIN2

best fit parameter values:
   r :    +0.988

This implies that our post-fit expectation for the signal yield is less than the SM Higgs prediction, especially for the spin 0 case, and consequently also our post-fit expectation for the separation between the two hypotheses is worse (with less events, it's harder to separate them).

While the approach with toys and CLs is statistically more sound (or at least more in line with the conventions used in particle physics), we can see how well the two models fit the data also by looking at the likelihood as function of x, while profiling r.

combine jcp_hww_2poi.root -M MultiDimFit -m 125.7 -P x --floatOtherPOI=1 --algo=grid -n JCP_X_SCAN --points=200
combine jcp_hww_2poi.root -M MultiDimFit -m 125.7 -P x --floatOtherPOI=1 --algo=grid -n JCP_X_SCAN_Asimov --points=200 -t -1 --expectSignal=1 --setPhysicsModelParameters x=0,r=1

From there you can see that the data prefers the scalar hypothesis (x=0), and in particular the events are even more scalar-like than what expected on average for a scalar (the likelihood approaches zero with a non-zero slope). You can also see that for larger x the a-priori expected curve grows faster for increasing x (i.e. it has a smaller that the a-priori sensitivity is better) but that is in part compensated by the fact that the observed starts faster. For a gaussian likelihood, a difference in 2*Δ log L of about 3 would correspond to a p-value of about 10%, which is in the same ballpark as the CLs number; a precise comparison cannot be made, since the two kinds of tests are different and since we're considering a variable at the boundary of the interval in which it is defined (and so it is not really expected to behave as a gaussian)

In all the results above, f_qq is not floated: it is a constant parameter set at zero. This means that everything we have done up to now tests gg→2⁺_m vs 0⁺.

Exercise for the reader

Redo the same tests for the fqq=1 case?
Hint: look at "combine --help" and search for "PhysicsModelParameter".

Answer (don't look until you've tried it!) Hide

<--/twistyPlugin twikiMakeVisibleInline-->

combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --generateExt=1 --generateNuis=0  --expectedFromGrid 0.5 -n jcp_hww_pre-fit_exp1  --setPhysicsModelParameters fqq=1 &> jcp_hww_pre-fit_exp1.log
combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --frequentist                     --expectedFromGrid 0.5 -n jcp_hww_post-fit_exp1 --setPhysicsModelParameters fqq=1  &> jcp_hww_post-fit_exp1.log
combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --frequentist                     -n jcp_hww_post-fit_obs --setPhysicsModelParameters fqq=1  &> jcp_hww_post-fit_obs1.log

<--/twistyPlugin-->

Plotting

For simplicity, let's work with the results of the post-fit observed run above which are saved in higgsCombinejcp_hww_post-fit_obs.HybridNew.mH125.7.8192.root.

To do this we need to use a macro that will extract the test-statistic values from the HypoTestResult and save them in a ROOT TTree:

root -l -q -b higgsCombinejcp_hww_post-fit_obs.HybridNew.mH125.7.8192.root "${CMSSW_BASE}/src/HiggsAnalysis/CombinedLimit/test/plotting/hypoTestResultTree.cxx(\"jcp_hww_post-fit_obs.qvals.root\",125.7,1,\"x\")" &> jcp_hww_post-fit_obs.qvals.log
cat jcp_hww_post-fit_obs.qvals.log

Attaching file higgsCombinejcp_hww_post-fit_obs.HybridNew.mH125.7.8192.root as _file0...
Processing /Users/adavid/workspace/CMSSW611-master/src/HiggsAnalysis/CombinedLimit/test/plotting/hypoTestResultTree.cxx("jcp_hww_post-fit_obs.qvals.root",125.7,1,"x")...
 - HypoTestResult_mh125.7_x1_2618602738
Saved test statistics distributions for 996 signal toys and 996 background toys to jcp_hww_post-fit_obs.qvals.root.

You can see that the macro found 996 toys in total, which correctly matches the parallel work of 6 processors (166*6=996 while 167*6=1002). Open the output file jcp_hww_post-fit_obs.qvals.root:

root -l jcp_hww_post-fit_obs.qvals.root

and let's explore:

%CODE{"cpp"}% q->Print() %ENDCODE%

******************************************************************************
*Tree    :q         : Test statistics                                        *
*Entries :     1993 : Total =           42889 bytes  File  Size =       8799 *
*        :          : Tree compression factor =   4.95                       *
******************************************************************************
*Br    0 :q         : q/F                                                    *
*Entries :     1993 : Total  Size=       8495 bytes  File Size  =       7539 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=   1.07     *
*............................................................................*
*Br    1 :mh        : mh/F                                                   *
*Entries :     1993 : Total  Size=       8500 bytes  File Size  =        143 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=  56.21     *
*............................................................................*
*Br    2 :weight    : weight/F                                               *
*Entries :     1993 : Total  Size=       8520 bytes  File Size  =        146 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=  55.08     *
*............................................................................*
*Br    3 :type      : type/I                                                 *
*Entries :     1993 : Total  Size=       8510 bytes  File Size  =        145 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=  55.45     *
*............................................................................*
*Br    4 :r         : r/F                                                    *
*Entries :     1993 : Total  Size=       8495 bytes  File Size  =        141 *
*Baskets :        1 : Basket Size=      32000 bytes  Compression=  57.00     *
*............................................................................*

The type encode which hypothesis the value of q corresponds to (observation has type zero). To see that, run:

%CODE{"cpp"}% q->Draw("type") gPad->SetLogy() %ENDCODE%

And finally, you can quickly visualise the results of the hypothesis test by running:

%CODE{"cpp"}% gPad->SetLogy(kFalse) q->Draw("-2*q>>hnull(50,-10,10)","type>0","goff") q->Draw("-2*q>>halt(50,-10,10)", "type<0","goff") q->Draw("-2*q>>hobs(50,-10,10)", "type==0","goff") hnull->SetLineColor(kBlue) halt->SetLineColor(kOrange) hobs->SetLineColor(kBlack) hobs->Scale(1e3) hobs->SetFillColor(kBlack) hnull->Draw() halt->Draw("same") hobs->Draw("same") %ENDCODE%

Note how the variable q needs to be multiplied by -2 in order to obtain the test-statistic -2*ln(L(A)/L(N)).

-- Mingshui Chen - 2016-11-24

Comments

<--/commentPlugin-->