Notes on Madgraph simulation of Higgs diphoton signal and background for NIU PHYS 474/790

You can find these notes online at http://www.niu.edu/spmartin/madgraph/madhiggs.html . You can also review how to start Madgraph at http://www.niu.edu/spmartin/madgraph . Notes on the tutorial are here: http://www.niu.edu/spmartin/madgraph/madtutor.html . More notes on Madgraph input syntax are here: http://www.niu.edu/spmartin/madgraph/madsyntax.html .

The "heft" model

In order to have the necessary Higgs-gluon-gluon and Higgs-photon-photon effective couplings, we will need to import the "heft" (Higgs Effective Field Theory) model. First, it is useful to copy the default version of this model to a version of your own. Then you can modify your own version without affecting anybody else's. To do so: Your private version of the heft model now lives in the directory
/xdata/$USER/madgraph/models/heft_$USER/
This directory contains information about the particles in the model (which in the case of heft are the same as the Standard Model) and about their interactions. By default, the mass of the Higgs boson is set to 125 GeV. You can change it by editing line 127 of the file
/xdata/$USER/madgraph/models/heft_$USER/parameters.py
in the obvious way. However, we will just keep the default here. Below, we will discuss how to change the Higgs coupling to photons or gluons. Start Madgraph as usual, and import the model as discussed in the previous notes:

The Higgs diphoton signal

Now we can ask Madgraph to generate the process of Higgs production followed by decay to two photons:
    generate p p > h, h > a a
    output FIRST_HIGGS_RUN
    launch FIRST_HIGGS_RUN
Note that "a" is the symbol for a photon. Let's choose to run Pythia and DELPHES this time, so we respond with:
    1
    3
    0
This brings us to the prompt that allows us to edit the card files. For our first runs, let's just keep the default run_card.dat. In particular, this means that we will be generating 10000 events at 13 TeV center of mass energy. However, we'll want to modify the plot_card.dat file, so select:
    9
Specifically, we want to make a plot for the two-photon invariant mass, because that's where we'll see the Higgs bump. As we learned last time, we want to add a class to the plot_card.dat file, but this time for photons. So insert a line, after line 83, that reads:

a 0 #Class number 3

Here, "a" is the symbol, and "0" is the code, for photon. Now, delete all of the lines between "Begin PlotDefs" and "End PlotDefs". In their place, put just three lines that read:

pt 1 1
pt 3 2
mij 3 2

This tells Madgraph that we want to plot the following things:

  • The transverse momentum (pT) of the jet (class 1) with the highest pT.
  • The pT's of the photons (class 3) with the highest two pTs.
  • The invariant mass of the pair of photons (class 3) with the two highest pTs.

Before saving, also edit the lines between "Begin PlotRange" and "End PlotRange" to make the histogram bin size and the range for the plots appropriate. Change the "pt" and "mij" entries to read:

pt 5 0 200 # bin size, min value, max value

and

mij 2 0 200

Now you can save and exit the file, and type

    0
to start the run.

When you inspect the Delphes plots they should look something like the following. First, the number of jets in each event looks like this:

picture

Recall that we only had photons in the final state of the process we generated. However, Pythia added "extra" jets from the initial state in the parton shower and hadronization, as occurs in the real world.

The number of photons in each event looks like this:

picture

In most events, we see two photons. However, there is a significant probability to not detect one or both of the photons. There are also a few "extra" photons that are generated by Pythia.

Next, look at the transverse momentum (pT) distributions of the leading and sub-leading photons:

picture picture

Note that these peak at slightly less than half of the mass of the Higgs boson, because they share its energy. Most of the photons have pT > 35 GeV. This will later help us cut down on background. Finally, look at the invariant mass histogram for the signal events:

picture

This has a nice peak centered at MH = 125 GeV, as it should. There are a few events at 50 GeV or lower; these occur when there is an "extra" photon that gets erroneously paired with one of the photons from the Higgs. We can see that the mass resolution, according to Delphes, is a few GeV. This corresponds at least roughly to the real-world mass resolution from ATLAS and CMS.

The diphoton background

Now we have to worry about the background. To look at the background alone, do:

    generate p p > a a / h
    output HIGGS_BACKGROUND_RUN
    launch HIGGS_BACKGROUND_RUN
Note that the syntax here tells Madgraph that we want to generate all Feynman diagrams that result in two photons, EXCEPT the ones that involve the Higgs.

Now do the same modifications of the plot_card.dat, and start the run, exactly as above. This time the plots should look something like the following.

The pT distributions of the leading and sub-leading photons are:

picture picture

Note that quite unlike the signal, these are very sharply peaked at small pT. At the generator level, there was a cut in the run_card.dat file which prevented the photon pT from being too small, otherwise it would keep blowing up at very small pT. (The reason for the generator level cut is to prevent this.)

Similarly, the invariant mass resolution of the background is peaked at small mass:

picture

The Higgs diphoton signal and background together

Next, let us combine the Higgs diphoton signal and the diphoton background in a single run. To do so, we type:
    generate p p > h, h > a a
    add process p p > a a / h
    output HIGGS_SB_RUN1
    launch HIGGS_SB_RUN1
Let's choose to run Pythia, and both PGS and Delphes:
    1
    2
    3
    0
For this run, we will want lots of events, so that the signal has a chance to stick out over the background. So, edit the run_card.dat:
    2
and on line 32 change the number of events from 10000 to 100000. Also, we will impose a stronger generator-level cut on the transverse momentum of the photons. So, change line 107 to read

35 = pta ! minimum pt for the photons

We also don't want Madevent to waste time producing events with small diphoton invariant mass, because those won't be anywhere near the Higgs peak anyway. So change line 168 to read:

90 = mmaa ! min invariant mass of gamma gamma pair

Next, edit the plot_card.dat:

    9
Here, add a new class 3 for photons, as before, and eliminate all other plots, except:

mij 3 2

Also, change the PlotRange on line 144 to read:

mij 2 100 200

so that the histogram bins will be 2 GeV wide, and run from 100 GeV to 200 GeV. (Remember that we already cut events with invariant mass less than 90 GeV, so it is not meaningful to look at plots with invariant mass less than about 100 GeV, taking into account the photon energy resolution.)

Here are the plots that result, for 13 TeV LHC with Delphes:

picture

and for 8 TeV LHC with Delphes:

picture

Sadly, it takes some imagination to claim that one sees the Higgs peak near 125 GeV in these plots.

If you look at the Results and Events database, you will see that Madgraph calculated the cross-section for the generated events as 5.44 pb = 5440 fb for the 13 TeV LHC, and 3.64 pb = 3640 fb for the 8 TeV LHC. Therefore, the integrated luminosity corresponding to the events that were generated is:

(100000 events)/(5440 fb) = 18.4 fb-1 for the 13 TeV LHC

(100000 events)/(3640 fb) = 27.5 fb-1 for the 8 TeV LHC

The 13 TeV LHC hasn't happened yet. The 8 TeV LHC actually collected about 20 fb-1, so our simulated event sample corresponds, naively, to more events than were actually collected. However, Madgraph's estimate of the production cross-section is not extremely accurate. Also, the ATLAS and CMS collaborations did a more advanced treatment of both the signal and background, chose their cuts more carefully, and did a sophisticated statistical analysis, combining also with the H > ZZ > 4 leptons signal.

Modifying the Higgs effective couplings

To see the Higgs diphoton peak more clearly, we can "cheat" by increasing the Higgs coupling to two photons, or by increasing the Higgs coupling to two gluons. The names for these couplings are, respectively, GC_1 and GC_13, and they are defined in the file:
/xdata/$USER/madgraph/models/heft_$USER/couplings.py
Increasing one or both of these Higgs couplings can be used to simulate what would occur if one used a more accurate treatment of the production cross-section. It will give us an enhanced peak.

So, let's edit this file and just double the Higgs-photon-photon coupling compared to the Madgraph heft value. To do this, change line 13 of the file from

value = '-(AH*complex(0,1))',

to

value = '-2.0*(AH*complex(0,1))',

This increases the signal amplitude by a factor of 2, and therefore increases the signal cross-section by a factor of 4.

Now, redoing everything as before, being sure to quit Madgraph and restart it and import the model heft_$USER again. This should give plots that look like the ones below.

13 TeV LHC with Delphes:

picture

8 TeV LHC with Delphes:

picture

This time the peaks are visible to the naked eye, thanks to our "cheat" of multiplying the signal cross-section by 4.

Doing your own detector cut-based analysis

In this class, we've relied on Madgraph's own plots. That's what you should use for the homework. However, for more sophisticated analyses at the research frontier, you probably want to impose your own cuts at the detector-level. There are several ways to do this:
  • Take the Pythia output file, and run it through your own detector simulator.
  • Take the PGS or Delphes output root files, and analyze them yourself, putting on your own cuts.
  • Take the PGS or Delphes LHCO files, and analyze them yourself, putting on your own cuts. If you want to learn to do this, a good place to start is: http://madgraph.phys.ucl.ac.be/Manual/lhco.html
The ROOT and LHCO files are available on the Results and Events Database by html, or you can deduce where to download the files directly by sftp or scp.