#+PROPERTY: header-args :exports both :output-dir results :session xs :kernel python3 * Init ** Required Modules #+NAME: e988e3f2-ad1f-49a3-ad60-bedba3863283 #+begin_src jupyter-python :exports both :tangle tangled/xs.py import numpy as np import matplotlib.pyplot as plt import monte_carlo #+end_src #+RESULTS: e988e3f2-ad1f-49a3-ad60-bedba3863283 ** Utilities #+NAME: 53548778-a4c1-461a-9b1f-0f401df12b08 #+BEGIN_SRC jupyter-python :exports both %run ../utility.py %load_ext autoreload %aimport monte_carlo %autoreload 1 #+END_SRC #+RESULTS: 53548778-a4c1-461a-9b1f-0f401df12b08 : The autoreload extension is already loaded. To reload it, use: : %reload_ext autoreload * Implementation #+NAME: 777a013b-6c20-44bd-b58b-6a7690c21c0e #+BEGIN_SRC jupyter-python :exports both :results raw drawer :exports code :tangle tangled/xs.py """ Implementation of the analytical cross section for q q_bar -> gamma gamma Author: Valentin Boettcher """ import numpy as np from scipy.constants import alpha # NOTE: a more elegant solution would be a decorator def energy_factor(charge, esp): """ Calculates the factor common to all other values in this module Arguments: esp -- center of momentum energy in GeV charge -- charge of the particle in units of the elementary charge """ return charge**4*(alpha/esp)**2/6 def diff_xs(θ, charge, esp): """ Calculates the differential cross section as a function of the azimuth angle θ in units of 1/GeV². Arguments: θ -- azimuth angle esp -- center of momentum energy in GeV charge -- charge of the particle in units of the elementary charge """ f = energy_factor(charge, esp) return f*((np.cos(θ)**2+1)/np.sin(θ)**2) def diff_xs_cosθ(cosθ, charge, esp): """ Calculates the differential cross section as a function of the cosine of the azimuth angle θ in units of 1/GeV². Arguments: cosθ -- cosine of the azimuth angle esp -- center of momentum energy in GeV charge -- charge of the particle in units of the elementary charge """ f = energy_factor(charge, esp) return f*((cosθ**2+1)/(1-cosθ**2)) def diff_xs_eta(η, charge, esp): """ Calculates the differential cross section as a function of the pseudo rapidity of the photons in units of 1/GeV^2. Arguments: η -- pseudo rapidity esp -- center of momentum energy in GeV charge -- charge of the particle in units of the elementary charge """ f = energy_factor(charge, esp) return f*(2*np.cosh(η)**2 - 1) def diff_xs_pt(pt, charge, esp): """ Calculates the differential cross section as a function of the transversal impulse of the photons in units of 1/GeV^2. Arguments: η -- transversal impulse esp -- center of momentum energy in GeV charge -- charge of the particle in units of the elementary charge """ f = energy_factor(charge, esp) return f*((esp/pt)**2/2 - 1) def total_xs_eta(η, charge, esp): """ Calculates the total cross section as a function of the pseudo rapidity of the photons in units of 1/GeV^2. If the rapditiy is specified as a tuple, it is interpreted as an interval. Otherwise the interval [-η, η] will be used. Arguments: η -- pseudo rapidity (tuple or number) esp -- center of momentum energy in GeV charge -- charge of the particle in units of the elementar charge """ f = energy_factor(charge, esp) if not isinstance(η, tuple): η = (-η, η) if len(η) != 2: raise ValueError('Invalid η cut.') def F(x): return np.tanh(x) - 2*x return 2*np.pi*f*(F(η[0]) - F(η[1])) #+END_SRC #+RESULTS: 777a013b-6c20-44bd-b58b-6a7690c21c0e * Calculations ** XS qq -> gamma gamma First, set up the input parameters. #+NAME: 7e62918a-2935-41ac-94e0-f0e7c3af8e0d #+BEGIN_SRC jupyter-python :exports both :results raw drawer η = 2.5 charge = 1/3 esp = 200 # GeV #+END_SRC #+RESULTS: 7e62918a-2935-41ac-94e0-f0e7c3af8e0d Set up the integration and plot intervals. #+begin_src jupyter-python :exports both :results raw drawer interval_η = [-η, η] interval = η_to_θ([-η, η]) interval_cosθ = np.cos(interval) interval_pt = η_to_pt([0, η], esp/2) plot_interval = [0.1, np.pi-.1] #+end_src #+RESULTS: *** Analytical Integratin And now calculate the cross section in picobarn. #+NAME: cf853fb6-d338-482e-bc55-bd9f8e796495 #+BEGIN_SRC jupyter-python :exports both :results raw file :file xs.tex xs_gev = total_xs_eta(η, charge, esp) xs_pb = gev_to_pb(xs_gev) tex_value(xs_pb, unit=r'\pico\barn', prefix=r'\sigma = ', prec=5, save=('results', 'xs.tex')) #+END_SRC #+RESULTS: cf853fb6-d338-482e-bc55-bd9f8e796495 : \(\sigma = \SI{0.05379}{\pico\barn}\) Compared to sherpa, it's pretty close. #+NAME: 81b5ed93-0312-45dc-beec-e2ba92e22626 #+BEGIN_SRC jupyter-python :exports both :results raw drawer sherpa = 0.0538009 xs_pb/sherpa #+END_SRC #+RESULTS: 81b5ed93-0312-45dc-beec-e2ba92e22626 : 0.9998585425137037 I had to set the runcard option ~EW_SCHEME: alpha0~ to use the pure QED coupling constant. *** Numerical Integration Plot our nice distribution: #+begin_src jupyter-python :exports both :results raw drawer plot_points = np.linspace(*plot_interval, 1000) fig, ax = set_up_plot() ax.plot(plot_points, gev_to_pb(diff_xs(plot_points, charge=charge, esp=esp))) ax.set_xlabel(r'$\theta$') ax.set_ylabel(r'$\frac{d\sigma}{d\Omega}$ [pb]') ax.axvline(interval[0], color='gray', linestyle='--') ax.axvline(interval[1], color='gray', linestyle='--', label=rf'$|\eta|={η}$') ax.legend() save_fig(fig, 'diff_xs', 'xs', size=[4, 4]) #+end_src #+RESULTS: [[file:./.ob-jupyter/d30ededeaa03958fae5b649f50f3c5c3e6ae4677.png]] Define the integrand. #+begin_src jupyter-python :exports both :results raw drawer def xs_pb_int(θ): return gev_to_pb(np.sin(θ)*diff_xs(θ, charge=charge, esp=esp)) #+end_src #+RESULTS: Plot the integrand. # TODO: remove duplication #+begin_src jupyter-python :exports both :results raw drawer fig, ax = set_up_plot() ax.plot(plot_points, xs_pb_int(plot_points)) ax.set_xlabel(r'$\theta$') ax.set_ylabel(r'$\sin(\theta)\cdot\frac{d\sigma}{d\theta}$ [pb]') ax.axvline(interval[0], color='gray', linestyle='--') ax.axvline(interval[1], color='gray', linestyle='--', label=rf'$|\eta|={η}$') ax.legend() save_fig(fig, 'xs_integrand', 'xs', size=[4, 4]) #+end_src #+RESULTS: [[file:./.ob-jupyter/78974a2e2315c72bd7ae8e4ac009b3d79cfe7001.png]] Intergrate σ with the mc method. #+begin_src jupyter-python :exports both :results raw drawer xs_pb_mc, xs_pb_mc_err = monte_carlo.integrate(xs_pb_int, interval, 10000) xs_pb_mc = xs_pb_mc*np.pi*2 xs_pb_mc, xs_pb_mc_err #+end_src #+RESULTS: | 0.05367831915813582 | 4.264036796520297e-05 | We gonna export that as tex. #+begin_src jupyter-python :exports both :results raw drawer tex_value(xs_pb_mc, unit=r'\pico\barn', prefix=r'\sigma = ', prec=5, save=('results', 'xs_mc.tex')) #+end_src #+RESULTS: : \(\sigma = \SI{0.05368}{\pico\barn}\) *** Sampling and Analysis Define the sample number. #+begin_src jupyter-python :exports both :results raw drawer sample_num = 1000 #+end_src #+RESULTS: Let's define a shortcut for our distribution. #+begin_src jupyter-python :exports both :results raw drawer def dist(x): return gev_to_pb(diff_xs_cosθ(x, charge, esp))*2*np.pi #+end_src #+RESULTS: Now we monte-carlo sample our distribution. We observe that the efficiency his very bad! #+begin_src jupyter-python :exports both :results raw drawer cosθ_sample, cosθ_efficiency = \ monte_carlo.sample_unweighted_array(sample_num, dist, interval_cosθ, report_efficiency=True) cosθ_efficiency #+end_src #+RESULTS: : 0.027318377817804704 Our distribution has a lot of variance, as can be seen by plotting it. #+begin_src jupyter-python :exports both :results raw drawer pts = np.linspace(*interval_cosθ, 100) fig, ax = set_up_plot() ax.plot(pts, dist(pts), label=r'$\frac{d\sigma}{d\Omega}$') #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/04d0c9300d134c04b087aef7bb0a1b6036038b64.png]] :END: We define a friendly and easy to integrate upper limit function. #+begin_src jupyter-python :exports both :results raw drawer upper_limit = dist(interval_cosθ[0]) \ /interval_cosθ[0]**2 upper_base = dist(0) def upper(x): return upper_base + upper_limit*x**2 def upper_int(x): return upper_base*x + upper_limit*x**3/3 ax.plot(pts, upper(pts), label='Upper bound') ax.legend() ax.set_xlabel(r'$\cos\theta$') ax.set_ylabel(r'$\frac{d\sigma}{d\Omega}$') save_fig(fig, 'upper_bound', 'xs_sampling', size=(4, 4)) fig #+end_src #+RESULTS: [[file:./.ob-jupyter/1a720f93049e88987bdddac861b1c3847501e271.png]] To increase our efficiency, we have to specify an upper bound. That is at least a little bit better. The numeric inversion is horribly inefficent. #+begin_src jupyter-python :exports both :results raw drawer cosθ_sample, cosθ_efficiency = \ monte_carlo.sample_unweighted_array(sample_num, dist, interval_cosθ, report_efficiency=True, upper_bound=[upper, upper_int]) cosθ_efficiency #+end_src #+RESULTS: : 0.07584353538369662 Nice! And now draw some histograms. We define an auxilliary method for convenience. #+begin_src jupyter-python :exports both :results raw drawer def draw_histo(points, xlabel, bins=20): fig, ax = set_up_plot() ax.hist(points, bins, histtype='step') ax.set_xlabel(xlabel) ax.set_xlim([points.min(), points.max()]) return fig, ax #+end_src #+RESULTS: The histogram for cosθ. #+begin_src jupyter-python :exports both :results raw drawer fig, _ = draw_histo(cosθ_sample, r'$\cos\theta$') save_fig(fig, 'histo_cos_theta', 'xs', size=(4,3)) #+end_src #+RESULTS: [[file:./.ob-jupyter/6c0c8783b68b394bd8539ce621e285b54f9a0d0c.png]] Now we define some utilities to draw real 4-impulse samples. #+begin_src jupyter-python :exports both :tangle tangled/xs.py def sample_impulses(sample_num, interval, charge, esp, seed=None): """Samples `sample_num` unweighted photon 4-impulses from the cross-section. :param sample_num: number of samples to take :param interval: cosθ interval to sample from :param charge: the charge of the quark :param esp: center of mass energy :param seed: the seed for the rng, optional, default is system time :returns: an array of 4 photon impulses :rtype: np.ndarray """ cosθ_sample = \ monte_carlo.sample_unweighted_array(sample_num, lambda x: diff_xs_cosθ(x, charge, esp), interval_cosθ) φ_sample = np.random.uniform(0, 1, sample_num) def make_impulse(esp, cosθ, φ): sinθ = np.sqrt(1-cosθ**2) return np.array([1, sinθ*np.cos(φ), sinθ*np.sin(φ), cosθ])*esp/2 impulses = np.array([make_impulse(esp, cosθ, φ) \ for cosθ, φ in np.array([cosθ_sample, φ_sample]).T]) return impulses #+end_src #+RESULTS: To generate histograms of other obeservables, we have to define them as functions on 4-impuleses. Using those to transform samples is analogous to transforming the distribution itself. #+begin_src jupyter-python :exports both :results raw drawer :tangle tangled/observables.py """This module defines some observables on arrays of 4-pulses.""" import numpy as np def p_t(p): """Transverse impulse :param p: array of 4-impulses """ return np.linalg.norm(p[:,1:3], axis=1) def η(p): """Pseudo rapidity. :param p: array of 4-impulses """ return np.arccosh(np.linalg.norm(p[:,1:], axis=1)/p_t(p))*np.sign(p[:, 3]) #+end_src #+RESULTS: Lets try it out. #+begin_src jupyter-python :exports both :results raw drawer impulse_sample = sample_impulses(2000, interval_cosθ, charge, esp) impulse_sample #+end_src #+RESULTS: : array([[100. , 28.29533668, 1.8580239 , -95.89536834], : [100. , 15.63107713, 8.96259547, 98.36331283], : [100. , 37.20714833, 36.40572302, -85.38296929], : ..., : [100. , 23.83872893, 21.09208613, 94.79893937], : [100. , 53.08852778, 37.82714439, 75.83347127], : [100. , 78.42901358, 58.54310277, 20.5327774 ]]) Now let's make a histogram of the η distribution. #+begin_src jupyter-python :exports both :results raw drawer η_sample = η(impulse_sample) draw_histo(η_sample, r'$\eta$') #+end_src #+RESULTS: :RESULTS: |
| | [[file:./.ob-jupyter/e93174871d533af3abc0a74dc6e4641aefc138ad.png]] :END: And the same for the p_t (transverse impulse) distribution. #+begin_src jupyter-python :exports both :results raw drawer p_t_sample = p_t(impulse_sample) draw_histo(p_t_sample, r'$p_T$') #+end_src #+RESULTS: :RESULTS: |
| | [[file:./.ob-jupyter/e63cd849824f57bdcad1ff45ef940def0f5c5934.png]] :END: