TWiki> Main Web>TWikiUsers>JeanetteLorenz>HistFitterTutorialOutsideAtlas>HistFitterAdvancedTutorialOutsideAtlas (2015-03-27, JeanetteLorenz) EditAttachPDF

Introduction

We discuss some advanced HistFitter questions and topics on this twiki page. This tutorial continues the introductionary HistFitter tutorial at HistFitterTutorialOutsideAtlas.

Advanced options in workspace building and fitting

Using Minos

We frequently check errors on fit parameters before and after fittting a model (e.g. in the so-called pull distributions). By default the error calculation is done by HESSE, which gives symmetric errors and is based on the second derivative of the loglikelihood at the minimum.

Alternatively, the error calculation can be run using MINOS, which gives asymmetric errors and is calculating errors using a 'hill climbing' algorithm that finds the points where the likelihood increases by +0.5 units in 1-D case. More information about fit parameter error calculation can be found in Minuit documentation.

HELP Note: Histfitter is using Minuit2 for the minimization.

The HESSE error calculation is only reliable if the likelihood is parabolical, but there are many cases when likelihoods are not everywhere parabolical, e.g. have kinks and double minimums. One needs to check NLL/PLL plots [example below] to see if the parabolic approximation is valid otherwise one needs to switch from the default HESSE error calculation to MINOS. This section illustrates on how to do this.

To activate the Minos calculation for only a subset of parameters use the option 

-m alpha1,alpha2,etc

If you want to run Minos for all parameters use: 

-m ALL

We illustrate the MINOS error calculation with the configuation file from HistFitterTutorialOutsideAtlas#Improving_the_limits_configuring

First, we create the workspace and perform a fit with the default error calculation performed by HESSE. We also draw likelihood (NLL) plots with -D "likelihood" option in order to see the effect: 

HistFitter.py -t -w -f -F excl -D "likelihood" -i analysis/tutorial/MyShapeFitExample.py

Now, we can draw profile likelihood curves (PLL) for the fit parameters using the Minos error calculation. To do so, on the command line, use the minos option -m alpha1,alpha2,etc or -m ALL in combination with the draw option -D likelihood. 

HistFitter.py -t -w -f -F excl -D "likelihood" -m "ALL" -i analysis/tutorial/MyShapeFitExample.py

After each step you can check out NLL/PLL plot stored in result folder, results/MyShapeFitExample/can_NLL__RooExpandedFitResult_afterFit_alpha* (And you will also see the plots displayed during running the commands above.)

Moreover we can compare the results in terms of the error calculation from first example using HESSE and from the second example using Minos. For fast comparison just look in the log file, and find the print out of the fit parameter initial/final values and errors.

HESSE (Example 1 result): 

  RooFitResult: minimized FCN value: -25.4406, estimated distance to minimum: 1.05749e-05
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Constant Parameter    Value     
  --------------------  ------------
  binWidth_obs_x_SR_metmeff2Jet_0    1.0000e+01
  binWidth_obs_x_SR_metmeff2Jet_1    1.0000e+01
  binWidth_obs_x_SR_metmeff2Jet_2    1.0000e+01
         nom_alpha_JES    0.0000e+00
  nom_alpha_KtScaleTop    0.0000e+00
   nom_alpha_KtScaleWZ    0.0000e+00
           nominalLumi    1.0000e+00

    Floating Parameter  InitialValue    FinalValue +/-  Error     GblCorr.
  --------------------  ------------  --------------------------  --------
                  Lumi    1.0000e+00    1.0000e+00 +/-  3.88e-02  0.008928
             alpha_JES    0.0000e+00   -9.1858e-02 +/-  9.18e-01  0.750468
      alpha_KtScaleTop    0.0000e+00    3.9326e-02 +/-  9.32e-01  0.383819
       alpha_KtScaleWZ    0.0000e+00    4.7297e-02 +/-  8.52e-01  0.469419
                mu_SIG    1.0000e+00    1.9798e-02 +/-  2.49e-01  0.778151

MINOS (Example 2 result): 

  RooFitResult: minimized FCN value: -25.4406, estimated distance to minimum: 1.05749e-05
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 MINOS=0 

    Constant Parameter    Value     
  --------------------  ------------
  binWidth_obs_x_SR_metmeff2Jet_0    1.0000e+01
  binWidth_obs_x_SR_metmeff2Jet_1    1.0000e+01
  binWidth_obs_x_SR_metmeff2Jet_2    1.0000e+01
         nom_alpha_JES    0.0000e+00
  nom_alpha_KtScaleTop    0.0000e+00
   nom_alpha_KtScaleWZ    0.0000e+00
           nominalLumi    1.0000e+00

    Floating Parameter  InitialValue    FinalValue (+HiError,-LoError)    GblCorr.
  --------------------  ------------  ----------------------------------  --------
                  Lumi    1.0000e+00    1.0000e+00 (+3.90e-02,-3.90e-02)  0.008928
             alpha_JES    0.0000e+00   -9.1858e-02 (+7.34e-01,-9.32e-01)  0.750468
      alpha_KtScaleTop    0.0000e+00    3.9326e-02 (+9.65e-01,-9.25e-01)  0.383819
       alpha_KtScaleWZ    0.0000e+00    4.7297e-02 (+9.30e-01,-6.97e-01)  0.469419
                mu_SIG    1.0000e+00    1.9798e-02 (+1.49e-01,-1.98e-02)  0.778151

Note the assymmetric errors in the second example. Also note that the error on the parameters alpha_JES and alpha_KtScaleWZ is now less than 1. This behaviour (which is sometimes called "profiling" ATLAS internally) means that the fit was able to reduce the uncertainties with respect to what you originally provided the fit as input. This behaviour can be desired and expected or also be too aggressive. Make thus sure you understand this behaviour whenever you see it!

Different types for systematic uncertainties

