bachelor_thesis/prog/python/qqgg/analytical_xs.org

Investigaton of Monte-Carlo Methods

• Init
  • Required Modules
  • Utilities
• Implementation
  • Center of Mass Frame
• Calculations
  • Analytical Integration
  • Numerical Integration
    • Integral over θ
    • Integration over η
    • Using VEGAS
    • Testing the Statistics
  • Sampling and Analysis
    • Sampling the cosθ cross section
    • Observables
    • Sampling the η cross section
    • Sampling with VEGAS
    • Some Histograms with Rivet
      • Init
      • Plot the Histos

Init

Required Modules

  import numpy as np
  import matplotlib.pyplot as plt
  import monte_carlo

Utilities

%run ../utility.py
%load_ext autoreload
%aimport monte_carlo
%autoreload 1

Implementation

Center of Mass Frame

  """
  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 * (np.tanh(η) ** 2 + 1)


  def diff_xs_p_t(p_t, charge, esp):
      """
      Calculates the differential cross section as a function of the
      transverse momentum (p_t) of the photons in units of 1/GeV^2.

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

      Arguments:
      p_t -- transverse momentum in GeV
      esp -- center of momentum energy in GeV
      charge -- charge of the particle in units of the elementary charge
      """

      f = energy_factor(charge, esp)
      sqrt_fact = np.sqrt(1 - (2 * p_t / esp) ** 2)
      return f / p_t * (1 / sqrt_fact + sqrt_fact)


  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]))

Analytical Integration

And now calculate the cross section in picobarn.

  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'))

\(\sigma = \SI{0.053793}{\pico\barn}\)

Lets plot the total xs as a function of η.

  fig, ax = set_up_plot()
  η_s = np.linspace(0, 3, 1000)
  ax.plot(η_s, gev_to_pb(total_xs_eta(η_s, charge, esp)))
  ax.set_xlabel(r'$\eta$')
  ax.set_ylabel(r'$\sigma$ [pb]')
  ax.set_xlim([0, max(η_s)])
  ax.set_ylim(0)
  save_fig(fig, 'total_xs', 'xs', size=[2.5, 2.5])

Compared to sherpa, it's pretty close.

  sherpa = 0.05380
  xs_pb - sherpa

-6.7112594623469635e-06

I had to set the runcard option EW_SCHEME: alpha0 to use the pure QED coupling constant.

Numerical Integration

Plot our nice distribution:

  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'$d\sigma/d\Omega$ [pb]')
  ax.set_xlim([plot_points.min(), plot_points.max()])
  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=[2.5, 2.5])

Define the integrand.

  def xs_pb_int(θ):
      return 2*np.pi*gev_to_pb(np.sin(θ)*diff_xs(θ, charge=charge, esp=esp))

  def xs_pb_int_η(η):
      return 2*np.pi*gev_to_pb(diff_xs_eta(η, charge, esp))

Plot the integrand. # TODO: remove duplication

  fig, ax = set_up_plot()
  ax.plot(plot_points, xs_pb_int(plot_points))
  ax.set_xlabel(r'$\theta$')
  ax.set_ylabel(r'$2\pi\cdot d\sigma/d\theta [pb]')
  ax.set_xlim([plot_points.min(), plot_points.max()])
  ax.axvline(interval[0], color='gray', linestyle='--')
  ax.axvline(interval[1], color='gray', linestyle='--', label=rf'$|\eta|={η}$')
  save_fig(fig, 'xs_integrand', 'xs', size=[3, 2.2])

Integral over θ

Intergrate σ with the mc method.

  xs_pb_res = monte_carlo.integrate(xs_pb_int, interval, epsilon=1e-3)
  xs_pb_res

IntegrationResult(result=0.054374269511586644, sigma=0.0009683114962275483, N=2123)

We gonna export that as tex.

  tex_value(*xs_pb_res.combined_result, unit=r'\pico\barn',
            prefix=r'\sigma = ', save=('results', 'xs_mc.tex'))
  tex_value(xs_pb_res.N, prefix=r'N = ', save=('results', 'xs_mc_N.tex'))

\(N = 2123\)

Integration over η

Plot the intgrand of the pseudo rap.

  fig, ax = set_up_plot()
  points = np.linspace(-4, 4, 1000)
  ax.set_xlim([-4, 4])
  ax.plot(points, xs_pb_int_η(points))
  ax.set_xlabel(r'$\eta$')
  ax.set_ylabel(r'$2\pi\cdot d\sigma/d\eta$ [pb]')
  ax.axvline(interval_η[0], color='gray', linestyle='--')
  ax.axvline(interval_η[1], color='gray', linestyle='--', label=rf'$|\eta|={η}$')
  save_fig(fig, 'xs_integrand_eta', 'xs', size=[3, 2])

  xs_pb_η = monte_carlo.integrate(xs_pb_int_η,
                                  interval_η, epsilon=1e-3)
  xs_pb_η

IntegrationResult(result=0.055151085713932124, sigma=0.000932443763615135, N=140)

As we see, the result is a little better if we use pseudo rapidity, because the differential cross section does not difverge anymore. But becase our η interval is covering the range where all the variance is occuring, the improvement is rather marginal.

And yet again export that as tex.

  tex_value(*xs_pb_η.combined_result, unit=r'\pico\barn', prefix=r'\sigma = ',
            save=('results', 'xs_mc_eta.tex'))
  tex_value(xs_pb_η.N, prefix=r'N = ', save=('results', 'xs_mc_eta_N.tex'))

\(N = 140\)

Using VEGAS

Now we use VEGAS on the θ parametrisation and see what happens.

  num_increments = 11
  xs_pb_vegas = monte_carlo.integrate_vegas(
      xs_pb_int,
      interval,
      num_increments=num_increments,
      alpha=1,
      epsilon=1e-3,
      acumulate=False,
      vegas_point_density=100,
  )
  xs_pb_vegas

VegasIntegrationResult(result=0.05367272259033393, sigma=0.00042389745443520844, N=275, increment_borders=array([0.16380276, 0.25021069, 0.37424685, 0.55594852, 0.83854962,
       1.27086143, 1.85730705, 2.29700734, 2.57637351, 2.76334644,
       2.89011895, 2.9777899 ]), vegas_iterations=12)

This is pretty good, although the variance reduction may be achieved partially by accumulating the results from all runns. Here this gives us one order of magnitude more than we wanted.

And export that as tex.

  tex_value(*xs_pb_vegas.combined_result, unit=r'\pico\barn',
            prefix=r'\sigma = ', save=('results', 'xs_mc_θ_vegas.tex'))
  tex_value(xs_pb_vegas.N, prefix=r'N = ', save=('results', 'xs_mc_θ_vegas_N.tex'))
  tex_value(num_increments, prefix=r'K = ', save=('results', 'xs_mc_θ_vegas_K.tex'))

\(K = 11\)

Surprisingly, acumulation, the result ain't much different. This depends, of course, on the iteration count.

  monte_carlo.integrate_vegas(
      xs_pb_int,
      interval,
      num_increments=num_increments,
      alpha=1,
      epsilon=1e-3,
      acumulate=True,
      vegas_point_density=100,
  )

VegasIntegrationResult(result=0.05405856822809377, sigma=0.0003046830887499666, N=275, increment_borders=array([0.16380276, 0.25854516, 0.39088443, 0.58108696, 0.86570559,
       1.29975004, 1.8564745 , 2.28724136, 2.56829826, 2.75438203,
       2.88442585, 2.9777899 ]), vegas_iterations=10)

Let's define some little helpers.

  """
  Some shorthands for common plotting tasks related to the investigation
  of monte-carlo methods in one rimension.

  Author: Valentin Boettcher <hiro at protagon.space>
  """

  import matplotlib.pyplot as plt
  import numpy as np
  from utility import *


  def plot_increments(ax, increment_borders, label=None, *args, **kwargs):
      """Plot the increment borders from a list.  The first and last one

      :param ax: the axis on which to draw
      :param list increment_borders: the borders of the increments
      :param str label: the label to apply to one of the vertical lines
      """

      ax.axvline(x=increment_borders[1], label=label, *args, **kwargs)

      for increment in increment_borders[1:-1]:
          ax.axvline(x=increment, *args, **kwargs)


  def plot_vegas_weighted_distribution(
      ax, points, dist, increment_borders, *args, **kwargs
  ):
      """Plot the distribution with VEGAS weights applied.

      :param ax: axis
      :param points: points
      :param dist: distribution
      :param increment_borders: increment borders
      """

      num_increments = increment_borders.size
      weighted_dist = dist.copy()

      for left_border, right_border in zip(increment_borders[:-1], increment_borders[1:]):
          length = right_border - left_border
          mask = (left_border <= points) & (points <= right_border)
          weighted_dist[mask] = dist[mask] * num_increments * length

      ax.plot(points, weighted_dist, *args, **kwargs)


  def plot_stratified_rho(ax, points, increment_borders, *args, **kwargs):
      """Plot the weighting distribution resulting from the increment
      borders.

      :param ax: axis
      :param points: points
      :param increment_borders: increment borders

      """

      num_increments = increment_borders.size
      ρ = np.empty_like(points)
      for left_border, right_border in zip(increment_borders[:-1], increment_borders[1:]):
          length = right_border - left_border
          mask = (left_border <= points) & (points <= right_border)
          ρ[mask] = 1 / (num_increments * length)

      ax.plot(points, ρ, *args, **kwargs)

And now we plot the integrand with the incremens.

  fig, ax = set_up_plot()
  ax.set_xlim(*interval)
  ax.set_xlabel(r"$\theta$")
  ax.set_ylabel(r"$2\pi\cdot d\sigma/d\theta$ [pb]")
  ax.set_ylim([0, 0.09])

  ax.plot(plot_points, xs_pb_int(plot_points), label="Distribution")

  plot_increments(
      ax,
      xs_pb_vegas.increment_borders,
      label="Increment Borders",
      color="gray",
      linestyle="--",
  )

  plot_vegas_weighted_distribution(
      ax,
      plot_points,
      xs_pb_int(plot_points),
      xs_pb_vegas.increment_borders,
      label="Weighted Distribution",
  )

  ax.legend(fontsize="small", loc="lower left")
  save_fig(fig, "xs_integrand_vegas", "xs", size=[5, 3])

Testing the Statistics

Let's battle test the statistics.

  num_runs = 1000
  num_within = 0

  for _ in range(num_runs):
      val, err = \
          monte_carlo.integrate(xs_pb_int, interval, epsilon=1e-3).combined_result
      if abs(xs_pb - val) <= err:
          num_within += 1

  num_within/num_runs

0.709

So we see: the standard deviation is sound.

Doing the same thing with VEGAS works as well.

  num_runs = 1000
  num_within = 0
  for _ in range(num_runs):
      val, err = \
          monte_carlo.integrate_vegas(xs_pb_int, interval,
                                      num_increments=10, alpha=1,
                                      epsilon=1e-3, acumulate=False,
                                      vegas_point_density=100).combined_result

      if abs(xs_pb - val) <= err:
          num_within += 1
  num_within/num_runs

0.66

Sampling and Analysis

Define the sample number.

  sample_num = 10000
  tex_value(
      sample_num, prefix="N = ", save=("results", "4imp-sample-size.tex"),
  )

\(N = 10000\)

Let's define shortcuts for our distributions. The 2π are just there for formal correctnes. Factors do not influecence the outcome.

  def dist_cosθ(x):
      return gev_to_pb(diff_xs_cosθ(x, charge, esp))

  def dist_η(x):
      return gev_to_pb(diff_xs_eta(x, charge, esp))

Sampling the cosθ cross section

Now we monte-carlo sample our distribution. We observe that the efficiency his very bad!

  cosθ_sample, cosθ_efficiency = \
      monte_carlo.sample_unweighted_array(sample_num, dist_cosθ,
                                          interval_cosθ, report_efficiency=True)
  cosθ_efficiency

0.027612227729295263

Let's save that.

  tex_value(
      cosθ_efficiency * 100,
      prefix=r"\mathfrak{e} = ",
      suffix=r"\%",
      save=("results", "naive_th_samp.tex"),
  )

\(\mathfrak{e} = 3\%\)

Our distribution has a lot of variance, as can be seen by plotting it.

  pts = np.linspace(*interval_cosθ, 100)
  fig, ax = set_up_plot()
  ax.plot(pts, dist_cosθ(pts))
  ax.set_xlabel(r'$\cos\theta$')
  ax.set_ylabel(r'$\frac{d\sigma}{d\Omega}$')

Text(0, 0.5, '$\\frac{d\\sigma}{d\\Omega}$')

We define a friendly and easy to integrate upper limit function.

  fig, ax = set_up_plot()
  upper_limit = dist_cosθ(interval_cosθ[0]) / interval_cosθ[0] ** 2
  upper_base = dist_cosθ(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.plot(pts, dist_cosθ(pts), label=r"$f_{\cos\theta}$")

  ax.legend(fontsize='small')
  ax.set_xlabel(r"$\cos\theta$")
  ax.set_ylabel(r"$\frac{d\sigma}{d\cos\theta}$ [pb]")
  save_fig(fig, "upper_bound", "xs_sampling", size=(3, 2.5))

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.

  cosθ_sample_tuned, cosθ_efficiency_tuned = \
      monte_carlo.sample_unweighted_array(sample_num, dist_cosθ,
                                          interval_cosθ, report_efficiency=True,
                                          upper_bound=[upper, upper_int])
  cosθ_efficiency_tuned

0.07879506593284073

<<cosθ-bare-eff>>

  tex_value(
      cosθ_efficiency_tuned * 100,
      prefix=r"\mathfrak{e} = ",
      suffix=r"\%",
      save=("results", "tuned_th_samp.tex"),
  )

\(\mathfrak{e} = 8\%\)

Nice! And now draw some histograms.

We define an auxilliary method for convenience.

  def draw_histo(points, xlabel, bins=50, range=None, **kwargs):
      heights, edges = np.histogram(points, bins, range=range, **kwargs)
      centers = (edges[1:] + edges[:-1]) / 2
      deviations = np.sqrt(heights)
      integral = heights @ (edges[1:] - edges[:-1])
      heights = heights / integral
      deviations = deviations / integral

      fig, ax = set_up_plot()
      ax.errorbar(centers, heights, deviations, linestyle="none", color="orange")
      ax.step(edges, [heights[0], *heights], color="#1f77b4")

      ax.set_xlabel(xlabel)
      ax.set_ylabel("Count")
      ax.set_xlim(range if range is not None else [points.min(), points.max()])
      return fig, ax

The histogram for cosθ.

fig, _ = draw_histo(cosθ_sample, r'$\cos\theta$')
save_fig(fig, 'histo_cos_theta', 'xs', size=(4,3))

Observables

Now we define some utilities to draw real 4-momentum samples.

  def sample_momenta(sample_num, interval, charge, esp, seed=None):
      """Samples `sample_num` unweighted photon 4-momenta 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 momenta

      :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_momentum(esp, cosθ, φ):
          sinθ = np.sqrt(1-cosθ**2)
          return np.array([1, sinθ*np.cos(φ), sinθ*np.sin(φ), cosθ])*esp/2

      momenta = np.array([make_momentum(esp, cosθ, φ) \
                           for cosθ, φ in np.array([cosθ_sample, φ_sample]).T])
      return momenta

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.

  """This module defines some observables on arrays of 4-pulses."""
  import numpy as np

  def p_t(p):
      """Transverse momentum

      :param p: array of 4-momenta
      """

      return np.linalg.norm(p[:,1:3], axis=1)

  def η(p):
      """Pseudo rapidity.

      :param p: array of 4-momenta
      """

      return np.arccosh(np.linalg.norm(p[:,1:], axis=1)/p_t(p))*np.sign(p[:, 3])

And import them.

  %aimport tangled.observables
  obs = tangled.observables

Lets try it out.

  momentum_sample = sample_momenta(sample_num, interval_cosθ, charge, esp)
  momentum_sample

array([[100.        ,  36.96504312,  56.43920882,  73.81193193],
       [100.        ,  47.94386939,   1.7718024 , -87.73964955],
       [100.        ,  43.23224845,   8.595659  , -89.76127974],
       ...,
       [100.        ,  40.84435458,  21.7798789 , -88.64183873],
       [100.        ,  64.07275183,   2.13440481, -76.74716144],
       [100.        ,  74.77621972,   2.38367531, -66.35386241]])

Now let's make a histogram of the η distribution.

  η_sample = obs.η(momentum_sample)
  fig, ax = draw_histo(η_sample, r'$\eta$', range=interval_η)
  save_fig(fig, 'histo_eta', 'xs_sampling', size=[3, 3])

And the same for the p_t (transverse momentum) distribution.

  p_t_sample = obs.p_t(momentum_sample)
  fig, ax = draw_histo(p_t_sample, r"$p_T$ [GeV]", range=interval_pt)
  save_fig(fig, "histo_pt", "xs_sampling", size=[3, 3])

That looks somewhat fishy, but it isn't.

  fig, ax = set_up_plot()
  points = np.linspace(interval_pt[0], interval_pt[1] - .01, 1000)
  ax.plot(points, gev_to_pb(diff_xs_p_t(points, charge, esp)))
  ax.set_xlabel(r'$p_T$')
  ax.set_xlim(interval_pt[0], interval_pt[1] + 1)
  ax.set_ylim([0, gev_to_pb(diff_xs_p_t(interval_pt[1] -.01, charge, esp))])
  ax.set_ylabel(r'$\frac{d\sigma}{dp_t}$ [pb]')
  save_fig(fig, 'diff_xs_p_t', 'xs_sampling', size=[4, 2])

this is strongly peaked at p_t=100GeV. (The jacobian goes like 1/x there!)

Sampling the η cross section

An again we see that the efficiency is way, way! better…

  η_sample, η_efficiency = monte_carlo.sample_unweighted_array(
      sample_num, dist_η, interval_η, report_efficiency=True
  )
  tex_value(
      η_efficiency * 100,
      prefix=r"\mathfrak{e} = ",
      suffix=r"\%",
      save=("results", "eta_eff.tex"),
  )

\(\mathfrak{e} = 41\%\)

<<η-eff>>

Let's draw a histogram to compare with the previous results.

  draw_histo(η_sample, r'$\eta$')

<Figure

size

432x288

with

1

Axes>

<matplotlib.axes._subplots.AxesSubplot

at

0x7f048b5dcd30>

Looks good to me :).

Sampling with VEGAS

To get the increments, we have to let VEGAS loose on our distribution. We throw away the integral, but keep the increments.

  K = 10
  increments = monte_carlo.integrate_vegas(
      dist_cosθ, interval_cosθ, num_increments=K, alpha=1, increment_epsilon=0.001
  ).increment_borders
  tex_value(
      K, prefix=r"K = ", save=("results", "vegas_samp_num_increments.tex"),
  )
  increments

array([-0.9866143 , -0.96987487, -0.93118512, -0.8386111 , -0.60462504,
       -0.00295135,  0.59955266,  0.83691601,  0.93048622,  0.96936191,
        0.9866143 ])

Visualizing the increment borders gives us the information we want.

  pts = np.linspace(*interval_cosθ, 100)
  fig, ax = set_up_plot()
  ax.plot(pts, dist_cosθ(pts))
  ax.set_xlabel(r'$\cos\theta$')
  ax.set_ylabel(r'$\frac{d\sigma}{d\Omega}$')
  ax.set_xlim(*interval_cosθ)
  plot_increments(ax, increments,
                  label='Increment Borderds', color='gray', linestyle='--')
  ax.legend()

<matplotlib.legend.Legend at 0x7f048b526550>

We can now plot the reweighted distribution to observe the variance reduction visually.

  pts = np.linspace(*interval_cosθ, 1000)
  fig, ax = set_up_plot()
  ax.plot(pts, dist_cosθ(pts), label="Distribution")
  plot_vegas_weighted_distribution(
      ax, pts, dist_cosθ(pts), increments, label="Weighted Distribution"
  )
  ax.set_xlabel(r"$\cos\theta$")
  ax.set_ylabel(r"$\frac{d\sigma}{d\cos\theta}$")
  ax.set_xlim(*interval_cosθ)
  plot_increments(
      ax, increments, label="Increment Borderds", color="gray", linestyle="--"
  )
  ax.legend(fontsize="small")
  save_fig(fig, "vegas_strat_dist", "xs_sampling", size=(3, 2.3))

I am batman! Let's plot the weighting distribution.

  pts = np.linspace(*interval_cosθ, 1000)
  fig, ax = set_up_plot()
  plot_stratified_rho(ax, pts, increments)
  ax.set_xlabel(r"$\cos\theta$")
  ax.set_ylabel(r"$\rho")
  ax.set_xlim(*interval_cosθ)
  save_fig(fig, "vegas_rho", "xs_sampling", size=(3, 2.3))

Now, draw a sample and look at the efficiency.

  cosθ_sample_strat, cosθ_efficiency_strat = \
      monte_carlo.sample_unweighted_array(sample_num, dist_cosθ,
                                          increment_borders=increments,
                                          report_efficiency=True)
  cosθ_efficiency_strat

0.59549

  tex_value(
      cosθ_efficiency_strat * 100,
      prefix=r"\mathfrak{e} = ",
      suffix=r"\%",
      save=("results", "strat_th_samp.tex"),
  )

\(\mathfrak{e} = 60\%\)

If we compare that to /hiro/bachelor_thesis/src/commit/160d9790c33ed7b609a7ae47ef752371f17c663d/prog/python/qqgg/cos%CE%B8-bare-eff, we can see the improvement :P. It is even better the /hiro/bachelor_thesis/src/commit/160d9790c33ed7b609a7ae47ef752371f17c663d/prog/python/qqgg/%CE%B7-eff. The histogram looks just the same.

fig, _ = draw_histo(cosθ_sample_strat, r'$\cos\theta$')
save_fig(fig, 'histo_cos_theta_strat', 'xs', size=(4,3))

Some Histograms with Rivet

Init

  import yoda

Plot the Histos

  def draw_yoda_histo(h, xlabel):
      edges = np.append(h.xMins(), h.xMax())
      heights = np.append(h.yVals(), h.yVals()[-1])
      centers = (edges[1:] + edges[:-1]) / 2

      fig, ax = set_up_plot()
      ax.errorbar(h.xVals(), h.yVals(), h.yErrs(), linestyle="none", color="orange")
      ax.step(edges, heights, color="#1f77b4", where="post")

      ax.set_xlabel(xlabel)
      ax.set_ylabel("Count")
      ax.set_xlim([h.xMin(), h.xMax()])
      return fig, ax

  yoda_file = yoda.read("../../runcards/qqgg/analysis/Analysis.yoda")
  sherpa_histos = {"pT": r"$p_T$ [GeV]", "eta": r"$\eta$", "cos_theta": r"$\cos\theta$"}

  for key, label in sherpa_histos.items():
      fig, ax = draw_yoda_histo(
          yoda_file["/MC_DIPHOTON_SIMPLE/" + key], r"Sherpa " + label
      )
      save_fig(fig, "histo_sherpa_" + key, "xs_sampling", size=(3, 3))