<jupyter_start><jupyter_text>Appendix 1: Adding systematic uncertainties<NAME>, <NAME>, <NAME>, and <NAME> 2018-2019 In this tutorial we'll explain how to add systematic uncertainties to the MadMiner workflow. Note that the treatment of systematic uncertainties changed substantially with `MadMiner v0.6`, including changes to the MadMiner file specification. Please don't use files from older MadMiner versions with systematic uncertainties.  Preparations Before you execute this notebook, make sure you have running installations of MadGraph, Pythia, and Delphes.<jupyter_code>import os
import logging
import numpy as np

from madminer.core import MadMiner
from madminer.lhe import LHEReader
from madminer.sampling import SampleAugmenter
from madminer import sampling
from madminer.ml import ScoreEstimator
from madminer.fisherinformation import FisherInformation, profile_information, project_information
from madminer.plotting import plot_systematics, plot_fisher_information_contours_2d
from particle import Particle<jupyter_output><empty_output><jupyter_text>Please enter here the environment variable pointing to your MG5 installation folder.<jupyter_code>mg_dir = os.getenv("MG_FOLDER_PATH")<jupyter_output><empty_output><jupyter_text>MadMiner uses the Python `logging` module to provide additional information and debugging output. You can choose how much of this output you want to see by switching the level in the following lines to `logging.DEBUG` or `logging.WARNING`.<jupyter_code># MadMiner output
logging.basicConfig(
    format="%(asctime)-5.5s %(name)-20.20s %(levelname)-7.7s %(message)s",
    datefmt="%H:%M",
    level=logging.DEBUG,
)

# Output of all other modules (e.g. matplotlib)
for key in logging.Logger.manager.loggerDict:
    if "madminer" not in key:
        logging.getLogger(key).setLevel(logging.WARNING)<jupyter_output><empty_output><jupyter_text>1. Parameters and benchmarks We'll just load the MadMiner setup from the first part of this tutorial:<jupyter_code>miner = MadMiner()
miner.load("data/setup.h5")<jupyter_output>16:19 madminer.utils.inter DEBUG   HDF5 file does not contain is_reference field.
16:19 madminer.core        INFO    Found 2 parameters:
16:19 madminer.core        INFO       CWL2 (LHA: dim6 2, maximal power in squared ME: (2,), range: (-20.0, 20.0))
16:19 madminer.core        INFO       CPWL2 (LHA: dim6 5, maximal power in squared ME: (2,), range: (-20.0, 20.0))
16:19 madminer.core        INFO    Found 6 benchmarks:
16:19 madminer.core        INFO       sm: CWL2 = 0.00e+00, CPWL2 = 0.00e+00
16:19 madminer.core        INFO       w: CWL2 = 15.20, CPWL2 = 0.10
16:19 madminer.core        INFO       neg_w: CWL2 = -1.54e+01, CPWL2 = 0.20
16:19 madminer.core        INFO       ww: CWL2 = 0.30, CPWL2 = 15.10
16:19 madminer.core        INFO       neg_ww: CWL2 = 0.40, CPWL2 = -1.53e+01
16:19 madminer.core        INFO       morphing_basis_vector_5: CWL2 = -1.58e+01, CPWL2 = -1.94e+01
16:19 madminer.core        INFO    Found morphing setup with 6 components
16:19 madminer.core        INFO    Did n[...]<jupyter_text>2. Set up systematics, save settings This is where things become interesting: We want to model systematic uncertainties. The main function is `add_systematics()`, the keyword `effect` determines how the effect of the nuisance parameters on the event weights is calculated. For `effect="norm"`, the nuisance parameter rescales thee overall cross section of one or multiple samples. For `effect="pdf"`, its effect is calculated with PDF variations. Finally, with `effect="scale"` scale variations are used.Here we consider three nuisance parameters: one for the signal normalization, one for the background normalization, and one for scale uncertainties (which we here assume to be correlated between signal and background).<jupyter_code>miner.add_systematics(effect="norm", systematic_name="signal_norm", norm_variation=1.1)
miner.add_systematics(effect="norm", systematic_name="bkg_norm", norm_variation=1.2)
miner.add_systematics(effect="scale", systematic_name="scales", scale="mu")<jupyter_output><empty_output><jupyter_text>Again, we save our setup:<jupyter_code>miner.save("data/setup_systematics.h5")<jupyter_output>16:19 madminer.core        INFO    Saving setup (including morphing) to data/setup_systematics.h5<jupyter_text>3. Run MadGraph Now it's time to run MadGraph. MadMiner will instruct MadGraph to use its built-in `systematics` tool to calculate how the event weights change under the scale variation.<jupyter_code>miner.run(
    sample_benchmark="sm",
    mg_directory=mg_dir,
    mg_process_directory="./mg_processes/signal_systematics",
    proc_card_file="cards/proc_card_signal.dat",
    param_card_template_file="cards/param_card_template.dat",
    run_card_file="cards/run_card_signal_small.dat",
    log_directory="logs/signal",
    python2_override=True,
    systematics=["signal_norm", "scales"],
)
miner.run(
    sample_benchmark="sm",
    mg_directory=mg_dir,
    mg_process_directory="./mg_processes/bkg_systematics",
    proc_card_file="cards/proc_card_background.dat",
    param_card_template_file="cards/param_card_template.dat",
    run_card_file="cards/run_card_background.dat",
    log_directory="logs/background",
    python2_override=True,
    systematics=["bkg_norm", "scales"],
)<jupyter_output>16:23 madminer.utils.inter INFO    Generating MadGraph process folder from cards/proc_card_background.dat at ./mg_processes/bkg_systematics
16:23 madminer.core        INFO    Run 0
16:23 madminer.core        INFO      Sampling from benchmark: sm
16:23 madminer.core        INFO      Original run card:       cards/run_card_background.dat
16:23 madminer.core        INFO      Original Pythia8 card:   None
16:23 madminer.core        INFO      Original config card:    None
16:23 madminer.core        INFO      Copied run card:         /madminer/cards/run_card_0.dat
16:23 madminer.core        INFO      Copied Pythia8 card:     None
16:23 madminer.core        INFO      Copied config card:      None
16:23 madminer.core        INFO      Param card:              /madminer/cards/param_card_0.dat
16:23 madminer.core        INFO      Reweight card:           /madminer/cards/reweight_card_0.dat
16:23 madminer.core        INFO      Log file:                run_0.log
16:23 madminer.core        INFO    C[...]<jupyter_text>4. Load events from LHE file When adding LHE or Delphes files, use the `systematics` keyword to list which systematic uncertainties apply to which sample:<jupyter_code>lhe = LHEReader("data/setup_systematics.h5")

lhe.add_sample(
    lhe_filename="mg_processes/signal_systematics/Events/run_01/unweighted_events.lhe.gz",
    sampled_from_benchmark="sm",
    is_background=False,
    k_factor=1.1,
    systematics=["signal_norm", "scales"],
)

lhe.add_sample(
    lhe_filename="mg_processes/bkg_systematics/Events/run_01/unweighted_events.lhe.gz",
    sampled_from_benchmark="sm",
    is_background=True,
    k_factor=1.1,
    systematics=["bkg_norm", "scales"],
)<jupyter_output><empty_output><jupyter_text>The next steps are unaffected by systematics.<jupyter_code># Partons giving rise to jets
particles = [
    *Particle.findall(lambda p: p.pdgid.is_quark),
    *Particle.findall(pdg_name="g"),
]

lhe.set_smearing(
    pdgids=[int(p.pdgid) for p in particles],
    energy_resolution_abs=0.0,
    energy_resolution_rel=0.1,
    pt_resolution_abs=None,
    pt_resolution_rel=None,
    eta_resolution_abs=0.1,
    eta_resolution_rel=0.0,
    phi_resolution_abs=0.1,
    phi_resolution_rel=0.0,
)

lhe.add_observable(
    "pt_j1",
    "j[0].pt",
    required=False,
    default=0.0,
)
lhe.add_observable(
    "delta_phi_jj",
    "j[0].deltaphi(j[1]) * (-1.0 + 2.0 * float(j[0].eta > j[1].eta))",
    required=True,
)
lhe.add_observable(
    "met",
    "met.pt",
    required=True,
)