Possible types for systematic uncertainties in HistFitter

We discussed already in the first part of this tutorial HistFitterTutorialOutsideAtlas#Systematic_uncertainties that HistFitter allows to implement systematic uncertainties as different types, depending on the effect the systematic uncertainty should have. In particular we introduced the types

  • overallSys if a systematic uncertainty only varies the scale of the nominal histogram
  • histoSys if the systematic uncertainty varies both the shape and the scale of a nominal histogram
  • normHistoSys if the systematic uncertainty only varies the shape of a nominal histogram

HistFitter provides much more possibilities to implement systematic uncertainties. The most common types (more types, in particular also user defined systematic uncertainties, exist) are:

Basic systematic methods in HistFactory
overallSys uncertainty of the global normalization, not affecting the shape
histoSys correlated uncertainty of shape and normalization
shapeSys uncertainty of statistical nature applied to a sum of samples, bin by bin
Additional systematic methods in HistFitter
overallNormSys overallSys constrained to conserve total event count in a list of region(s)
normHistoSys histoSys constrained to conserve total event count in a list of region(s)
normHistoSysOneSide one-sided normHistoSys uncertainty built from tree-based or weight-based inputs
normHistoSysOneSideSym symmetrized normHistoSysOneSide
overallHistoSys factorized normalization shape and uncertainty; described with overallSys and histoSys respectively
overallNormHistoSys overallHistoSys in which the shape uncertainty is modeled with a normHistoSys and the global normalization uncertainty is modeled with an overallSys
shapeStat shapeSys applied to an individual sample

