MC: Understand the generator status

Contents

MC: Understand the generator status¶

The interplay of generator status and simulator status is discussed in the chaptor on the simulator status.

import awkward as ak
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import uproot

from lcio_checks.mc.simulation import add_simulation_info
from lcio_checks.util import config, load_or_make

f = uproot.open(f"{config['data_dir']}/P2f_z_eehiq.root")["MyLCTuple"]
mc = f.arrays(filter_name="mc*", entry_stop=-1)
mc = add_simulation_info(mc)

def validate_my_explanation(mc):
    mc4 = mc[mc.mcgst == 4]
    assert all(ak.num(mc4.mcgst) == 2)
    assert all(mc4.mcpdg[:, 0] == 11)
    assert all(mc4.mcpdg[:, 1] == -11)
    assert all(ak.flatten(mc4.mcmox) == 0.8750143)
    assert all(ak.flatten(mc4.mcmoy) == 0)
    assert all(np.abs(ak.flatten(mc4.mcmoz)) == 125)
    assert ak.min(mc4.mcene) == ak.max(mc4.mcene)
    assert abs(mc4.mcene[0][0] - 125) < 0.1  # E_e^2 = (M_e^2 + p_z_e^2)
    assert all(ak.flatten(mc4.mcvtx) == 0)
    assert all(ak.flatten(mc4.mcvty) == 0)
    assert all(ak.flatten(mc4.mcvtx == mc4.mcepx))
    assert all(ak.flatten(mc4.mcvty == mc4.mcepy))
    assert all(ak.flatten(mc4.mcvtz == mc4.mcepz))


validate_my_explanation(mc)

from lcio_checks.mc.generation import (
    count_mc_per_generator_status,
    vertex_in_z_per_generator_status,
)

fig1 = count_mc_per_generator_status(mc)
fig2 = vertex_in_z_per_generator_status(mc)

../_images/MC_generator_status_3_0.png

../_images/MC_generator_status_3_1.png

created_in_simu = np.unique(ak.flatten(mc[mc.isCreatedInSimulation].mcgst))
assert len(created_in_simu) == 1
assert created_in_simu[0] == 0  # Generator status 0

Meaning of the generator status¶

Adapted from LCIO docs

0: Created in simulation
1: Undecayed particle, stable in the generator
2: Particle decayed in the generator
3: Documentation line (used for overlay?)
4: Beam parameters

Particle PDGs (per generator status)¶

@load_or_make(["pdg_counts.csv"])
def load_pdg_counts(mc):
    pdg_counts = {}
    for pdg in np.unique(ak.flatten(mc.mcpdg)):
        tmp_mc = mc[mc.mcpdg == pdg]
        pdg_counts[pdg] = {}
        for gen_status in range(5):
            pdg_counts[pdg][gen_status] = ak.sum(
                ak.num(tmp_mc[tmp_mc.mcgst == gen_status].mcpdg)
            )
    df = pd.DataFrame(pdg_counts)
    df.index.name = "generator status"
    return (df,)