lhe.add_cut("(a[0] + a[1]).m > 122.0")
lhe.add_cut("(a[0] + a[1]).m < 128.0")
lhe.add_cut("pt_j1 > 30.0")
lhe.analyse_samples()
lhe.save("data/lhe_data_systematics.h5")<jupyter_output><empty_output><jupyter_text>A look at distributions The function `plot_systematics()` makes it easy to check the effect of the various nuisance parameters on a distribution:<jupyter_code>?plot_systematics
_ = plot_systematics(
    filename="data/lhe_data_systematics.h5",
    theta=np.array([0.0, 0.0]),
    observable="pt_j1",
    obs_label="$p_{T,j}$",
    obs_range=(30.0, 400.0),
    ratio_range=(0.7, 1.4),
)
_ = plot_systematics(
    filename="data/lhe_data_systematics.h5",
    theta=np.array([0.0, 0.0]),
    observable="delta_phi_jj",
    obs_label=r"$\Delta\phi$",
    obs_range=(-3.1, 3.1),
    ratio_range=(0.7, 1.4),
)<jupyter_output><empty_output><jupyter_text>5. Sampling Let's generate training data for the SALLY method.<jupyter_code>sampler = SampleAugmenter("data/lhe_data_systematics.h5", include_nuisance_parameters=True)

x, theta, t_xz, _ = sampler.sample_train_local(
    theta=sampling.benchmark("sm"),
    n_samples=100000,
    folder="./data/samples",
    filename="train_score_systematics",
)<jupyter_output><empty_output><jupyter_text>6. Training The SALLY estimator will learn the score in terms of both the physics parameters and the nuisance parameters. In our case, its output will therefore be a vector with five components.<jupyter_code>estimator = ScoreEstimator(n_hidden=(50,))

estimator.train(
    method="sally",
    x="data/samples/x_train_score_systematics.npy",
    t_xz="data/samples/t_xz_train_score_systematics.npy",
)

estimator.save("models/sally_systematics")<jupyter_output><empty_output><jupyter_text>7. Fisher information The Fisher information is now also a 5x5 matrix. Constraint terms on the nuisance parameters, representing our prior knowledge on their values, can be calculated with `FisherInformation.nuisance_constraint_information()`.<jupyter_code>fisher = FisherInformation("data/lhe_data_systematics.h5")
?fisher.full_information
info_sally, _ = fisher.full_information(
    theta=np.zeros(5),
    model_file="models/sally_systematics",
    luminosity=1000000.0,
    include_xsec_info=False,
)
info_nuisance = fisher.nuisance_constraint_information()
info = info_sally + info_nuisance
for row in info:
    print(" ".join(["{:6.2f}".format(entry) for entry in row]))<jupyter_output><empty_output><jupyter_text>Fisher contours and profiled information When looking at a subset of these five parameter, we can either ignore the other parameters ("project"), or conservatively pick the most pessimistic value of them ("profile" them). MadMiner provides the functions `profile_information()` and `project_information()` for this purpose. Let's do that for the parameter space of the two physics parameters:<jupyter_code>info_proj = project_information(info, [0, 1])
info_prof = profile_information(info, [0, 1])
_ = plot_fisher_information_contours_2d(
    [info_proj, info_prof],
    inline_labels=["Projected", "Profiled"],
    xrange=(-0.1, 0.1),
    yrange=(-0.1, 0.1),
)<jupyter_output><empty_output><jupyter_text>For illustration, we can look at just one physics parameter and one nuisance parameter:<jupyter_code>info_demo = project_information(info, [0, 4])

info_demo_profiled = np.zeros((2, 2))
info_demo_profiled[0, 0] = profile_information(info_demo, [0])
_ = plot_fisher_information_contours_2d(
    [info_demo, info_demo_profiled],
    inline_labels=[None, "Profiled"],
    xrange=(-0.2, 0.2),
    yrange=(-0.6, 0.6),
    xlabel="Parameter of interest",
    ylabel="Nuisance parameter",
)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10616080751827359, 'fraction_numerical': 0.03228332753219631, 'mean_word_length': 4.698066435299951, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 65, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '634599', 'license': 'MIT', 'n_tokens_mistral': 4108, 'n_tokens_neox': 3795, 'n_words': 936, 'path': 'examples/tutorial_particle_physics/A1_systematic_uncertainties.ipynb', 'repo_name': 'diana-hep/madminer', 'size': '11517'}