something is still wrong with the PDF

prog/python/qqgg/monte_carlo.py

@@ -43,9 +43,7 @@ class IntegrationResult:
         return self.result, self.sigma
 
 
-def integrate(
-    f, interval, epsilon=0.01, seed=None, cut_fn=None, **kwargs
-) -> IntegrationResult:
+def integrate(f, interval, epsilon=0.01, seed=None, **kwargs) -> IntegrationResult:
     """Monte-Carlo integrates the function `f` over an interval.
 
     :param f: function of one variable, kwargs are passed to it

prog/python/qqgg/parton_density_function_stuff.org

@@ -22,17 +22,15 @@
 #+END_SRC
 
 #+RESULTS:
-: The autoreload extension is already loaded. To reload it, use:
-:   %reload_ext autoreload
 
 ** Global Config
 #+begin_src jupyter-python :exports both :results raw drawer
 η = 2.4
-e_proton = 100  # GeV
+e_proton = 6500  # GeV
 interval_η = [-η, η]
 interval = η_to_θ([-η, η])
 interval_cosθ = np.cos(interval)
-num_samples = 3_000
+num_samples = 10_000
 #+end_src
 
 #+RESULTS:
@@ -172,6 +170,18 @@ We begin by implementing the same sermon for the lab frame.
   @vectorize([float64(float64, float64, float64)], nopython=True)
   def averaged_tchanel_q2(e_proton, x_1, x_2):
       return 2 * x_1 * x_2 * e_proton ** 2
+
+
+  def cut_pT_from_eta(greater_than=0):
+      def cut(e_proton, η, x1, x2):
+          cosθ = np.cos(η_to_θ(η))
+          _, _, p1, p2 = momenta(e_proton, x1, x2, cosθ)
+          return (
+              np.sqrt((p1[0][1:3] ** 2).sum()) > greater_than
+              and np.sqrt((p2[0][1:3] ** 2).sum()) > greater_than
+          )
+
+      return cut
 #+end_src
 
 #+RESULTS:
@@ -181,7 +191,7 @@ We begin by implementing the same sermon for the lab frame.
   from numba.extending import get_cython_function_address
 
 
-  def get_xs_distribution_with_pdf(xs, q, e_hadron, quarks=None, pdf=None):
+  def get_xs_distribution_with_pdf(xs, q, e_hadron, quarks=None, pdf=None, cut=None):
       """Creates a function that takes an event (type np.ndarray) of the
       form [cosθ, impulse fractions of quarks in hadron 1, impulse
       fractions of quarks in hadron 2] and returns the differential
@@ -196,9 +206,11 @@ We begin by implementing the same sermon for the lab frame.
       :param quarks: the constituent quarks np.ndarray of the form [[id, charge], ...],
       the default is a proton
       :param pdf: the PDF to use, the default is "NNPDF31_lo_as_0118"
+      :param cut: cut function with signature (energy hadron, cosθ, x_1,
+      x_2) to return 0, when the event does not fit the cut
+
       :returns: differential cross section summed over flavors and weighted with the pdfs
       :rtype: function