(This table is extracted from the HistFitter paper: http://arxiv.org/abs/1410.1280 )

Which systematics type should be used?

Given these rather large amount of systematics types, one could get confused which type should be used in a certain situation. Alex Koutsman made a nice flowchart to illustrate the decision process.

HistFitter_systematics.png

"Normalized" systematic uncertainties and the transfer factor approach

A basic ingredient in many analyses is the extrapolation from control to signal or validation regions for estimating a major background in a semi-data-driven way. The idea is to fix a background (model) in a fit to data in control regions that were specifically designed for the specific backgrounds. Using then Monte Carlo simulations, this background model is extrapolated to signal regions in order to obtain a background estimate in them. This extrapolation from control to signal regions is only possible if the kinematic variables used in the extrapolation are well described by Monte Carlo simulation. This assumption can be checked in validation regions kinematically close to the signal regions, but usually containing a lower signal contamination.

This analysis strategy is illustrated in the two plots taken from the HistFitter paper (http://arxiv.org/abs/1410.1280) - look into this paper if you are not familiar with this concept of extrapolation!

cartoon_CRVRSR.png

A schematic view of an analysis strategy with multiple control, validation and signal regions. All regions can have single- or multiple bins, as illustrated by the dashed lines. The extrapolation from the control to the signal regions is verified in the validation regions that lie in the extrapolation phase space.

overview-1.png

A typical analysis strategy flow with HistFitter.

The concept of transfer factors

The extrapolation from control to signal (validation) regions itself is described by so-called transfer factors, where we need one transfer factor for every sample (i.e. background process) and every extrapolation. The transfer factor from a control to a signal (validation) region is defined as ratio of the Monte Carlo yields in the signal (validation) region to the yield in the control region.

Mathematical description

Show Hide

For the mathematical description we cite here the HistFitter paper:

" A key underlying concept of the fitting procedure is the implicit use of ratios of expected event counts, called transfer factors, or TFs, for each simulated background process between each SR and CR.

The normalized background predictions used in fit are:

\begin{eqnarray} N_p({\rm CR}, {\rm est.}) &=& \mu_{p}\times {\rm MC}_p({\rm CR},{\rm init.}) \,, \nonumber \\ N_p({\rm SR}, {\rm est.}) &=& \mu_{p}\times {\rm MC}_p({\rm SR},{\rm init.}) \,, \nonumber (1) \end{eqnarray}

where $N_p({\rm CR}, {\rm est.})$ ($N_p({\rm SR}, {\rm est.})$) is the CR (SR) background estimate for each simulated physics processes $p$ considered in the analysis, ${\rm MC}_p({\rm CR}, {\rm init.})$ (${\rm MC}_p({\rm SR}, {\rm init.})$) is the initial estimate as obtained from the MC simulation, and $\mu_{p}$ is the normalization factor as obtained in the fit to data.

In the fit, the background estimate(s) is (are) typically driven by the statistics in the CR(s). Define $N_p({\rm CR}, {\rm fit})$ as the fitted number of events for process $p$ in the CR. Then equivalently:

\begin{eqnarray} N_p({\rm SR}, {\rm est.}) &=& \mu_{p} \times {\rm MC}_p({\rm SR},{\rm init.}) \nonumber \\ &\equiv& N_p({\rm CR}, {\rm fit}) \times \left [\frac{{\rm MC}_p({\rm SR}, {\rm init.})}{{\rm MC}_p({\rm CR}, {\rm init.})}\right ] \,. \nonumber (2) \end{eqnarray}

The ratio appearing in the square brackets of Eq. (2) is defined as the transfer factor TF. The two notation are equivalent in terms of the fit: there is a normalization quantity, derived from the fit to data, multiplied by a fixed constant. In other words, albeit the fit uses Eq. (1) internally, $N_p({\rm SR}, {\rm est.})$ can always be interpreted as a TF multiplied by the fitted number of background events in the CR.

An important feature of the TF approach is that systematic uncertainties on the predicted background processes can be (partially) canceled in the extrapolation; a virtue of the ratio of MC estimates. The total uncertainty on the number of background events in the SR is then a combination of the statistical uncertainties in the CR(s) and the residual systematic uncertainties of the extrapolation. For this reason, CRs are sometimes defined by somewhat loose selection criteria, in order to increase CR data event statistics without significantly increasing residual uncertainties in the TFs, which in turn reduces the extrapolation uncertainties to the SR. "

The usage of "normalized" systematic uncertainties

If using the transfer factor approach systematic uncertainties on the predicated background processes might be partially cancel, if having a similar impact in control and signal (validation) regions. This well-known virtue of Monte Carlo ratios can help in reducing the total uncertainty on background estimates in the signal regions, under the assumption that the control regions were carefully chosen.

Exactly this transfer factor approach is meant when using systematic uncertainties with Norm in the name.

To use these types, we first need to explain to HistFitter which regions should be taken as normalization regions. The normalization are used to normalize the systematic uncertainties to. Let's take e.g. an up variation of a nominal histogram. We could e.g. define the nominal histogram as normalization region. When doing so, we would normalize our up histogram such that it has the same integral as our nominal histogram. Effectively this means, that the up variation will only differ in shape to the nominal histogram, but no longer has a different normalization (== scale).

Usually, we take as normalization regions all the control regions that the analysis has, but the choice really depends on your needs - you should thus discuss this with your reviewers.

Getting a bit more concrete - let's have again a look at the example configuration file analysis/tutorial/MyConfigExample.py which already reflects a quite complex analysis.

Open this file: 

nedit analysis/tutorial/MyConfigExample.py &

There are three systematic uncertainties implemented as histoSys: JES, KtScaleWZ, KtScaleTop

First we want to see the impact of these systematic uncertainties in the signal region SS. In order to see this run a background-only fit with the validation regions turned on (doValidation=True, line 50 in the configuration file): 

HistFitter.py -t -w -f analysis/tutorial/MyConfigExample.py

Create from this a systematic table: 

SysTable.py -w results/MyConfigExample/BkgOnly_combined_NormalMeasurement_model_afterFit.root -c SS_metmeff2Jet -o systable_SS.tex

In the configuration file analysis/tutorial/MyConfigExample.py you find the following lines: 

topSample.setNormRegions([("SLWR","nJet"),("SLTR","nJet")])
wzSample.setNormRegions([("SLWR","nJet"),("SLTR","nJet")])
bgSample.setNormRegions([("SLWR","nJet"),("SLTR","nJet")])

This tells HistFitter to use the two control regions SLWR and SLTR as normalization regions, to be used for the backgrounds Top, WZ and BG.

Now switch in the configuration file the three uncertainties JES, KtScaleWZ, KtScaleTop to the type overallNormHistoSys. (You will see that these lines already exist in the configuration file as commented lines, uncomment those and comment the original ones.)

Running again: 

HistFitter.py -t -w -f analysis/tutorial/MyConfigExample.py
SysTable.py -w results/MyConfigExample/BkgOnly_combined_NormalMeasurement_model_afterFit.root -c SS_metmeff2Jet -o systable_SS_2.tex

Now compare the two tables: 

kompare systable_SS.tex  systable_SS_2.tex

Indeed you find that the three uncertainties JES, KtScaleWZ, KtScaleTop are reduced if using the systematics type overallNormHistoSys with normalization regions.

Behind the scenes: internal construction of "normalized" systematic uncertainties

Show Hide

The construction of the normalized histograms contains multiple steps and this thus not necessarily very intuitive. For this reason, we illustrate specifically for the systematics type overallNormHistoSys on how the original input histograms are modified.

Various normalization steps are carried out as follows:

  • Temporary histograms are built containing the sum of the yields in all normalization regions (this is done for the nominal, the up and the down histogram).
  • First normalization step: Up and down temporary histograms are scaled to the yield in the nominal temporary histogram.
  • The obtained scaling factor is applied to the original up and down histograms.
  • If the systematics type contains an overall in the name - second normalization step: Up and down histograms are scaled to the yields in the nominal histogram. The resulting histograms are internally used with an histoSys type, the scaling factor with a overallSys type.

Warning, important Caution: Care must be taken when applying Norm type systematics, as only samples that carry a normalization factor mu_... in the fit can be used. Otherwise the uncertainty on the total event count or normalization is not taken into account for this sample, which can lead to underestimation of the sample uncertainty in the signal region.

Exercise:

In the data file data/MyConfigExample.root we have stored both the original, unnormalized histograms and the normalized histograms. Open this file: 

root data/MyConfigExample.root

And list the histograms with .ls. Histograms like hWZNom_SLTR_obs_nJet or hWZJESHigh_SLTR_obs_nJet are the unnormalized, original histograms (syntax: h=+ sample + Nom or systematics + region + =_obs_ + variable). Histograms like hTopJESHigh_SLWRSLTRNorm or hTopNom_SLWRSLTRNorm (so with a region string identifying the normalization region SLWRSLTR and an additional Norm) are the temporary histograms used in the first bullet point above. Finally, histograms like hWZJESHigh_SLTR_obs_nJetNorm are the final, normalized histograms that are used further (they differ from the name of the original histogram only by a Norm at the end).

Compare some of these histograms to see the effect of using overallNormHistoSys.

Additional tools for cross-checks

Visualizing systematic uncertainties versus the nominal histogram in a specific region

Visualizing systematic variations as input by the user is often a very helpful procedure to immediately see when systematic uncertainties look unexpected with respect to a nominal histogram. Weird systematic variations may lead to unstable fits or too a large, unexpected reduction of the uncertainties originally put .

For these reasons, we introduce in this example two very simple (and in fact not yet very sophisticated) macros helping the user to visualize the impact of a systematic uncertainty with respect to a nominal histogram.

Warning, important As these macros are very simple, they are tuned very much to this specific example. You can easily rewrite them to your use case. On the longer time-scale we hope to integrate this functionality into the main HistFitter scripts. Checks as we present in this exercise are usually quite helpful in debugging some model if a fit is failing.

Outgoing from the example HistFitterAdvancedTutorialOutsideAtlas#The_usage_of_normalized_systemat above, we illustrate the effect that using overallNormHistoSys has on the size of systematic uncertainties. The two macros scripts/plotSyst.C and scripts/plotUpDown.C are already updated to the example HistFitterAdvancedTutorialOutsideAtlas#The_usage_of_normalized_systemat above, such that they work really 'out-of-the-box'. (Don't hesitate to open both files to understand what the macros are doing and - eventually - to modify them according your needs.)

Just run now: 

root scripts/plotSyst.C

This will create two plots, as you can also see in the calls in scripts/plotSyst.C. The one plot, c_hTopJESHigh_SS_obs_metmeff2Jet_JES_metmeff2Jet.eps, is showing the nominal Top background histogram in the signal region SS together with the up and down variations caused by the JES systematic uncertainty as originally input by the user or obtained from -t. The second plot, c_hTopJESHigh_SS_obs_metmeff2JetNorm_JES_metmeff2Jet.eps, is showing again the nominal histogram with the two up and down JES variations, but this time the JES uncertainty was normalized as described in the example above using the control regions SLWR and SLTR as normalization regions.

  • not using a normalized systematic type, but just histoSys:

c_hTopJESHigh_SS_obs_metmeff2Jet_JES_metmeff2Jet.png

  • using overallNormHistoSys

c_hTopJESHigh_SS_obs_metmeff2JetNorm_JES_metmeff2Jet.png

Comparing these two plots you clearly see the impact of reducing systematic uncertainties when using the transfer factor approach, so when using a normalized systematic uncertainties type.

How to debug failed fits

Rebuilding hypothesis tests outside of -p

We have seen in HistFitterTutorialOutsideAtlas#Details_on_the_output_of_p that various fits are performed when running hypothesis tests with -p, among those:

  • a free fit of the background + signal model
  • a fit with the signal strength fixed to 1
  • a fit with the signal strength fixed to 0

If one or multiple of these fits fail, no or no reliable CLs value can be calculated. Typical problems that lead to failures are:

  • the background + signal model is very complicated and there are (fit) instabilities related to systematic uncertainties,
  • the signal contribution is very large in comparison to a very small background prediction and we have very little observed data around.

When the calculation of the CLs value is failing it is mandatory for the user to check first the log file to understand which of the fits done in -p is really failing. Afterwards, the user can try to rebuild this fit, using the functionalities HistFitter provides for -f. For technical reasons -f can provide a little bit more verbose output than -p, which is helpful for debugging any problematic fit. In addition, the user can make before- and after-fit plots and correlations matrices which can also give some hints.

We don't have really complicated background+signal models in the tutorial and we also tried to make sure that none of our examples is failing. (Although this depends a little bit on the actual Root version used.) Therefore, in the following we will try to reproduce the fits done in the example 'Details on the output pf -p' to illustrate on how to rebuild fits outside of -p. Of course, you can also take one of your private problematic models and try to debug this along the lines described in this exercise. (Your tutor is happy to help you with any questions regarding problems appearing there.)

If you have not done so yesterday, please create the log file of the example HistFitterTutorialOutsideAtlas#Details_on_the_output_of_p: 

HistFitter.py -p analysis/tutorial/MyUserAnalysis.py 2>&1 | tee out.log

We will use this log file to check the results we get in the fits rebuilt.

The option -C in HistFitter

This option allows to fix certain parameters to a values provided by the user in a fit done with -f. The syntax is 

-C parameter1:value1,parameter2:value2,... etc.

Also see the examples below.

Reproducing a fit of -p with the signal strength fixed to 0.

For this just run: 

HistFitter.py -f -C "mu_Sig:0." analysis/tutorial/MyUserAnalysis.py

The fit result of this is: 

  RooFitResult: minimized FCN value: -0.273073, estimated distance to minimum: 0.00027612
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Constant Parameter    Value     
  --------------------  ------------
  binWidth_obs_x_UserRegion_cuts_0    1.0000e+00
  binWidth_obs_x_UserRegion_cuts_1    1.0000e+00
                mu_Sig    0.0000e+00
         nom_alpha_cor    0.0000e+00
         nom_alpha_ucb    0.0000e+00
         nom_alpha_ucs    0.0000e+00
  nom_gamma_stat_UserRegion_cuts_bin_0    1.0000e+00
           nominalLumi    1.0000e+00

    Floating Parameter  InitialValue    FinalValue +/-  Error     GblCorr.
  --------------------  ------------  --------------------------  --------
                  Lumi    1.0000e+00    1.0000e+00 +/-  3.88e-02  0.000000
             alpha_cor    0.0000e+00    1.5405e-01 +/-  9.69e-01  0.122812
             alpha_ucb    0.0000e+00    2.4486e-01 +/-  9.00e-01  0.185254
             alpha_ucs    0.0000e+00    0.0000e+00 +/-  9.93e-01  0.000000
  gamma_stat_UserRegion_cuts_bin_0    1.0000e+00    1.0614e+00 +/-  1.98e-01  0.190423

Comparing to out.log, you see that this fit result is pretty similar to the fit done in b) and c), although not identical. The reason for the small differences are small differences in the fit options that are used in -f and in -p. Most importantly, the luminosity was fixed and did not change in -f. If you are interested in the other more subtle differences, don't hesitate to ask your tutor. For the common practice, the two fit results are close enough.

Reproducing a fit of -p with the signal strength set to 1.

Run: 

HistFitter.py -f -C "mu_Sig:1." analysis/tutorial/MyUserAnalysis.py

You should get the following result: 

  RooFitResult: minimized FCN value: -0.35802, estimated distance to minimum: 9.17589e-06
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Constant Parameter    Value     
  --------------------  ------------
  binWidth_obs_x_UserRegion_cuts_0    1.0000e+00
  binWidth_obs_x_UserRegion_cuts_1    1.0000e+00
                mu_Sig    1.0000e+00
         nom_alpha_cor    0.0000e+00
         nom_alpha_ucb    0.0000e+00
         nom_alpha_ucs    0.0000e+00
  nom_gamma_stat_UserRegion_cuts_bin_0    1.0000e+00
           nominalLumi    1.0000e+00

    Floating Parameter  InitialValue    FinalValue +/-  Error     GblCorr.
  --------------------  ------------  --------------------------  --------
                  Lumi    1.0000e+00    9.9878e-01 +/-  3.88e-02  0.038996
             alpha_cor    0.0000e+00   -2.1561e-01 +/-  9.67e-01  0.230151
             alpha_ucb    0.0000e+00   -1.7437e-01 +/-  9.85e-01  0.196699
             alpha_ucs    0.0000e+00   -8.2217e-02 +/-  9.93e-01  0.099913
  gamma_stat_UserRegion_cuts_bin_0    1.0000e+00    9.0940e-01 +/-  1.94e-01  0.282557

Comparing to out.log you see that this is pretty similar to the result obtained in e).

Help The fit in f) uses an Asimov set and also a signal strength of 1. Also this fit can be rebuild, essentially using the HistFitter option -a, but it is a little bit more tricky. Ask your tutor in case you're interested.

Help The other fits in -p are free fits of the background+signal model and can be easily rebuild by using the standard -f option without using the option -C.

Help For all these examples you can also create plots and correlation matrices etc. using the options -d and -D.

Modifying the scan range of the upper limit scan

Although the scan range in an upper limit scan is automatically determined when running -l there are cases where the user wants to modify the scan range. In this example we illustrate on how to do this easily. We will work on a typical example where the user does not get reasonable results back without modifying the scan range.

Rewriting the config file such that we can give a signal point via the command line.

In this example we work again with the config file analysis/tutorial/MySimpleChannelConfig.py which we used in the tutorial part 1 to produce a contour plot.

Let's open this file with your favorite text editor, e.g. 

nedit analysis/tutorial/MySimpleChannelConfig.py &

Scrolling down to line 82 you see all the different signal models that were tested in -p yesterday: 

sigSamples = ["SU_680_310_0_10","SU_440_145_0_10","SU_200_160_0_10","SU_440_340_0_10","SU_440_100_0_10","SU_120_130_0_10","SU_600_280_0_10","SU_320_115_0_10","SU_360_175_0_10","SU_920_310_0_10","SU_280_205_0_10","SU_1080_340_0_10","SU_40_280_0_10","SU_760_160_0_10","SU_200_115_0_10","SU_280_280_0_10","SU_40_160_0_10","SU_520_280_0_10","SU_120_220_0_10","SU_680_220_0_10","SU_40_115_0_10","SU_920_190_0_10","SU_320_130_0_10","SU_440_280_0_10","SU_360_100_0_10","SU_120_160_0_10","SU_1080_190_0_10","SU_840_250_0_10","SU_120_100_0_10","SU_120_340_0_10","SU_840_280_0_10","SU_80_115_0_10","SU_840_130_0_10","SU_320_175_0_10","SU_120_205_0_10","SU_520_100_0_10","SU_400_130_0_10","SU_360_310_0_10","SU_160_115_0_10","SU_1000_310_0_10","SU_40_220_0_10","SU_440_130_0_10","SU_1000_190_0_10","SU_80_220_0_10","SU_840_160_0_10","SU_120_145_0_10","SU_440_175_0_10","SU_360_280_0_10","SU_320_145_0_10","SU_400_160_0_10","SU_1000_340_0_10","SU_600_310_0_10","SU_320_190_0_10","SU_840_310_0_10","SU_200_220_0_10","SU_440_205_0_10",
"SU_360_205_0_10","SU_120_280_0_10","SU_1080_130_0_10","SU_160_145_0_10","SU_520_250_0_10","SU_840_100_0_10","SU_160_220_0_10","SU_120_190_0_10","SU_40_205_0_10","SU_280_250_0_10","SU_80_145_0_10","SU_200_175_0_10","SU_840_190_0_10","SU_240_145_0_10","SU_160_205_0_10","SU_400_115_0_10","SU_440_250_0_10","SU_600_340_0_10","SU_80_100_0_10","SU_520_190_0_10","SU_1160_190_0_10","SU_80_130_0_10","SU_400_190_0_10","SU_400_175_0_10","SU_600_130_0_10","SU_1080_100_0_10","SU_200_340_0_10","SU_1160_310_0_10","SU_440_160_0_10","SU_240_220_0_10","SU_200_100_0_10","SU_240_130_0_10","SU_360_130_0_10","SU_1000_250_0_10","SU_920_130_0_10","SU_240_190_0_10","SU_520_340_0_10","SU_40_175_0_10","SU_240_100_0_10","SU_400_145_0_10","SU_40_145_0_10","SU_240_205_0_10","SU_1080_280_0_10","SU_600_250_0_10","SU_360_145_0_10","SU_520_130_0_10","SU_1000_130_0_10","SU_440_310_0_10","SU_600_160_0_10","SU_920_280_0_10","SU_760_280_0_10","SU_280_190_0_10","SU_280_175_0_10","SU_120_310_0_10","SU_440_220_0_10","SU_1000_220_0_10","SU_1160_250_0_10",
"SU_400_205_0_10","SU_160_160_0_10","SU_1000_280_0_10","SU_1000_160_0_10","SU_400_100_0_10","SU_760_190_0_10","SU_680_160_0_10","SU_840_220_0_10","SU_360_340_0_10","SU_1080_220_0_10","SU_360_250_0_10","SU_760_130_0_10","SU_440_115_0_10","SU_240_160_0_10","SU_200_310_0_10","SU_200_145_0_10","SU_600_220_0_10","SU_280_130_0_10","SU_520_220_0_10","SU_1080_160_0_10","SU_40_190_0_10","SU_1160_160_0_10","SU_280_310_0_10","SU_920_160_0_10","SU_80_190_0_10","SU_40_310_0_10","SU_1160_130_0_10","SU_40_250_0_10","SU_40_100_0_10","SU_400_220_0_10","SU_40_340_0_10","SU_1000_100_0_10","SU_120_175_0_10","SU_280_220_0_10","SU_760_340_0_10","SU_240_115_0_10","SU_440_190_0_10","SU_1160_340_0_10","SU_600_100_0_10","SU_200_250_0_10","SU_280_145_0_10","SU_200_190_0_10","SU_200_205_0_10","SU_760_250_0_10","SU_120_250_0_10","SU_80_175_0_10","SU_40_130_0_10","SU_920_250_0_10","SU_80_160_0_10","SU_240_175_0_10","SU_280_100_0_10","SU_1080_310_0_10","SU_920_340_0_10","SU_120_115_0_10","SU_1160_100_0_10","SU_280_340_0_10",
"SU_1160_220_0_10","SU_200_130_0_10","SU_160_175_0_10","SU_360_220_0_10"]

Hard-coding all the signal models may work if you work with only O(100) signal models, but it gets quite cumbersome and the config file quite to difficult to read if you use O(1000) models. We can do better.

Below line 82 you find a commented part: 

#if not 'sigSamples' in dir():
#    sigSamples=["SU_680_310_0_10"]

(which is much shorter).

Uncomment these two lines and comment line 82.

When you are now running -l like: 

HistFitter.py -l analysis/tutorial/MySimpleChannelConfig.py

you will only evaluate the upper limit for the signal model SU_680_310_0_10 (and not for all the other points that you run over yesterday with '-p'). (If you have not done the exercise 'Part 5: Making a complete exclusion contour plot (*)' before, please go back to this exercise first before running the command.)

Run now over a different signal model using the -g option and safe the output into some log file: 

HistFitter.py -l -g SU_360_220_0_10 analysis/tutorial/MySimpleChannelConfig.py 2>&1 | tee out.log

(You could also run over multiple signal models by giving a comma-separated list, e.g. SU_160_175_0_10,SU_360_220_0_10.)

You will find some worrisome output: 

<INFO> HypoTestTool: The computed upper limit is: 14.4325 +/- 0
<INFO> HypoTestTool:  expected limit (median) 0
<INFO> HypoTestTool:  expected limit (-1 sig) 0
<INFO> HypoTestTool:  expected limit (+1 sig) 0
<INFO> HypoTestTool:  expected limit (-2 sig) 0
<INFO> HypoTestTool:  expected limit (+2 sig) 0
TROOT::Append:0: RuntimeWarning: Replacing existing TH1: CLs_observed (Potential memory leak).
Error in <TGraphPainter::PaintGraph>: illegal number of points (0)
<INFO> HypoTestTool:  writing result plot to results/upperlimit_cls_poi_SU_360_220_0_10_Asym_CLs_grid_ts3.root.eps
Error in <TGraphPainter::PaintGraph>: illegal number of points (0)
Info in <TCanvas::Print>: eps file results/upperlimit_cls_poi_SU_360_220_0_10_Asym_CLs_grid_ts3.root.eps has been created
<ERROR> ConfigMgrCPP: All fits seem to have failed - cannot compute upper limit!
<INFO> ConfigMgrCPP: Now storing HypoTestInverterResult <debug_SU_360_220_0_10>
<INFO> ConfigMgrCPP: Done. Stored upper limit in file <results/MySimpleChannelAnalysis_upperlimit.root>
<INFO> HistFitter: Leaving HistFitter... Bye!

This tells you that the upper limit scan failed and you have not gotten any result. In the next step we will study the log file to understand what happened.

Reading the log file

It case that a upper limit scan is not successful it is helpful to study the output log file first before debugging the problem. In this exercise we explain the different parts of the log file and how the upper scan is set up in HistFitter internally.

Let's have a look at the log file: 

less out.log

You see the following parts: (also search for the line indicated by using /):

a) At the beginning you have a lot of messages from the ConfigMgr, like: 

   <INFO> ConfigManager: analysisName: MySimpleChannelAnalysis

This tells you that HistFitter is reading in the configuration that you provided.

b) Starting from the following line: 

   <INFO> HypoTestTool: >>> Running HypoTestInverter on the workspace combined

The part of '-l' itself starts. You see there a lot of messages indicating the configuration of '-l' and that HistFitter copies the background and the background+signal model from an object called 'ModelConfig'. (The ModelConfig object holds the PDF with associated information as e.g. the parameter of interest in a RooFit workspace. For more information on this consult a RooFit/RooStats tutorial.)

c) The next you see is that a free fit of the background + signal model is performed: 

   Info in <StandardHypoTestInvDemo>:  Doing a first fit to the observed data

... and further below ... 

   <INFO> HypoTestTool: StandardHypoTestInvDemo - Best Fit value : mu_SIG = 0.00105946 +/- 0.721623

The tool has thus determined the best fitted value of the signal strength. Based on this we determine the scan range for a first upper limit scan as: 

   <INFO> HypoTestTool: Doing a fixed scan  in interval : 0 , 14.4325

The lower value 0 just corresponds to a signal strength of 0, and the upper value is calculated as (fitted value of the signal strength + 20 * the error on the fitted signal strength)

d) You see in the following a number of hypothesis tests performed, starting with a signal strength of 0, and going up to the upper value of the scan range in equidistant steps. The number of steps was given by you in the config file, line 26: 

   configMgr.nPoints=20       # number of values scanned of signal-strength for upper-limit determination of signal strength.

At the end of all these hypothesis tests you see (again in the log file): 

   <INFO> HypoTestTool: Time to perform limit scan 
   <INFO> HypoTestTool: >>> Done running HypoTestInverter on the workspace combined
   [#0] WARNING:Eval -- HypoTestInverterResult::FindInterpolatedLimit - not enough points to get the inverted interval - return 1
   [#0] ERROR:Eval -- HypoTestInverterResult::GetLimitDistribution not  enough points -  return 0

This tells you that it was not possible to evaluate the upper limit in this first scan. In this particular case the reason is that none of the tested signal strengths resulted in a CLs value above 0.05. (We need to have at least one CLs value above and one below 0.05 in the upper limit scan such that we can determine the signal strength matching the crossing point of the upper limit curve with the line at 0.05.)

To understand this a bit better, go back to the beginning of the log file and search for 'CLs ='. You will find the results of all hypothesis tests done.

e) After the first upper limit scan, a second upper limit scan is started. You find it when looking for the second occurrence of 

   <INFO> HypoTestTool: >>> Running HypoTestInverter on the workspace combined

The purpose of this second upper limit scan is to re-evaluate the upper limit obtained in first scan more precisely by restricting the scan range to [0, 1.2 * (expected upper limit of first scan + 2 sigma)]

In this case, the second scan is failing completely, as the upper limit in the first scan could not be evaluated.

f) The ERROR messages at the end of the log file want to tell you that there was a problem: 

   <ERROR> ConfigMgrCPP: All fits seem to have failed - cannot compute upper limit!

You should not use the result of this scan for any physics interpretations. Therefore, the result is stored with a special name: debug_SU_360_220_0_10 to give you an hint that you should use the output only to debug the result.

Recovering the failed upper limit scan.

Obviously, the scan range that was automatically determined in a free fit of the signal+background model was much too large. HistFitter provides the possibility to fix the scan range of the second upper limit scan by hand.

Open again the configuration file: 

   nedit analysis/tutorial/MySimpleChannelConfig.py &

and uncomment line 29: 

   ## set scan range for the upper limit
   configMgr.scanRange = (0., 1.)

Now rerun 

   HistFitter.py -l -g SU_360_220_0_10 analysis/tutorial/MySimpleChannelConfig.py 2>&1 | tee out2.log

This results in some meaningful result: 

   <INFO> HypoTestTool: The computed upper limit is: 0.363536 +/- 0
   <INFO> HypoTestTool:  expected limit (median) 0.373768
   <INFO> HypoTestTool:  expected limit (-1 sig) 0.257714
   <INFO> HypoTestTool:  expected limit (+1 sig) 0.562129
   <INFO> HypoTestTool:  expected limit (-2 sig) 0.190629
   <INFO> HypoTestTool:  expected limit (+2 sig) 0.832839

HELP Note: if the scan range set by the user is too small, i.e. the smaller than the expected upper limit + 2 sigma, the user will see a WARNING at the end of the output.

Further information on the discovery fit

The general concepts of discovery fit are discussed on the main histfitter tutorial page HistFitterTutorialOutsideAtlas#Model_independent_upper_limits_t. Here we will run a discovery hypothesis test for the tutorial example configuration file, analysis/tutorial/MyOneBinExample.py. This is a simple counting experiment, with input from TTrees (one bin, one region), and it was already used to run "exclusion fits" in Part 1 of the main tutorial HistFitterTutorialOutsideAtlas. If you open the configuration file - analysis/tutorial/MyOneBinExample.py, you can see the settings we use for discovery and exclusion fits, e.g. 

configMgr.nTOYs=5000  # number of toys when doing frequentist calculator
configMgr.calculatorType=0 # 2=asymptotic calculator, 0=frequentist calculator

There are two blocks of model configurations, first for the "Discovery fit" (also called model independent signal limit fit) and second for the "Exclusion fit" (also called model independent signal limit fit), you can specify the fit type using -F option, -F 'disc' for discovery fit and -F 'excl' for exclusion fit. For the discovery fit we have only one signal region, and as its model independent input signal is a "dummy" signal model without signal uncertainties and with a signal expectation of 1 event. You can put more than one control region but only one signal region. This is essential in order to stay model independent (otherwise you are alrey including an assumption on how the specific signal model behaves in different signal regions). For the same reason signal contamination in control regions are ignored in the discovery fit setup.

The first step is to create a workspace, 

HistFitter.py -t -w -f -F disc analysis/tutorial/MyOneBinExample.py

The option -F disc sets the quantity myFitType to FitType.Discovery, which activates the following setup for the signal sample, 

srBin.addDiscoverySamples(["Discovery"],[1.],[0.],[10000.],[kMagenta])

1) 'dummy' signal count is set to 1, 2) no systematic is applied on signal

Second step is to run a discovery hypothesis test: 

HistFitter.py -F disc -z analysis/tutorial/MyOneBinExample.py

The result should come out as: 

Results HypoTestCalculator_result: 
 - Null p-value = 0.6132 +/- 0.00688746
 - Significance = -0.287669 +/- 0.0179936 sigma
 - Number of Alt toys: 2500
 - Number of Null toys: 5000
 - Test statistic evaluated on data: -5.33661e-05
 - CL_b: 0.6132 +/- 0.00688746
 - CL_s+b: 0.6796 +/- 0.00933261
 - CL_s: 1.10828 +/- 0.019662

In the output one can see that there are 5000 toys generated for the null hypothesis, and 2500 toys for the alternate hypothesis.

Further important options or features in HistFitter

In this section we discuss further important options that are available in HistFitter and simplify life. wink

Build histograms from user input

This functionality is very useful if you just want to quickly define or build some histogram, e.g. to test what happens if the observed data were at that specific value and not at that other value...

You see an example for this again in analysis/tutorial/MyShapeFitExample.py: 

dataSample = Sample("Data",kBlack)
dataSample.setData()
dataSample.buildHisto([1.,1.,5.,15.,4.,0.],"SR","metmeff2Jet",0.1,0.1)
#dataSample.buildStatErrors([1.,1.,2.4,3.9,2.,0.],"SR","metmeff2Jet")

The first command 

dataSample.buildHisto([1.,1.,5.,15.,4.,0.],"SR","metmeff2Jet",0.1,0.1)

builds a data histogram with 6 bins with the bin content [1.,1.,5.,15.,4.,0.] in the region SR for the distribution metmeff2Jet. The last two numbers indicate the lower bin edge of the lowest bin and the bin width.

The second command, which is commented, 

dataSample.buildStatErrors([1.,1.,2.4,3.9,2.,0.],"SR","metmeff2Jet")

would build statistical errors on the histogram created by buildHisto, but you not need this for data histograms. It's just put for illustration here.

When you run 

HistFitter.py -t -F excl analysis/tutorial/MyShapeFitExample.py

and open the file 

root data/MyShapeFitExample.root

You can have a look how the histogram built looks like 

hData_SR_obs_metmeff2Jet->Draw()

This usual option has been used in this tutorial quite often, as it allows us to avoid using real ATLAS data.

Blinding of signal, control and validation regions

When you do not want to use real data an alternative to the example above of building some fake data is to blind the data. Blinding the data means that the data will be assumed to be equal to the total background estimate before the fit.

This possibility is usually used when you design an analysis and are only interested in looking at the expected sensitivity, but do not want to bias yourself with the knowledge about the real data.

Looking again at analysis/tutorial/MyShapeFitExample.py you find the lines: 

configMgr.blindSR = False # Blind the SRs (default is False)
configMgr.blindCR = False # Blind the CRs (default is False)
configMgr.blindVR = False # Blind the VRs (default is False)

In this case our configuration file has only one signal region, but no control or validation regions. It thus does only make sense to blind the signal region. Try this by modifying the following line: 

configMgr.blindSR = True # Blind the SRs (default is False)

Running now 

HistFitter.py -t -w -F excl -f -m "all" analysis/tutorial/MyShapeFitExample.py

and producing a yields table: 

YieldsTable.py -s Top,WZ,SM_GG_onestepCC_425_385_345 -c SR_metmeff2Jet -w results/MyShapeFitExample/Exclusion_SM_GG_onestepCC_425_385_345_combined_NormalMeasurement_model_afterFit.root -o MyYieldsTable.tex

Looking at MyYieldsTable.tex you find indeed that the observed data corresponds to the yields of the Top + WZ backgrounds.

You will notice that the blinding did only use the background yields, but not the signal yield. If you want to add the signal yield as well to the 'fake' observed data, activate the option: 

configMgr.useSignalInBlindedData = True

Exercise: convince yourself that now indeed the 'fake' observed data corresponds to Top + WZ + signal.

Recycling already built histograms

This very useful option allows you to reuse histograms that you have already built some time ago and to only build the histograms you really don't have from trees using -t.

To illustrate how it works please run again 

HistFitter.py -t -F excl analysis/tutorial/MyShapeFitExample.py

This will - as you know - produce the .root file data/MyShapeFitExample.root. Rename this file, e.g. as data/MyShapeFitExample_template.root

Now uncomment the following lines in the configuration file analysis/tutorial/MyShapeFitExample.py: 

#activate using of background histogram cache file to speed up processes
configMgr.useCacheToTreeFallback = True # enable the fallback to trees
configMgr.useHistBackupCacheFile = True # enable the use of an alternate data file
configMgr.histBackupCacheFile =  "data/MyShapeFitExample_template.root" # the data file of your previous fit (= the backup cache file)

To illustrate the benefits of these lines, let's now run over another signal model, e.g. SM_GG_onestepCC_505_305_105: 

HistFitter.py -w -F excl -g SM_GG_onestepCC_505_305_105 analysis/tutorial/MyShapeFitExample.py

Warning, important Note: it's very important to run only -w and not -t! If running -t you will reproduce all histograms, like this only the histograms you don't have in data/MyShapeFitExample_template.root.

Lines like this: 

<INFO> ConfigManager:   Sample: SM_GG_onestepCC_505_305_105
<INFO> PrepareHistos: Could not get histogram <hSM_GG_onestepCC_505_305_105Nom_SR_obs_metmeff2Jet> from backupCacheFile data/MyShapeFitExample_template.root, trying cacheFile
<INFO> PrepareHistos: Could not find histogram <hSM_GG_onestepCC_505_305_105Nom_SR_obs_metmeff2Jet> in data/MyShapeFitExample.root, trying from tree

inform you that for the signal point SM_GG_onestepCC_505_305_105 no histogram could be found in the file data/MyShapeFitExample_template.root and neither in the file data/MyShapeFitExample.root. The required histogram will thus be produced from tree.

This option is very helpful if creating tons of different workspaces with different signal models in a very fast and efficient way.

Setting sample specific weights

In some situations it is necessary to give a specific background or signal sample a specific weight - e.g. when you want to test just a tenth of the signal model that you have implemented or if you have derived some reweighting procedure for some background and want to apply it to just that background etc.

For such cases we have the option 

sigSample.addSampleSpecificWeight("0.001") #set whatever weight you want, but it needs to be given as string

Looking at our configuration file analysis/tutorial/MyShapeFitExample.py you can find this line and uncomment it to study its effect on the result. (You certainly know by now how to do this smile )

Final note

There are also further options available in HistFitter and in any case the package is growing. So don't hesistate to check the code itself and to send us your questions.

-- JeanetteLorenz - 2015-03-27

Edit | Attach | Watch | Print version | History: r1 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r1 - 2015-03-27 - JeanetteLorenz
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright &© 2008-2024 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
or Ideas, requests, problems regarding TWiki? use Discourse or Send feedback