(df,) = load_pdg_counts(mc, redo=True)
df = df.transpose()

	-20433	-20423	-20413	-20323	-20313	-20213	-10433	-10431	-10423	-10421	-10413	-10411	-10323	-10321	-10313	-10311	-10213	-10211	-4232	-4224	-4222	-4212	-4132	-4122	-4114	-4112	-3334	-3324	-3322	-3314	-3312	-3224	-3222	-3214	-3212	-3122	-3114	-3112	-2224	-2214	-2212	-2114	-2112	-1114	-533	-531	-523	-521	-435	-433	-431	-425	-423	-421	-415	-413	-411	-325	-323	-321	-315	-313	-311	-215	-213	-211	-16	-15	-14	-13	-12	-11	-5	-4	-3	-2	-1	1	2	3	4	5	11	12	13	14	16	21	22	23	91	92	94	111	113	115	130	211	213	215	221	223	225	310	311	313	315	321	323	325	331	333	335	411	413	415	421	423	425	431	433	435	443	511	521	1114	2112	2114	2212	2214	2224	3112	3114	3122	3212	3214	3222	3224	3312	3322	3324	4112	4114	4122	4132	4214	4224	4232	4324	10111	10113	10211	10213	10221	10223	10311	10313	10321	10323	10331	10333	10411	10413	10421	10423	10431	10433	20113	20213	20223	20313	20323	20333	20413	20423	20433	1000010020	1000010030	1000010040	1000010050	1000010060	1000010070	1000020030	1000020040	1000020050	1000020060	1000020070	1000020080	1000020090	1000030040	1000030050	1000030060	1000030070	1000030080	1000030090	1000030100	1000030110	1000040050	1000040060	1000040070	1000040080	1000040090	1000040100	1000040110	1000040120	1000050070	1000050080	1000050090	1000050100	1000050110	1000050120	1000050130	1000050140	1000060090	1000060100	1000060110	1000060120	1000060130	1000060140	1000070100	1000070110	1000070120	1000070130	1000070140	1000070150	1000070160	1000070170	1000080130	1000080140	1000080150	1000080160	1000080170	1000080180	1000080190	1000080200	1000090170	1000090180	1000090190	1000090200	1000090210	1000090220	1000090230	1000100180	1000100190	1000100200	1000100210	1000100220	1000100230	1000100240	1000100250	1000110210	1000110220	1000110230	1000110240	1000110250	1000110260	1000120210	1000120220	1000120230	1000120240	1000120250	1000120260	1000120270	1000120280	1000130220	1000130230	1000130240	1000130250	1000130260	1000130270	1000130280	1000130290	1000130300	1000140240	1000140250	1000140260	1000140270	1000140280	1000140290	1000140300	1000140310	1000140320	1000140330	1000140340	1000150260	1000150300	1000150310	1000150320	1000150330	1000150340	1000150350	1000150360	1000150370	1000160320	1000160330	1000160340	1000160350	1000160360	1000160370	1000160380	1000170350	1000170360	1000170370	1000170380	1000180340	1000180350	1000180360	1000180370	1000180380	1000180390	1000180400	1000180420	1000180430	1000190380	1000190390	1000190400	1000190410	1000190420	1000190430	1000200400	1000200410	1000200420	1000200430	1000200440	1000200450	1000200460	1000210430	1000210440	1000210450	1000210460	1000210470	1000210480	1000210500	1000220440	1000220450	1000220460	1000220470	1000220480	1000220490	1000220500	1000220510	1000220520	1000230470	1000230480	1000230490	1000230500	1000230510	1000230520	1000240480	1000240490	1000240500	1000240510	1000240520	1000240530	1000240540	1000240550	1000240560	1000240580	1000250510	1000250520	1000250521	1000250530	1000250540	1000250550	1000250570	1000250580	1000250590	1000260520	1000260530	1000260540	1000260550	1000260560	1000260570	1000260580	1000260590	1000260600	1000260610	1000270550	1000270560	1000270570	1000270580	1000270590	1000270600	1000270610	1000270620	1000280560	1000280570	1000280580	1000280590	1000280600	1000280610	1000280620	1000280640	1000290580	1000290590	1000290610	1000290620	1000290630	1000290650	1000320730	1000501100	1000501150	1000521140	1000521150	1000521180	1000551240	1000601330	1000601390	1000611350	1000621420	1000661480	1000661490	1000661530	1000681610	1000691580	1000691600	1000701590	1000701660	1000711650	1000721610	1000721660	1000731660	1000741660	1000741680	1000751730
generator status
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	1	0	2	0	0	14	0	0	0	0	168	0	99	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	364	0	0	0	0	0	16857	0	0	13244	15521	406	352578	0	0	0	0	0	0	0	0	0	0	380842	662	13173	15442	0	0	1220943	0	0	0	0	14696	34	0	302	15102	0	0	230	2	0	370	0	0	0	593	0	0	103	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	62535	0	51733	0	0	124	0	403	118	0	72	0	1	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	6154	2834	25	7	13	9	1324	18788	38	144	70	6	1	1	100	255	577	187	35	18	2	2	38	288	2108	1197	237	1	6	6	21	54	493	598	39	12	2	7	118	464	2011	334	115	2	5	17	105	443	258	2	4	3	32	105	514	59	57	11	3	1	33	60	37	17	7	3	3	13	264	180	213	34	15	9	18	108	236	53	31	3	5	10	42	320	155	158	6	1	1	1	11	27	58	116	9	7	2	3	6	11	22	113	30	42	8	4	5	2	2	3	15	9	26	4	7	3	3	18	17	25	19	13	1	1	28	5	10	4	3	3	13	16	17	7	22	1	1	2	8	5	2	2	2	4	7	19	12	11	2	2	3	3	10	4	2	2	1	12	6	26	16	23	10	6	2	2	13	4	23	7	6	1	12	9	47	30	35	5	6	1	1	1	6	9	1	41	13	13	3	3	1	3	5	53	32	56	11	15	7	4	1	6	6	32	7	18	3	8	2	10	13	21	12	17	10	8	1	1	1	4	1	7	10	2	2	1	1	1	1	1	1	1	1	1	1	2	1	1	2	1	1	5	1	1	1	1	2	1	3
1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2384	0	1586	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	7862	0	0	0	0	0	98993	2	0	63	85	83	1159686	0	0	0	0	0	0	0	0	0	0	1160100	79	64	84	2	0	325719	0	0	0	0	0	0	0	5625	98968	0	0	0	0	0	0	0	0	0	7899	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1615	0	2355	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
2	9	11	13	98	96	600	7	2	22	18	31	13	826	69	577	69	2015	846	2	1	2	1	3	36	2	2	0	1	31	1	21	18	162	26	139	470	23	105	139	130	0	140	0	95	1	1	1	1	6	50	75	37	221	486	28	247	210	179	2250	0	144	1696	5146	1011	12155	0	0	2	0	0	0	143940	2	811	2710	26064	16787	16567	26302	2692	811	2	143501	0	0	0	0	13893	122280	680	11904	28284	35369	83174	16661	910	0	0	12185	1082	6771	16961	947	5650	5610	1727	114	0	2107	191	779	322	27	187	269	22	522	241	37	70	51	9	2	1	1	98	0	114	0	130	126	98	20	488	128	19	138	15	18	32	2	2	5	29	1	1	2	2	1	793	1663	838	1944	686	1971	69	566	82	830	23	156	23	23	16	38	2	6	461	620	466	69	82	11	16	16	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	43199	2	811	3408	26367	16017	15579	26236	3120	811	2	43199	0	0	0	0	17297	93126	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
4	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	43199	0	0	0	0	0	0	0	0	0	0	43199	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0

