bachelor_thesis/prog/python/qqgg/analytical_xs.org

#+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 <hiro@protagon.space>
  """

  import numpy as np

  # 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*/(137.036*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².

      Here dΩ=sinθdθdφ

      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².

      Here dΩ=d(cosθ)dφ

      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.

      This is actually the crossection dσ/(dφdη).

      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)*2*np.exp(-η)/np.cosh(η)**2

  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
:RESULTS:
# [goto error]
:   File "<ipython-input-39-3fe28bd90183>", line 20
:     return charge**4*/(137.036*esp)**2/6
:                      ^
: SyntaxError: invalid syntax
:END:

* 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=6, save=('results', 'xs.tex'))
#+END_SRC

#+RESULTS: cf853fb6-d338-482e-bc55-bd9f8e796495
: \(\sigma = \SI{0.053793}{\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.05380
  xs_pb - sherpa
#+END_SRC

#+RESULTS: 81b5ed93-0312-45dc-beec-e2ba92e22626
: -6.710540074485183e-06

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=[3, 3])
#+end_src

#+RESULTS:
:RESULTS:
: /usr/lib/python3.8/site-packages/tikzplotlib/_axes.py:211: MatplotlibDeprecationWarning: Passing the minor parameter of get_xticks() positionally is deprecated since Matplotlib 3.2; the parameter will become keyword-only two minor releases later.
:   data, "minor x", obj.get_xticks("minor"), obj.get_xticklabels("minor")
: /usr/lib/python3.8/site-packages/tikzplotlib/_axes.py:216: MatplotlibDeprecationWarning: Passing the minor parameter of get_yticks() positionally is deprecated since Matplotlib 3.2; the parameter will become keyword-only two minor releases later.
:   data, "minor y", obj.get_yticks("minor"), obj.get_yticklabels("minor")
[[file:./.ob-jupyter/84554de6897392b423848ccff74c3b1bdbbac799.png]]
:END:

Define the integrand.
#+begin_src jupyter-python :exports both :results raw drawer
  def xs_pb_int(θ):
      return 2*np.pi*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/9e547bdeaa79bb956057b552090b4ab6305a20e6.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
  xs_pb_mc, xs_pb_mc_err
#+end_src

#+RESULTS:
| 0.053599995094203025 | 0.0002656256326580591 |

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 = ', err=xs_pb_mc_err, save=('results', 'xs_mc.tex'))
#+end_src

#+RESULTS:
: \(\sigma = \SI{0.05360\pm 0.00027}{\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.026054076452916193

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:
| <matplotlib.lines.Line2D | at | 0x7f3dc1880fd0> |
[[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.07947335025380711

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/d473cf9d22d8fe203293e6d17a92497aad3109d3.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.        ,  16.36721437,   5.49983339,  98.49805138],
:        [100.        ,  58.57273425,  71.41789832, -38.32386465],
:        [100.        ,  73.44354101,  23.28263462, -63.74923693],
:        ...,
:        [100.        ,  86.35169559,  12.11748171,  48.95458411],
:        [100.        ,  58.83596982,   7.2563454 , -80.53368306],
:        [100.        ,  55.49634462,  66.91946554, -49.41599807]])

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:
| <Figure | size | 432x288 | with | 1 | Axes> | <matplotlib.axes._subplots.AxesSubplot | at | 0x7f3dc20caf10> |
[[file:./.ob-jupyter/70b0373104835c034a1444836807c357c2f8aacb.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$ [GeV]')
#+end_src

#+RESULTS:
:RESULTS:
| <Figure | size | 432x288 | with | 1 | Axes> | <matplotlib.axes._subplots.AxesSubplot | at | 0x7f3dc17120a0> |
[[file:./.ob-jupyter/0329079132169b82385536d1b707f5fe884f70aa.png]]
:END: