{
"cells": [
{
"cell_type": "markdown",
"id": "f40b7065-5763-4abc-91ff-40437cf86517",
"metadata": {},
"source": [
"# Simple $e^+e^- \\to \\mu^+ \\mu^- $ Monte Carlo Event Generator example"
]
},
{
"cell_type": "markdown",
"id": "215df39b-69fc-4f74-b882-3c96f5526783",
"metadata": {},
"source": [
"Differential cross section:\n",
"$$\n",
"\\frac{d \\sigma}{ d \\Omega} = \\frac{\\alpha^2}{4 s } ( 1+ \\cos^2(\\theta))\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "f7124c5f-9f86-4a0a-ac08-6e30cc55e01d",
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"import random\n",
"\n",
"import hist\n",
"import numpy as np\n",
"\n",
"import pylhe"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "f2553a00-1d72-4d4b-ba6d-c552ca0a7b59",
"metadata": {},
"outputs": [],
"source": [
"alpha = 1 / 127.4 # alpha QED\n",
"aqcd = 0.1075 # alpha QCD\n",
"EB = 209\n",
"s = (2 * EB) ** 2\n",
"theta_min = 0\n",
"theta_max = math.pi\n",
"phi_min = 0\n",
"phi_max = 2 * math.pi"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "7242869e-7fbb-4eee-9d49-8a9d63724e35",
"metadata": {},
"outputs": [],
"source": [
"# https://equation-database.readthedocs.io/en/latest/_autosummary/equation_database.isbn_9780471887416.html#equation_database.isbn_9780471887416.equation_6_32\n",
"\n",
"\n",
"def dsigma(s, theta, _phi):\n",
" return (math.cos(theta) ** 2 + 1) / 4 * alpha**2 / s"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "147dfe17-a5d9-4b68-a07f-6f8d8aa390a1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Estimated integral: 1.4779947196668872e-09\n",
"Real integral: 1.477056777521761e-09\n"
]
}
],
"source": [
"# Monte Carlo integration\n",
"\n",
"\n",
"def monte_carlo_integration(\n",
" func, s, theta_min, theta_max, phi_min, phi_max, num_samples\n",
"):\n",
" theta_samples = [random.uniform(theta_min, theta_max) for _ in range(num_samples)]\n",
" phi_samples = [random.uniform(phi_min, phi_max) for _ in range(num_samples)]\n",
" func_values = [\n",
" (phi_max - phi_min)\n",
" * (theta_max - theta_min)\n",
" * func(s, theta, phi)\n",
" * math.sin(theta)\n",
" for theta, phi in zip(theta_samples, phi_samples)\n",
" ]\n",
" maximum = np.max(func_values)\n",
" integral = np.mean(func_values)\n",
" return integral, maximum\n",
"\n",
"\n",
"# Parameters\n",
"\n",
"num_samples = 1_000_000\n",
"\n",
"# Perform integration\n",
"result, maximum = monte_carlo_integration(\n",
" dsigma, s, theta_min, theta_max, phi_min, phi_max, num_samples\n",
")\n",
"print(f\"Estimated integral: {result}\")\n",
"\n",
"# https://equation-database.readthedocs.io/en/latest/_autosummary/equation_database.isbn_9780471887416.html#equation_database.isbn_9780471887416.equation_6_33\n",
"sigma = 4 * math.pi / 3 * alpha**2 / s\n",
"print(f\"Real integral: {sigma}\")"
]
},
{
"cell_type": "markdown",
"id": "692b7fad-a519-49b8-9dbc-d7d8c27eaba9",
"metadata": {},
"source": [
"## Weighted Events\n",
"\n",
"The information about the distribution is carried by the weights of the events."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "5c26e5e3-8e1e-470f-96a6-376a6e294bf9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" 5 0 1.8882916063e-15 1.7472400000e+05 7.8492935636e-03 1.0750000000e-01\n",
" 11 -1 0 0 0 0 0.00000000e+00 0.00000000e+00 2.09000000e+02 2.09000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
" -11 -1 0 0 0 0 0.00000000e+00 0.00000000e+00 -2.09000000e+02 2.09000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
" 22 2 1 2 0 0 0.00000000e+00 0.00000000e+00 0.00000000e+00 4.18000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
" 13 1 3 3 0 0 -7.52255283e+01 -1.61873499e+02 -1.08715639e+02 2.09000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
" -13 1 3 3 0 0 7.52255283e+01 1.61873499e+02 1.08715639e+02 2.09000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
"\n"
]
}
],
"source": [
"num_samples = 1_000_000\n",
"integ = 0.0\n",
"lheevents = []\n",
"for _i in range(num_samples):\n",
" theta = random.uniform(theta_min, theta_max)\n",
" phi = random.uniform(phi_min, phi_max)\n",
" sig = (\n",
" (phi_max - phi_min)\n",
" * (theta_max - theta_min)\n",
" * dsigma(s, theta, phi)\n",
" * math.sin(theta)\n",
" / num_samples\n",
" )\n",
" integ += sig\n",
" # Fill the LHE event\n",
" e = pylhe.LHEEvent(\n",
" eventinfo=pylhe.LHEEventInfo(\n",
" nparticles=5,\n",
" pid=0,\n",
" weight=sig, # The individual weight per event\n",
" scale=s,\n",
" aqed=alpha,\n",
" aqcd=aqcd,\n",
" ),\n",
" particles=[\n",
" pylhe.LHEParticle(\n",
" id=11,\n",
" status=-1,\n",
" mother1=0,\n",
" mother2=0,\n",
" color1=0,\n",
" color2=0,\n",
" px=0.0,\n",
" py=0.0,\n",
" pz=EB,\n",
" e=EB,\n",
" m=0.0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" pylhe.LHEParticle(\n",
" id=-11,\n",
" status=-1,\n",
" mother1=0,\n",
" mother2=0,\n",
" color1=0,\n",
" color2=0,\n",
" px=0.0,\n",
" py=0.0,\n",
" pz=-EB,\n",
" e=EB,\n",
" m=0.0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" pylhe.LHEParticle(\n",
" id=22,\n",
" status=2,\n",
" mother1=1,\n",
" mother2=2,\n",
" color1=0,\n",
" color2=0,\n",
" px=0.0,\n",
" py=0.0,\n",
" pz=EB - EB,\n",
" e=EB + EB,\n",
" m=0.0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" pylhe.LHEParticle(\n",
" id=13,\n",
" status=1,\n",
" mother1=3,\n",
" mother2=3,\n",
" color1=0,\n",
" color2=0,\n",
" px=EB * math.sin(theta) * math.cos(phi),\n",
" py=EB * math.sin(theta) * math.sin(phi),\n",
" pz=EB * math.cos(theta),\n",
" e=EB,\n",
" m=0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" pylhe.LHEParticle(\n",
" id=-13,\n",
" status=1,\n",
" mother1=3,\n",
" mother2=3,\n",
" color1=0,\n",
" color2=0,\n",
" px=-EB * math.sin(theta) * math.cos(phi),\n",
" py=-EB * math.sin(theta) * math.sin(phi),\n",
" pz=-EB * math.cos(theta),\n",
" e=EB,\n",
" m=0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" ],\n",
" )\n",
" lheevents.append(e)\n",
"# Fill the LHE init\n",
"lheinit = pylhe.LHEInit(\n",
" initInfo=pylhe.LHEInitInfo(\n",
" beamA=11,\n",
" beamB=-11,\n",
" energyA=EB,\n",
" energyB=EB,\n",
" PDFgroupA=0,\n",
" PDFgroupB=0,\n",
" PDFsetA=0,\n",
" PDFsetB=0,\n",
" weightingStrategy=3,\n",
" numProcesses=1,\n",
" ),\n",
" procInfo=[\n",
" pylhe.LHEProcInfo(\n",
" xSection=integ,\n",
" error=0,\n",
" unitWeight=1,\n",
" procId=1,\n",
" )\n",
" ],\n",
" weightgroup={},\n",
" LHEVersion=3,\n",
")\n",
"lhef = pylhe.LHEFile(lheinit, lheevents)\n",
"print(lhef.events[0].tolhe(rwgt=False, weights=False))\n",
"lhef.tofile(\"weighted.lhe.gz\", rwgt=False, weights=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2fca49b0-2d93-476a-9b6b-835b72a24bb0",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "c94b9451-e7cc-4f40-9c1e-63dd8e180595",
"metadata": {},
"outputs": [],
"source": [
"arr = pylhe.to_awkward(\n",
" pylhe.LHEFile.fromfile(\"weighted.lhe.gz\", with_attributes=True).event\n",
")\n",
"hist1 = hist.Hist.new.Reg(20, 0, math.pi).Weight()\n",
"hist1.fill(arr.particles.vector[:, -1].theta, weight=arr.eventinfo.weight)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b43498cd-1fe7-49c1-a24c-5b66cf4e257f",
"metadata": {},
"outputs": [],
"source": [
"hist1.plot()"
]
},
{
"cell_type": "markdown",
"id": "77d19d7a-07c7-4da3-81c6-31e716f3b15d",
"metadata": {},
"source": [
"## Unweighted Events\n",
"\n",
"The information about the distribution is carried by the distribution/number of the events."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "a43da3e7-5263-4557-8024-cb08ae00f01b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" 5 0 1.4779947197e-15 1.7472400000e+05 7.8492935636e-03 1.0750000000e-01\n",
" 11 -1 0 0 0 0 0.00000000e+00 0.00000000e+00 2.09000000e+02 2.09000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
" -11 -1 0 0 0 0 0.00000000e+00 0.00000000e+00 -2.09000000e+02 2.09000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
" 22 2 1 2 0 0 0.00000000e+00 0.00000000e+00 0.00000000e+00 4.18000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
" 13 1 3 3 0 0 1.31998175e+02 -7.87144118e+01 -1.41638706e+02 2.09000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
" -13 1 3 3 0 0 -1.31998175e+02 7.87144118e+01 1.41638706e+02 2.09000000e+02 0.00000000e+00 0.0000e+00 9.0000e+00\n",
"\n"
]
}
],
"source": [
"num_samples = 1_000_000\n",
"integ = 0.0\n",
"lheevents = []\n",
"while len(lheevents) < num_samples:\n",
" theta = random.uniform(theta_min, theta_max)\n",
" phi = random.uniform(phi_min, phi_max)\n",
" sig = (\n",
" (phi_max - phi_min)\n",
" * (theta_max - theta_min)\n",
" * dsigma(s, theta, phi)\n",
" * math.sin(theta)\n",
" )\n",
" if sig / maximum > random.uniform(\n",
" 0, 1\n",
" ): # Pick events randomly according to their contribution\n",
" e = pylhe.LHEEvent(\n",
" eventinfo=pylhe.LHEEventInfo(\n",
" nparticles=5,\n",
" pid=0,\n",
" weight=result / num_samples, # The same weight per event\n",
" scale=s,\n",
" aqed=alpha,\n",
" aqcd=aqcd,\n",
" ),\n",
" particles=[\n",
" pylhe.LHEParticle(\n",
" id=11,\n",
" status=-1,\n",
" mother1=0,\n",
" mother2=0,\n",
" color1=0,\n",
" color2=0,\n",
" px=0.0,\n",
" py=0.0,\n",
" pz=EB,\n",
" e=EB,\n",
" m=0.0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" pylhe.LHEParticle(\n",
" id=-11,\n",
" status=-1,\n",
" mother1=0,\n",
" mother2=0,\n",
" color1=0,\n",
" color2=0,\n",
" px=0.0,\n",
" py=0.0,\n",
" pz=-EB,\n",
" e=EB,\n",
" m=0.0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" pylhe.LHEParticle(\n",
" id=22,\n",
" status=2,\n",
" mother1=1,\n",
" mother2=2,\n",
" color1=0,\n",
" color2=0,\n",
" px=0.0,\n",
" py=0.0,\n",
" pz=EB - EB,\n",
" e=EB + EB,\n",
" m=0.0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" pylhe.LHEParticle(\n",
" id=13,\n",
" status=1,\n",
" mother1=3,\n",
" mother2=3,\n",
" color1=0,\n",
" color2=0,\n",
" px=EB * math.sin(theta) * math.cos(phi),\n",
" py=EB * math.sin(theta) * math.sin(phi),\n",
" pz=EB * math.cos(theta),\n",
" e=EB,\n",
" m=0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" pylhe.LHEParticle(\n",
" id=-13,\n",
" status=1,\n",
" mother1=3,\n",
" mother2=3,\n",
" color1=0,\n",
" color2=0,\n",
" px=-EB * math.sin(theta) * math.cos(phi),\n",
" py=-EB * math.sin(theta) * math.sin(phi),\n",
" pz=-EB * math.cos(theta),\n",
" e=EB,\n",
" m=0,\n",
" lifetime=0,\n",
" spin=9.0,\n",
" ),\n",
" ],\n",
" )\n",
" lheevents.append(e)\n",
"lheinit = pylhe.LHEInit(\n",
" initInfo=pylhe.LHEInitInfo(\n",
" beamA=11,\n",
" beamB=-11,\n",
" energyA=EB,\n",
" energyB=EB,\n",
" PDFgroupA=0,\n",
" PDFgroupB=0,\n",
" PDFsetA=0,\n",
" PDFsetB=0,\n",
" weightingStrategy=3,\n",
" numProcesses=1,\n",
" ),\n",
" procInfo=[\n",
" pylhe.LHEProcInfo(\n",
" xSection=result,\n",
" error=0,\n",
" unitWeight=1,\n",
" procId=1,\n",
" )\n",
" ],\n",
" weightgroup={},\n",
" LHEVersion=3,\n",
")\n",
"lhef = pylhe.LHEFile(lheinit, lheevents)\n",
"print(lhef.events[0].tolhe(rwgt=False, weights=False))\n",
"lhef.tofile(\"unweighted.lhe\", rwgt=False, weights=False)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "e840fba3-e797-497c-9ad8-e48451928614",
"metadata": {},
"outputs": [
{
"ename": "AttributeError",
"evalue": "'LHEFile' object has no attribute 'event'",
"output_type": "error",
"traceback": [
"\u001b[31m---------------------------------------------------------------------------\u001b[39m",
"\u001b[31mAttributeError\u001b[39m Traceback (most recent call last)",
"\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[11]\u001b[39m\u001b[32m, line 1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m arr = pylhe.to_awkward(\u001b[43mpylhe\u001b[49m\u001b[43m.\u001b[49m\u001b[43mLHEFile\u001b[49m\u001b[43m.\u001b[49m\u001b[43mfromfile\u001b[49m\u001b[43m(\u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43munweighted.lhe\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mwith_attributes\u001b[49m\u001b[43m=\u001b[49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\u001b[43m.\u001b[49m\u001b[43mevent\u001b[49m)\n\u001b[32m 2\u001b[39m hist2 = hist.Hist.new.Reg(\u001b[32m20\u001b[39m, \u001b[32m0\u001b[39m, math.pi).Weight()\n\u001b[32m 3\u001b[39m hist2.fill(arr.particles.vector[:, -\u001b[32m1\u001b[39m].theta, weight=arr.eventinfo.weight)\n",
"\u001b[31mAttributeError\u001b[39m: 'LHEFile' object has no attribute 'event'"
]
}
],
"source": [
"arr = pylhe.to_awkward(\n",
" pylhe.LHEFile.fromfile(\"unweighted.lhe\", with_attributes=True).event\n",
")\n",
"hist2 = hist.Hist.new.Reg(20, 0, math.pi).Weight()\n",
"hist2.fill(arr.particles.vector[:, -1].theta, weight=arr.eventinfo.weight)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c9ab56a1-50bc-405a-b7e5-c7d6b61881de",
"metadata": {},
"outputs": [],
"source": [
"hist2.plot()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c9be8352-6ab8-45c2-a4ac-fa458ee41652",
"metadata": {},
"outputs": [],
"source": [
"hist1.plot_ratio(hist2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "932f8f9b-b7d9-4a6e-9f34-f2f90aeb86a3",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}