Energy in the event per generator status¶

def _energy_per_generator_status(mc, per_event=True):
    if per_event:
        max_bin = ak.max(ak.sum(mc.mcene, axis=1))
    else:
        max_bin = ak.max(mc.mcene)
    bins = np.arange(0, max_bin, 2)

    fig, ax = plt.subplots(figsize=(7, 5))
    for i, gen_status in enumerate(range(5)):
        x = mc[mc.mcgst == gen_status].mcene
        if per_event:
            _x = ak.sum(x, axis=1)
        else:
            _x = ak.flatten(x)
        ax.hist(
            _x, bins=bins, label=f"{gen_status=}: {ak.sum(ak.num(x)):> 9}", alpha=0.5
        )
    if per_event:
        ax.set_xlabel("energy per event [GeV]")
    else:
        ax.set_xlabel("energy per MCP [GeV]")
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.axvline(250, color="black", label="250 GeV")
    # for i in range(2, 8):
    #     ax.axvline(i * 250, color="black")
    ax.legend(loc="upper left")
    return (fig,)


@load_or_make(["event_energy.png"])
def event_energy_per_generator_status(mc):
    return _energy_per_generator_status(mc, per_event=True)


event_energy_per_generator_status(mc, redo=True);

../_images/MC_generator_status_9_0.png

@load_or_make(["mcp_energy.png"])
def mcp_energy_per_generator_status(mc):
    return _energy_per_generator_status(mc, per_event=False)


mcp_energy_per_generator_status(mc, redo=True);

../_images/MC_generator_status_10_0.png

Appendix: Some more distributions per generator status¶

@load_or_make(["mcepz_log.png"])
def end_point_z_per_generator_status(mc):
    vtx_z = np.abs(mc.mcepz)
    bins1 = np.arange(0, ak.max(vtx_z), 100)
    fig, ax = plt.subplots(figsize=(7, 5))
    for i, gen_status in enumerate(range(5)):
        x = ak.flatten(vtx_z[mc.mcgst == gen_status])
        ax.hist(x, bins=bins1, label=f"{gen_status=}: {len(x):> 9}", alpha=0.5)
    ax.set_xlabel("z [mm]")
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.legend()
    return (fig,)


end_point_z_per_generator_status(mc);

../_images/MC_generator_status_12_0.png

def var_per_generator_status(mc, var="mcepz"):
    vals = getattr(mc, var)
    bins1 = np.linspace(-1, 1, 100)
    bins2 = np.arange(0, ak.max(np.abs(vals)), 100)
    fig, [ax1, ax2] = plt.subplots(ncols=2, figsize=(10, 5))
    for i, gen_status in enumerate(range(5)):
        kw = dict(color=f"C{i}", alpha=0.5)
        x = ak.flatten(vals[mc.mcgst == gen_status])
        ax1.hist(
            x, bins=bins1, label=f"{gen_status}: {np.sum(np.abs(x) < 1):> 9}", **kw
        )
        ax2.hist(np.abs(x), bins=bins2, label=f"{gen_status}: {len(x):> 9}", **kw)
    ax1.set_title(var)
    ax1.set_xlabel("energy [GeV]")
    # ax.set_xscale("log")
    # ax.set_yscale("log")
    ax1.legend(title="generatorStatus")
    ax2.legend(title="generatorStatus")
    ax2.set_xscale("log")
    ax2.set_yscale("log")
    return (fig,)


def per_generator_status(var):
    x = getattr(mc, var)
    if all(x[0][0] == ak.flatten(x)):
        print(f"{var} only takes single value: {x[0][0]}")
        return

    @load_or_make([f"per_generator_status_{var}.png"])
    def var_plot(var, mc):
        return var_per_generator_status(mc, var)

    return var_plot(var, mc)


# for var in mc.fields:
#     per_generator_status(var);

per_generator_status("mcmoz")
per_generator_status("mcepz")
per_generator_status("mcepy")
per_generator_status("mcmoy")
per_generator_status("mcvtx")
per_generator_status("mcvty")
per_generator_status("mcvtz");

../_images/MC_generator_status_13_0.png

../_images/MC_generator_status_13_1.png

../_images/MC_generator_status_13_2.png

../_images/MC_generator_status_13_3.png

../_images/MC_generator_status_13_4.png

../_images/MC_generator_status_13_5.png

../_images/MC_generator_status_13_6.png

previous

Reproduce the cuts on generator level

next

The generator status 4