-
       """
 
       pdf = pdf or lhapdf.mkPDF("NNPDF31_lo_as_0118", 0)
@@ -213,6 +225,10 @@ We begin by implementing the same sermon for the lab frame.
 
       # @jit(float64(float64[4])) Unfortunately that does not work as yet!
       def distribution(event: np.ndarray) -> float:
+
+          if cut and not cut(e_hadron, *event):
+              return 0
+
           cosθ, x_1, x_2 = event
 
           q2_value = q(e_hadron, x_1, x_2)
@@ -275,8 +291,8 @@ Let's plot it for some random values 😃.
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7fca5ee5dbb0> |
-[[file:./.ob-jupyter/f12b49e327ad8ec823397b3bbe910e0d021f4dcd.png]]
+| <matplotlib.lines.Line2D | at | 0x7f0bd18b0700> |
+[[file:./.ob-jupyter/4b7773815fe4943f422b5943cc67e48b7a75cc23.png]]
 :END:
 
 Having set both x to the same value, we get a symmetric distribution as expected.
@@ -292,7 +308,7 @@ But first we have to find an upper bound, which is expensive!
 #+end_src
 
 #+RESULTS:
-: 2171.468698483163
+: 5721.40648474465
 
 Beware!, this is darn slow, becaus the efficiency is soooo low.
 #+begin_src jupyter-python :exports both :results raw drawer
@@ -309,16 +325,13 @@ Beware!, this is darn slow, becaus the efficiency is soooo low.
 #+end_src
 
 #+RESULTS:
-:RESULTS:
-: sample_unweighted_array
 : 0.0004240744427426794
-:END:
 
 ** Switching Horses: Sampling η
 We set up a new distribution.
 #+begin_src jupyter-python :exports both :results raw drawer
   dist_η, x_limits = get_xs_distribution_with_pdf(
-      diff_xs_η, averaged_tchanel_q2, e_proton
+      diff_xs_η, averaged_tchanel_q2, e_proton, cut=cut_pT_from_eta(greater_than=20)
   )
 #+end_src
 
@@ -336,8 +349,8 @@ Plotting it, we can see that the variance is reduced.
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7fca5b671280> |
-[[file:./.ob-jupyter/8fcbc6b5ce43840c893bc34aff9178d3c9cbf5fe.png]]
+| <matplotlib.lines.Line2D | at | 0x7f0bcf65abe0> |
+[[file:./.ob-jupyter/5f2e010bf22bb8e157d5327258268cf8a465510d.png]]
 :END:
 
 Lets plot how the pdf looks.
@@ -351,7 +364,7 @@ Lets plot how the pdf looks.
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7fca5bb7bd60> |
+| <matplotlib.lines.Line2D | at | 0x7f0bd1831100> |
 [[file:./.ob-jupyter/b92f0c4b2c9f2195ae14444748fcdb7708d81c19.png]]
 :END:
 
@@ -368,19 +381,19 @@ figure out the cpu mapping.
       interval=intervals_η,
       proc="auto",
       report_efficiency=True,
-      #cache="cache/pdf/huge",
+      cache="cache/pdf/huge",
   )
   result
 #+end_src
 
 #+RESULTS:
-: array([[ 1.3399508 ,  0.02307023,  0.06686396],
-:        [ 2.32026227,  0.03462208,  0.42094429],
-:        [-2.33108558,  0.01082153,  0.01427812],
+: array([[-1.43205911,  0.34762794,  0.01802916],
+:        [-2.13582434,  0.02851933,  0.07415878],
+:        [ 1.79962451,  0.01835771,  0.2725591 ],
 :        ...,
-:        [-1.95937336,  0.18571132,  0.0883617 ],
-:        [ 1.0941634 ,  0.20425188,  0.12512649],
-:        [-1.94326146,  0.03546833,  0.16993543]])
+:        [ 1.1782848 ,  0.1113013 ,  0.01554592],
+:        [-1.27942179,  0.12550189,  0.05389335],
+:        [ 2.22491947,  0.04076429,  0.13610809]])
 
 
 
@@ -397,6 +410,6 @@ Let's look at a histogramm of eta samples.
 
 #+RESULTS:
 :RESULTS:
-| <Figure | size | 432x288 | with | 1 | Axes> | <matplotlib.axes._subplots.AxesSubplot | at | 0x7fca5b460f10> |
-[[file:./.ob-jupyter/a09a3d7a5a19938f065660ea42fe47e9ecd04bba.png]]
+| <Figure | size | 432x288 | with | 1 | Axes> | <matplotlib.axes._subplots.AxesSubplot | at | 0x7f0bcd5a1be0> |
+[[file:./.ob-jupyter/b55eed7143c116ac0471a510bd62174a42e1ac31.png]]
 :END:

prog/python/qqgg/tangled/pdf.py

@@ -130,10 +130,22 @@ def diff_xs_η(e_proton, charge, η, x_1, x_2):
 def averaged_tchanel_q2(e_proton, x_1, x_2):
     return 2 * x_1 * x_2 * e_proton ** 2
 
+
+def cut_pT_from_eta(greater_than=0):
+    def cut(e_proton, η, x1, x2):
+        cosθ = np.cos(η_to_θ(η))
+        _, _, p1, p2 = momenta(e_proton, x1, x2, cosθ)
+        return (
+            np.sqrt((p1[0][1:3] ** 2).sum()) > greater_than
+            and np.sqrt((p2[0][1:3] ** 2).sum()) > greater_than
+        )
+
+    return cut
+
 from numba.extending import get_cython_function_address
 
 
-def get_xs_distribution_with_pdf(xs, q, e_hadron, quarks=None, pdf=None):
+def get_xs_distribution_with_pdf(xs, q, e_hadron, quarks=None, pdf=None, cut=None):
     """Creates a function that takes an event (type np.ndarray) of the
     form [cosθ, impulse fractions of quarks in hadron 1, impulse
     fractions of quarks in hadron 2] and returns the differential
@@ -148,9 +160,11 @@ def get_xs_distribution_with_pdf(xs, q, e_hadron, quarks=None, pdf=None):
     :param quarks: the constituent quarks np.ndarray of the form [[id, charge], ...],
     the default is a proton
     :param pdf: the PDF to use, the default is "NNPDF31_lo_as_0118"
+    :param cut: cut function with signature (energy hadron, cosθ, x_1,
+    x_2) to return 0, when the event does not fit the cut
+
     :returns: differential cross section summed over flavors and weighted with the pdfs
     :rtype: function
-
     """
 
     pdf = pdf or lhapdf.mkPDF("NNPDF31_lo_as_0118", 0)
@@ -165,6 +179,10 @@ def get_xs_distribution_with_pdf(xs, q, e_hadron, quarks=None, pdf=None):
 
     # @jit(float64(float64[4])) Unfortunately that does not work as yet!
     def distribution(event: np.ndarray) -> float:
+
+        if cut and not cut(e_hadron, *event):
+            return 0
+
         cosθ, x_1, x_2 = event
 
         q2_value = q(e_hadron, x_1, x_2)

prog/runcards/pp/Sherpa.yaml

@@ -1,8 +1,8 @@
 # general options
 OUTPUT: 3
 
-# no events
-EVENTS: 100000  # TODO: Increase!
+# event count
+EVENTS: 10000  # TODO: Increase!
 
 # save results
 GENERATE_RESULT_DIRECTORY: true
@@ -39,9 +39,10 @@ FRAGMENTATION: None
 SHOWER: None
 MI_HANDLER: None
 
-# cut to $\abs{\eta} \leq 2.5$
+# cut to $\abs{\eta} \leq 2.5$, and $p_T \geq 20$
 SELECTORS:
-- [PT, 22, 20, 6500]
+- [Eta, 22, -2.5, 2.5]
+- [PT, 22, 20, 13000]
 
 # no transverse impulses
 BEAM_REMNANTS: false

prog/runcards/pp/makefile

@@ -12,4 +12,4 @@ ROOT_DIR:=$(shell dirname $(realpath $(firstword $(MAKEFILE_LIST))))
 analysis: Sherpa.yaml
 	RIVET_ANALYSIS_PATH="$(ROOT_DIR)/../../analysis/qqgg_proton" mpirun --use-hwthread-cpus --use-hwthread-cpus Sherpa
 	rivet-mkhtml *.yoda -o analysis
-	cp Analysis.yoda analysis
+	cp *.yoda analysis

prog/runcards/qqgg/Sherpa.yaml

@@ -1,7 +1,7 @@
 # general options
 OUTPUT: 3
 
-# no events
+# event count
 EVENTS: 1000000
 
 # save results