analyze the partonic xs, implement cut for "second" photon

prog/analysis/qqgg_partonic/.gitignore

@@ -0,0 +1 @@
+*.so

prog/analysis/qqgg_partonic/MC_DIPHOTON_PARTONIC.cc

@@ -0,0 +1,74 @@
+// -*- C++ -*-
+#include "Rivet/Analysis.hh"
+#include "Rivet/Projections/DressedLeptons.hh"
+#include "Rivet/Projections/FastJets.hh"
+#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Projections/MissingMomentum.hh"
+#include "Rivet/Projections/PromptFinalState.hh"
+
+namespace Rivet {
+
+/// @brief Generate some simple histograms of the diphoton process. TODO: more
+class MC_DIPHOTON_PARTONIC : public Analysis {
+public:
+  /// Constructor
+  DEFAULT_RIVET_ANALYSIS_CTOR(MC_DIPHOTON_PARTONIC);
+
+  /// @name Analysis methods
+  //@{
+
+  /// Book histograms and initialise projections before the run
+  void init() {
+    // Initialise and register projections
+
+    // The basic final-state projection:
+    // all final-state particles within
+    // the given eta acceptance
+    // We hardcode this TODO: make this configurable
+    FinalState fs;
+    declare(fs, "FS");
+
+    // cut has been made in sherpa
+    IdentifiedFinalState ifs {};
+    ifs.acceptId(PID::PHOTON);
+    declare(ifs, "IFS");
+
+    auto energy = info().energies()[0].first;
+    book(_h_pT, "pT", 50, 16.3, energy);
+    book(_h_eta, "eta", 50, -2.5, 2.5);
+    book(_h_cos_theta, "cos_theta", 50, -.986, .986);
+  }
+
+  /// Perform the per-event analysis
+  void analyze(const Event &event) {
+    const double weight = 1.0;
+
+    Particles photons = apply<IdentifiedFinalState>(event, "IFS").particles();
+
+    // they are both the same, so we take the first
+    for (const auto &photon : photons) {
+      _h_pT->fill(photon.pT(), weight);
+      _h_eta->fill(photon.eta(), weight);
+      _h_cos_theta->fill(cos(photon.theta()), weight);
+    }
+  }
+
+  //@}
+
+  void finalize() {
+    normalize(_h_pT);
+    normalize(_h_eta);
+    normalize(_h_cos_theta);
+  }
+
+  /// @name Histograms
+  //@{
+  Histo1DPtr _h_pT;
+  Histo1DPtr _h_eta;
+  Histo1DPtr _h_cos_theta;
+  //@}
+};
+
+DECLARE_RIVET_PLUGIN(MC_DIPHOTON_PARTONIC);
+
+} // namespace Rivet

prog/analysis/qqgg_partonic/MC_DIPHOTON_PARTONIC.info

@@ -0,0 +1,25 @@
+# -*- yaml -*-
+Name: MC_DIPHOTON_PARTONIC
+Year: 2020
+Summary: Simple analysis for a bare d d_bar anihilation into two photons.
+Status: UNVALIDATED NOTREENTRY NOHEPDATA SINGLEWEIGHT UNPHYSICAL
+Authors:
+ - Valentin Boettcher <hiro@protagon.space>
+#RunInfo: <Describe event types, cuts, and other general generator config tips.>
+Beams: [1, -1]
+Energies: [[3250,6500]]
+#Luminosity_fb: 139.0
+Description:
+  'A brief description of what is measured and what it is useful for
+  in terms of MC testing, tuning, reinterpretation, etc. Use LaTeX
+  for maths, e.g. $\pT > 50\;\GeV$.'
+ValidationInfo:
+  'A description of the process used to validate the Rivet code against
+  the original experimental analysis. Cut-flow tables and similar information
+  are welcome'
+#ReleaseTests:
+# - $A my-hepmc-prefix :MODE=some_rivet_flag
+Keywords: []
+ToDo:
+ - Implement the analysis, test it, remove this ToDo, and mark as VALIDATED :-)
+ - add more documentation in here

prog/analysis/qqgg_partonic/MC_DIPHOTON_PARTONIC.plot

@@ -0,0 +1,8 @@
+BEGIN PLOT /MC_DIPHOTON_PARTONIC/d01-x01-y01
+Title=[Insert title for histogram d01-x01-y01 here]
+XLabel=[Insert $x$-axis label for histogram d01-x01-y01 here]
+YLabel=[Insert $y$-axis label for histogram d01-x01-y01 here]
+# + any additional plot settings you might like, see make-plots documentation
+END PLOT
+
+# ... add more histograms as you need them ...

prog/analysis/qqgg_partonic/makefile

@@ -0,0 +1,2 @@
+RivetMcDiphotonSimple.so: MC_DIPHOTON_PARTONIC.cc
+	rivet-build RivetMcDiphotonSimple.so MC_DIPHOTON_PARTONIC.cc

prog/python/qqgg/lab_xs/lab_xs.c

@@ -1044,6 +1044,7 @@ static const char __pyx_k_lab_xs[] = "lab_xs";
 static const char __pyx_k_e_proton[] = "e_proton";
 static const char __pyx_k_tanh_eta[] = "tanh_eta";
 static const char __pyx_k_lab_xs_pyx[] = "lab_xs.pyx";
+static const char __pyx_k_second_eta[] = "second_eta";
 static const char __pyx_k_diff_xs_eta[] = "diff_xs_eta";
 static const char __pyx_k_cline_in_traceback[] = "cline_in_traceback";
 static const char __pyx_k_averaged_tchanel_q2[] = "averaged_tchanel_q2";
@@ -1059,6 +1060,7 @@ static PyObject *__pyx_n_s_main;
 static PyObject *__pyx_n_s_name;
 static PyObject *__pyx_n_s_pT;
 static PyObject *__pyx_n_s_rap;
+static PyObject *__pyx_n_s_second_eta;
 static PyObject *__pyx_n_s_tanh_eta;
 static PyObject *__pyx_n_s_test;
 static PyObject *__pyx_n_s_x_1;
@@ -1066,12 +1068,15 @@ static PyObject *__pyx_n_s_x_2;
 static PyObject *__pyx_pf_6lab_xs_diff_xs_eta(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_e_proton, double __pyx_v_charge, double __pyx_v_eta, double __pyx_v_x_1, double __pyx_v_x_2); /* proto */
 static PyObject *__pyx_pf_6lab_xs_2pT(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_e_proton, double __pyx_v_eta, double __pyx_v_x_1, double __pyx_v_x_2); /* proto */
 static PyObject *__pyx_pf_6lab_xs_4averaged_tchanel_q2(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_e_proton, double __pyx_v_x_1, double __pyx_v_x_2); /* proto */
+static PyObject *__pyx_pf_6lab_xs_6second_eta(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_eta, double __pyx_v_x_1, double __pyx_v_x_2); /* proto */
 static PyObject *__pyx_tuple_;
 static PyObject *__pyx_tuple__3;
 static PyObject *__pyx_tuple__5;
+static PyObject *__pyx_tuple__7;
 static PyObject *__pyx_codeobj__2;
 static PyObject *__pyx_codeobj__4;
 static PyObject *__pyx_codeobj__6;
+static PyObject *__pyx_codeobj__8;
 /* Late includes */
 
 /* "lab_xs.pyx":3
@@ -1355,7 +1360,7 @@ static PyObject *__pyx_pf_6lab_xs_2pT(CYTHON_UNUSED PyObject *__pyx_self, double
   /* "lab_xs.pyx":30
  *         * x_2
  *         / (x_1 + x_2 - (x_1 - x_2) * tanh_eta)
- *         * sqrt((1 - tanh_eta ** 2))             # <<<<<<<<<<<<<<
+ *         * sqrt(1 - tanh_eta ** 2)             # <<<<<<<<<<<<<<
  *     )
  * 
  */
@@ -1389,6 +1394,7 @@ static PyObject *__pyx_pf_6lab_xs_2pT(CYTHON_UNUSED PyObject *__pyx_self, double
  * 
  * def averaged_tchanel_q2(double e_proton, double x_1, double x_2):             # <<<<<<<<<<<<<<
  *     return 2 * x_1 * x_2 * e_proton ** 2
+ * 
  */
 
 /* Python wrapper */
@@ -1474,6 +1480,8 @@ static PyObject *__pyx_pf_6lab_xs_4averaged_tchanel_q2(CYTHON_UNUSED PyObject *_
  * 
  * def averaged_tchanel_q2(double e_proton, double x_1, double x_2):
  *     return 2 * x_1 * x_2 * e_proton ** 2             # <<<<<<<<<<<<<<
+ * 
+ * def second_eta(double eta, double x_1, double x_2):
  */
   __Pyx_XDECREF(__pyx_r);
   __pyx_t_1 = PyFloat_FromDouble((((2.0 * __pyx_v_x_1) * __pyx_v_x_2) * pow(__pyx_v_e_proton, 2.0))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 35, __pyx_L1_error)
@@ -1487,6 +1495,7 @@ static PyObject *__pyx_pf_6lab_xs_4averaged_tchanel_q2(CYTHON_UNUSED PyObject *_
  * 
  * def averaged_tchanel_q2(double e_proton, double x_1, double x_2):             # <<<<<<<<<<<<<<
  *     return 2 * x_1 * x_2 * e_proton ** 2
+ * 
  */
 
   /* function exit code */
@@ -1500,6 +1509,133 @@ static PyObject *__pyx_pf_6lab_xs_4averaged_tchanel_q2(CYTHON_UNUSED PyObject *_
   return __pyx_r;
 }
 
+/* "lab_xs.pyx":37
+ *     return 2 * x_1 * x_2 * e_proton ** 2
+ * 
+ * def second_eta(double eta, double x_1, double x_2):             # <<<<<<<<<<<<<<
+ *     cdef double rap = atanh((x_1 - x_2) / (x_1 + x_2))
+ *     return -(-eta + 2 * rap)
+ */
+
+/* Python wrapper */
+static PyObject *__pyx_pw_6lab_xs_7second_eta(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
+static PyMethodDef __pyx_mdef_6lab_xs_7second_eta = {"second_eta", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_6lab_xs_7second_eta, METH_VARARGS|METH_KEYWORDS, 0};
+static PyObject *__pyx_pw_6lab_xs_7second_eta(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
+  double __pyx_v_eta;
+  double __pyx_v_x_1;
+  double __pyx_v_x_2;
+  PyObject *__pyx_r = 0;
+  __Pyx_RefNannyDeclarations
+  __Pyx_RefNannySetupContext("second_eta (wrapper)", 0);
+  {
+    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_eta,&__pyx_n_s_x_1,&__pyx_n_s_x_2,0};
+    PyObject* values[3] = {0,0,0};
+    if (unlikely(__pyx_kwds)) {
+      Py_ssize_t kw_args;
+      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
+      switch (pos_args) {
+        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
+        CYTHON_FALLTHROUGH;
+        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
+        CYTHON_FALLTHROUGH;
+        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
+        CYTHON_FALLTHROUGH;
+        case  0: break;
+        default: goto __pyx_L5_argtuple_error;
+      }
+      kw_args = PyDict_Size(__pyx_kwds);
+      switch (pos_args) {
+        case  0:
+        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_eta)) != 0)) kw_args--;
+        else goto __pyx_L5_argtuple_error;
+        CYTHON_FALLTHROUGH;
+        case  1:
+        if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_x_1)) != 0)) kw_args--;
+        else {
+          __Pyx_RaiseArgtupleInvalid("second_eta", 1, 3, 3, 1); __PYX_ERR(0, 37, __pyx_L3_error)
+        }
+        CYTHON_FALLTHROUGH;
+        case  2:
+        if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_x_2)) != 0)) kw_args--;
+        else {
+          __Pyx_RaiseArgtupleInvalid("second_eta", 1, 3, 3, 2); __PYX_ERR(0, 37, __pyx_L3_error)
+        }
+      }
+      if (unlikely(kw_args > 0)) {
+        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "second_eta") < 0)) __PYX_ERR(0, 37, __pyx_L3_error)
+      }
+    } else if (PyTuple_GET_SIZE(__pyx_args) != 3) {
+      goto __pyx_L5_argtuple_error;
+    } else {
+      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
+      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
+      values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
+    }
+    __pyx_v_eta = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_eta == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 37, __pyx_L3_error)
+    __pyx_v_x_1 = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_x_1 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 37, __pyx_L3_error)
+    __pyx_v_x_2 = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_x_2 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 37, __pyx_L3_error)
+  }
+  goto __pyx_L4_argument_unpacking_done;
+  __pyx_L5_argtuple_error:;
+  __Pyx_RaiseArgtupleInvalid("second_eta", 1, 3, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 37, __pyx_L3_error)
+  __pyx_L3_error:;
+  __Pyx_AddTraceback("lab_xs.second_eta", __pyx_clineno, __pyx_lineno, __pyx_filename);
+  __Pyx_RefNannyFinishContext();
+  return NULL;
+  __pyx_L4_argument_unpacking_done:;
+  __pyx_r = __pyx_pf_6lab_xs_6second_eta(__pyx_self, __pyx_v_eta, __pyx_v_x_1, __pyx_v_x_2);
+
+  /* function exit code */
+  __Pyx_RefNannyFinishContext();
+  return __pyx_r;
+}
+
+static PyObject *__pyx_pf_6lab_xs_6second_eta(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_eta, double __pyx_v_x_1, double __pyx_v_x_2) {
+  double __pyx_v_rap;
+  PyObject *__pyx_r = NULL;
+  __Pyx_RefNannyDeclarations
+  PyObject *__pyx_t_1 = NULL;
+  __Pyx_RefNannySetupContext("second_eta", 0);
+
+  /* "lab_xs.pyx":38
+ * 
+ * def second_eta(double eta, double x_1, double x_2):
+ *     cdef double rap = atanh((x_1 - x_2) / (x_1 + x_2))             # <<<<<<<<<<<<<<
+ *     return -(-eta + 2 * rap)
+ */
+  __pyx_v_rap = atanh(((__pyx_v_x_1 - __pyx_v_x_2) / (__pyx_v_x_1 + __pyx_v_x_2)));
+
+  /* "lab_xs.pyx":39
+ * def second_eta(double eta, double x_1, double x_2):
+ *     cdef double rap = atanh((x_1 - x_2) / (x_1 + x_2))
+ *     return -(-eta + 2 * rap)             # <<<<<<<<<<<<<<
+ */
+  __Pyx_XDECREF(__pyx_r);
+  __pyx_t_1 = PyFloat_FromDouble((-((-__pyx_v_eta) + (2.0 * __pyx_v_rap)))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 39, __pyx_L1_error)
+  __Pyx_GOTREF(__pyx_t_1);
+  __pyx_r = __pyx_t_1;
+  __pyx_t_1 = 0;
+  goto __pyx_L0;
+
+  /* "lab_xs.pyx":37
+ *     return 2 * x_1 * x_2 * e_proton ** 2
+ * 
+ * def second_eta(double eta, double x_1, double x_2):             # <<<<<<<<<<<<<<
+ *     cdef double rap = atanh((x_1 - x_2) / (x_1 + x_2))
+ *     return -(-eta + 2 * rap)
+ */
+
+  /* function exit code */
+  __pyx_L1_error:;
+  __Pyx_XDECREF(__pyx_t_1);
+  __Pyx_AddTraceback("lab_xs.second_eta", __pyx_clineno, __pyx_lineno, __pyx_filename);
+  __pyx_r = NULL;
+  __pyx_L0:;
+  __Pyx_XGIVEREF(__pyx_r);
+  __Pyx_RefNannyFinishContext();
+  return __pyx_r;
+}
+
 static PyMethodDef __pyx_methods[] = {
   {0, 0, 0, 0}
 };
@@ -1558,6 +1694,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
   {&__pyx_n_s_name, __pyx_k_name, sizeof(__pyx_k_name), 0, 0, 1, 1},
   {&__pyx_n_s_pT, __pyx_k_pT, sizeof(__pyx_k_pT), 0, 0, 1, 1},
   {&__pyx_n_s_rap, __pyx_k_rap, sizeof(__pyx_k_rap), 0, 0, 1, 1},
+  {&__pyx_n_s_second_eta, __pyx_k_second_eta, sizeof(__pyx_k_second_eta), 0, 0, 1, 1},
   {&__pyx_n_s_tanh_eta, __pyx_k_tanh_eta, sizeof(__pyx_k_tanh_eta), 0, 0, 1, 1},
   {&__pyx_n_s_test, __pyx_k_test, sizeof(__pyx_k_test), 0, 0, 1, 1},
   {&__pyx_n_s_x_1, __pyx_k_x_1, sizeof(__pyx_k_x_1), 0, 0, 1, 1},
@@ -1601,11 +1738,24 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) {
  * 
  * def averaged_tchanel_q2(double e_proton, double x_1, double x_2):             # <<<<<<<<<<<<<<
  *     return 2 * x_1 * x_2 * e_proton ** 2
+ * 
  */
   __pyx_tuple__5 = PyTuple_Pack(3, __pyx_n_s_e_proton, __pyx_n_s_x_1, __pyx_n_s_x_2); if (unlikely(!__pyx_tuple__5)) __PYX_ERR(0, 34, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_tuple__5);
   __Pyx_GIVEREF(__pyx_tuple__5);
   __pyx_codeobj__6 = (PyObject*)__Pyx_PyCode_New(3, 0, 3, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__5, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_lab_xs_pyx, __pyx_n_s_averaged_tchanel_q2, 34, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__6)) __PYX_ERR(0, 34, __pyx_L1_error)
+
+  /* "lab_xs.pyx":37
+ *     return 2 * x_1 * x_2 * e_proton ** 2
+ * 
+ * def second_eta(double eta, double x_1, double x_2):             # <<<<<<<<<<<<<<
+ *     cdef double rap = atanh((x_1 - x_2) / (x_1 + x_2))
+ *     return -(-eta + 2 * rap)
+ */
+  __pyx_tuple__7 = PyTuple_Pack(4, __pyx_n_s_eta, __pyx_n_s_x_1, __pyx_n_s_x_2, __pyx_n_s_rap); if (unlikely(!__pyx_tuple__7)) __PYX_ERR(0, 37, __pyx_L1_error)
+  __Pyx_GOTREF(__pyx_tuple__7);
+  __Pyx_GIVEREF(__pyx_tuple__7);
+  __pyx_codeobj__8 = (PyObject*)__Pyx_PyCode_New(3, 0, 4, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__7, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_lab_xs_pyx, __pyx_n_s_second_eta, 37, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__8)) __PYX_ERR(0, 37, __pyx_L1_error)
   __Pyx_RefNannyFinishContext();
   return 0;
   __pyx_L1_error:;
@@ -1910,12 +2060,25 @@ if (!__Pyx_RefNanny) {
  * 
  * def averaged_tchanel_q2(double e_proton, double x_1, double x_2):             # <<<<<<<<<<<<<<
  *     return 2 * x_1 * x_2 * e_proton ** 2
+ * 
  */
   __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6lab_xs_5averaged_tchanel_q2, NULL, __pyx_n_s_lab_xs); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 34, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   if (PyDict_SetItem(__pyx_d, __pyx_n_s_averaged_tchanel_q2, __pyx_t_1) < 0) __PYX_ERR(0, 34, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 
+  /* "lab_xs.pyx":37
+ *     return 2 * x_1 * x_2 * e_proton ** 2
+ * 
+ * def second_eta(double eta, double x_1, double x_2):             # <<<<<<<<<<<<<<
+ *     cdef double rap = atanh((x_1 - x_2) / (x_1 + x_2))
+ *     return -(-eta + 2 * rap)
+ */
+  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6lab_xs_7second_eta, NULL, __pyx_n_s_lab_xs); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 37, __pyx_L1_error)
+  __Pyx_GOTREF(__pyx_t_1);
+  if (PyDict_SetItem(__pyx_d, __pyx_n_s_second_eta, __pyx_t_1) < 0) __PYX_ERR(0, 37, __pyx_L1_error)
+  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+
   /* "lab_xs.pyx":1
  * from libc.math cimport tanh, atanh, sqrt             # <<<<<<<<<<<<<<
  * 

prog/python/qqgg/lab_xs/lab_xs.html

@@ -634,7 +634,7 @@ static PyObject *__pyx_pf_6lab_xs_2pT(CYTHON_UNUSED PyObject *__pyx_self, double
 <pre class="cython line score-0">&#xA0;<span class="">27</span>:         <span class="o">*</span> <span class="n">x_1</span></pre>
 <pre class="cython line score-0">&#xA0;<span class="">28</span>:         <span class="o">*</span> <span class="n">x_2</span></pre>
 <pre class="cython line score-0">&#xA0;<span class="">29</span>:         <span class="o">/</span> <span class="p">(</span><span class="n">x_1</span> <span class="o">+</span> <span class="n">x_2</span> <span class="o">-</span> <span class="p">(</span><span class="n">x_1</span> <span class="o">-</span> <span class="n">x_2</span><span class="p">)</span> <span class="o">*</span> <span class="n">tanh_eta</span><span class="p">)</span></pre>
-<pre class="cython line score-5" onclick="(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)">+<span class="">30</span>:         <span class="o">*</span> <span class="n">sqrt</span><span class="p">((</span><span class="mf">1</span> <span class="o">-</span> <span class="n">tanh_eta</span> <span class="o">**</span> <span class="mf">2</span><span class="p">))</span></pre>
+<pre class="cython line score-5" onclick="(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)">+<span class="">30</span>:         <span class="o">*</span> <span class="n">sqrt</span><span class="p">(</span><span class="mf">1</span> <span class="o">-</span> <span class="n">tanh_eta</span> <span class="o">**</span> <span class="mf">2</span><span class="p">)</span></pre>
 <pre class='cython code score-5 '>  __pyx_t_1 = <span class='py_c_api'>PyFloat_FromDouble</span>((((((2.0 * __pyx_v_e_proton) * __pyx_v_x_1) * __pyx_v_x_2) / ((__pyx_v_x_1 + __pyx_v_x_2) - ((__pyx_v_x_1 - __pyx_v_x_2) * __pyx_v_tanh_eta))) * sqrt((1.0 - pow(__pyx_v_tanh_eta, 2.0)))));<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 30, __pyx_L1_error)</span>
   <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);
   __pyx_r = __pyx_t_1;
@@ -643,8 +643,8 @@ static PyObject *__pyx_pf_6lab_xs_2pT(CYTHON_UNUSED PyObject *__pyx_self, double
 </pre><pre class="cython line score-0">&#xA0;<span class="">31</span>:     <span class="p">)</span></pre>
 <pre class="cython line score-0">&#xA0;<span class="">32</span>: </pre>
 <pre class="cython line score-0">&#xA0;<span class="">33</span>: </pre>
-<pre class="cython line score-74" onclick="(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)">+<span class="">34</span>: <span class="k">def</span> <span class="nf">averaged_tchanel_q2</span><span class="p">(</span><span class="n">double</span> <span class="n">e_proton</span><span class="p">,</span> <span class="n">double</span> <span class="n">x_1</span><span class="p">,</span> <span class="n">double</span> <span class="n">x_2</span><span class="p">):</span></pre>
-<pre class='cython code score-74 '>/* Python wrapper */
+<pre class="cython line score-76" onclick="(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)">+<span class="">34</span>: <span class="k">def</span> <span class="nf">averaged_tchanel_q2</span><span class="p">(</span><span class="n">double</span> <span class="n">e_proton</span><span class="p">,</span> <span class="n">double</span> <span class="n">x_1</span><span class="p">,</span> <span class="n">double</span> <span class="n">x_2</span><span class="p">):</span></pre>
+<pre class='cython code score-76 '>/* Python wrapper */
 static PyObject *__pyx_pw_6lab_xs_5averaged_tchanel_q2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
 static PyMethodDef __pyx_mdef_6lab_xs_5averaged_tchanel_q2 = {"averaged_tchanel_q2", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_6lab_xs_5averaged_tchanel_q2, METH_VARARGS|METH_KEYWORDS, 0};
 static PyObject *__pyx_pw_6lab_xs_5averaged_tchanel_q2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
@@ -741,6 +741,7 @@ static PyObject *__pyx_pf_6lab_xs_4averaged_tchanel_q2(CYTHON_UNUSED PyObject *_
   <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);
   if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_averaged_tchanel_q2, __pyx_t_1) &lt; 0) <span class='error_goto'>__PYX_ERR(0, 34, __pyx_L1_error)</span>
   <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;
+  __pyx_codeobj__6 = (PyObject*)<span class='pyx_c_api'>__Pyx_PyCode_New</span>(3, 0, 3, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__5, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_lab_xs_pyx, __pyx_n_s_averaged_tchanel_q2, 34, __pyx_empty_bytes);<span class='error_goto'> if (unlikely(!__pyx_codeobj__6)) __PYX_ERR(0, 34, __pyx_L1_error)</span>
 </pre><pre class="cython line score-6" onclick="(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)">+<span class="">35</span>:     <span class="k">return</span> <span class="mf">2</span> <span class="o">*</span> <span class="n">x_1</span> <span class="o">*</span> <span class="n">x_2</span> <span class="o">*</span> <span class="n">e_proton</span> <span class="o">**</span> <span class="mf">2</span></pre>
 <pre class='cython code score-6 '>  <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_r);
   __pyx_t_1 = <span class='py_c_api'>PyFloat_FromDouble</span>((((2.0 * __pyx_v_x_1) * __pyx_v_x_2) * pow(__pyx_v_e_proton, 2.0)));<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 35, __pyx_L1_error)</span>
@@ -748,4 +749,113 @@ static PyObject *__pyx_pf_6lab_xs_4averaged_tchanel_q2(CYTHON_UNUSED PyObject *_
   __pyx_r = __pyx_t_1;
   __pyx_t_1 = 0;
   goto __pyx_L0;
+</pre><pre class="cython line score-0">&#xA0;<span class="">36</span>: </pre>
+<pre class="cython line score-74" onclick="(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)">+<span class="">37</span>: <span class="k">def</span> <span class="nf">second_eta</span><span class="p">(</span><span class="n">double</span> <span class="n">eta</span><span class="p">,</span> <span class="n">double</span> <span class="n">x_1</span><span class="p">,</span> <span class="n">double</span> <span class="n">x_2</span><span class="p">):</span></pre>
+<pre class='cython code score-74 '>/* Python wrapper */
+static PyObject *__pyx_pw_6lab_xs_7second_eta(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
+static PyMethodDef __pyx_mdef_6lab_xs_7second_eta = {"second_eta", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_6lab_xs_7second_eta, METH_VARARGS|METH_KEYWORDS, 0};
+static PyObject *__pyx_pw_6lab_xs_7second_eta(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
+  double __pyx_v_eta;
+  double __pyx_v_x_1;
+  double __pyx_v_x_2;
+  PyObject *__pyx_r = 0;
+  <span class='refnanny'>__Pyx_RefNannyDeclarations</span>
+  <span class='refnanny'>__Pyx_RefNannySetupContext</span>("second_eta (wrapper)", 0);
+  {
+    static PyObject **__pyx_pyargnames[] = {&amp;__pyx_n_s_eta,&amp;__pyx_n_s_x_1,&amp;__pyx_n_s_x_2,0};
+    PyObject* values[3] = {0,0,0};
+    if (unlikely(__pyx_kwds)) {
+      Py_ssize_t kw_args;
+      const Py_ssize_t pos_args = <span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args);
+      switch (pos_args) {
+        case  3: values[2] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 2);
+        CYTHON_FALLTHROUGH;
+        case  2: values[1] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 1);
+        CYTHON_FALLTHROUGH;
+        case  1: values[0] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 0);
+        CYTHON_FALLTHROUGH;
+        case  0: break;
+        default: goto __pyx_L5_argtuple_error;
+      }
+      kw_args = <span class='py_c_api'>PyDict_Size</span>(__pyx_kwds);
+      switch (pos_args) {
+        case  0:
+        if (likely((values[0] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_eta)) != 0)) kw_args--;
+        else goto __pyx_L5_argtuple_error;
+        CYTHON_FALLTHROUGH;
+        case  1:
+        if (likely((values[1] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_x_1)) != 0)) kw_args--;
+        else {
+          <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>("second_eta", 1, 3, 3, 1); <span class='error_goto'>__PYX_ERR(0, 37, __pyx_L3_error)</span>
+        }
+        CYTHON_FALLTHROUGH;
+        case  2:
+        if (likely((values[2] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_x_2)) != 0)) kw_args--;
+        else {
+          <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>("second_eta", 1, 3, 3, 2); <span class='error_goto'>__PYX_ERR(0, 37, __pyx_L3_error)</span>
+        }
+      }
+      if (unlikely(kw_args &gt; 0)) {
+        if (unlikely(<span class='pyx_c_api'>__Pyx_ParseOptionalKeywords</span>(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "second_eta") &lt; 0)) <span class='error_goto'>__PYX_ERR(0, 37, __pyx_L3_error)</span>
+      }
+    } else if (<span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args) != 3) {
+      goto __pyx_L5_argtuple_error;
+    } else {
+      values[0] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 0);
+      values[1] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 1);
+      values[2] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 2);
+    }
+    __pyx_v_eta = __pyx_<span class='py_c_api'>PyFloat_AsDouble</span>(values[0]); if (unlikely((__pyx_v_eta == (double)-1) &amp;&amp; <span class='py_c_api'>PyErr_Occurred</span>())) <span class='error_goto'>__PYX_ERR(0, 37, __pyx_L3_error)</span>
+    __pyx_v_x_1 = __pyx_<span class='py_c_api'>PyFloat_AsDouble</span>(values[1]); if (unlikely((__pyx_v_x_1 == (double)-1) &amp;&amp; <span class='py_c_api'>PyErr_Occurred</span>())) <span class='error_goto'>__PYX_ERR(0, 37, __pyx_L3_error)</span>
+    __pyx_v_x_2 = __pyx_<span class='py_c_api'>PyFloat_AsDouble</span>(values[2]); if (unlikely((__pyx_v_x_2 == (double)-1) &amp;&amp; <span class='py_c_api'>PyErr_Occurred</span>())) <span class='error_goto'>__PYX_ERR(0, 37, __pyx_L3_error)</span>
+  }
+  goto __pyx_L4_argument_unpacking_done;
+  __pyx_L5_argtuple_error:;
+  <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>("second_eta", 1, 3, 3, <span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args)); <span class='error_goto'>__PYX_ERR(0, 37, __pyx_L3_error)</span>
+  __pyx_L3_error:;
+  <span class='pyx_c_api'>__Pyx_AddTraceback</span>("lab_xs.second_eta", __pyx_clineno, __pyx_lineno, __pyx_filename);
+  <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();
+  return NULL;
+  __pyx_L4_argument_unpacking_done:;
+  __pyx_r = __pyx_pf_6lab_xs_6second_eta(__pyx_self, __pyx_v_eta, __pyx_v_x_1, __pyx_v_x_2);
+
+  /* function exit code */
+  <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();
+  return __pyx_r;
+}
+
+static PyObject *__pyx_pf_6lab_xs_6second_eta(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_eta, double __pyx_v_x_1, double __pyx_v_x_2) {
+  double __pyx_v_rap;
+  PyObject *__pyx_r = NULL;
+  <span class='refnanny'>__Pyx_RefNannyDeclarations</span>
+  <span class='refnanny'>__Pyx_RefNannySetupContext</span>("second_eta", 0);
+/* … */
+  /* function exit code */
+  __pyx_L1_error:;
+  <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_t_1);
+  <span class='pyx_c_api'>__Pyx_AddTraceback</span>("lab_xs.second_eta", __pyx_clineno, __pyx_lineno, __pyx_filename);
+  __pyx_r = NULL;
+  __pyx_L0:;
+  <span class='refnanny'>__Pyx_XGIVEREF</span>(__pyx_r);
+  <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();
+  return __pyx_r;
+}
+/* … */
+  __pyx_tuple__7 = <span class='py_c_api'>PyTuple_Pack</span>(4, __pyx_n_s_eta, __pyx_n_s_x_1, __pyx_n_s_x_2, __pyx_n_s_rap);<span class='error_goto'> if (unlikely(!__pyx_tuple__7)) __PYX_ERR(0, 37, __pyx_L1_error)</span>
+  <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_tuple__7);
+  <span class='refnanny'>__Pyx_GIVEREF</span>(__pyx_tuple__7);
+/* … */
+  __pyx_t_1 = PyCFunction_NewEx(&amp;__pyx_mdef_6lab_xs_7second_eta, NULL, __pyx_n_s_lab_xs);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 37, __pyx_L1_error)</span>
+  <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);
+  if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_second_eta, __pyx_t_1) &lt; 0) <span class='error_goto'>__PYX_ERR(0, 37, __pyx_L1_error)</span>
+  <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;
+</pre><pre class="cython line score-0" onclick="(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)">+<span class="">38</span>:     <span class="k">cdef</span> <span class="kt">double</span> <span class="nf">rap</span> <span class="o">=</span> <span class="n">atanh</span><span class="p">((</span><span class="n">x_1</span> <span class="o">-</span> <span class="n">x_2</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">x_1</span> <span class="o">+</span> <span class="n">x_2</span><span class="p">))</span></pre>
+<pre class='cython code score-0 '>  __pyx_v_rap = atanh(((__pyx_v_x_1 - __pyx_v_x_2) / (__pyx_v_x_1 + __pyx_v_x_2)));
+</pre><pre class="cython line score-6" onclick="(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)">+<span class="">39</span>:     <span class="k">return</span> <span class="o">-</span><span class="p">(</span><span class="o">-</span><span class="n">eta</span> <span class="o">+</span> <span class="mf">2</span> <span class="o">*</span> <span class="n">rap</span><span class="p">)</span></pre>
+<pre class='cython code score-6 '>  <span class='pyx_macro_api'>__Pyx_XDECREF</span>(__pyx_r);
+  __pyx_t_1 = <span class='py_c_api'>PyFloat_FromDouble</span>((-((-__pyx_v_eta) + (2.0 * __pyx_v_rap))));<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 39, __pyx_L1_error)</span>
+  <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);
+  __pyx_r = __pyx_t_1;
+  __pyx_t_1 = 0;
+  goto __pyx_L0;
 </pre></div></body></html>

prog/python/qqgg/lab_xs/lab_xs.pyx

@@ -34,3 +34,6 @@ def pT(double e_proton, double eta, double x_1, double x_2):
 def averaged_tchanel_q2(double e_proton, double x_1, double x_2):
     return 2 * x_1 * x_2 * e_proton ** 2
 
+def second_eta(double eta, double x_1, double x_2):
+    cdef double rap = atanh((x_1 - x_2) / (x_1 + x_2))
+    return -(-eta + 2 * rap)

prog/python/qqgg/parton_density_function_stuff.org

@@ -27,12 +27,13 @@
 
 ** Global Config
 #+begin_src jupyter-python :exports both :results raw drawer
-η = 1
-e_proton = 6500  # GeV
-interval_η = [-η, η]
-interval = η_to_θ([-η, η])
-interval_cosθ = np.cos(interval)
-num_samples = 10_000
+  η = 1
+  min_pT = 200
+  e_proton = 6500  # GeV
+  interval_η = [-η, η]
+  interval = η_to_θ([-η, η])
+  interval_cosθ = np.cos(interval)
+  num_samples = 10_000
 #+end_src
 
 #+RESULTS:
@@ -53,6 +54,7 @@ We begin by implementing the same sermon for the lab frame.
   from numba import jit, vectorize, float64, boolean
   import lab_xs.lab_xs as c_xs
 
+
   @vectorize([float64(float64, float64, float64, float64)], nopython=True)
   def energy_factor(e_proton, charge, x_1, x_2):
       """Calculates the factor common to all other values in this module.
@@ -140,15 +142,84 @@ We begin by implementing the same sermon for the lab frame.
       return f * ((np.tanh(η - rap)) ** 2 + 1)
 
 
+  class Cut:
+      def __init__(self):
+          self._other = None
+          self._current_comb = self._call
+
+          self._greater_than = 0
+          self._lower_than = np.inf
+
+      def __gt__(self, greater_than):
+          self._greater_than = greater_than
+
+          return self
+
+      def __lt__(self, lower_than):
+          self._lower_than = lower_than
+
+          return self
+
+      def _or_comb(self, event):
+          return self._call(event) or self._other(event)
+
+      def _and_comb(self, event):
+          return self._call(event) and self._other(event)
+
+      def _call(self, event):
+          return self._greater_than < self._calculate(event) < self._lower_than
+
+      def _calculate(self, event):
+          raise NotImplementedError('"_calulate" must be implemented.')
+
+      def __call__(self, event):
+          return self._current_comb(event)
+
+      def __and__(self, other):
+          self._other = other
+          self._current_comb = self._and_comb
+
+          return self
+
+      def __or__(self, other):
+          self._other = other
+          self._current_comb = self._or_comb
+
+          return self
+
+      def apply(self, function):
+          @wraps(function)
+          def wrapper(event):
+              if self(event):
+                  return function(event)
+
+              return 0
+
+          return wrapper
+
+
   @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_η(greater_than=0):
-      def cut(e_proton, η, x_1, x_2):
-          return c_xs.pT(e_proton, η, x_1, x_2) > greater_than
 
-      return cut
+  class CutpT(Cut):
+      def __init__(self):
+          super().__init__()
+
+      def _calculate(self, event):
+          e_hadron, eta, x_1, x_2 = event
+          return c_xs.pT(e_hadron, eta, x_1, x_2)
+
+
+  class CutOtherEta(Cut):
+      def __init__(self):
+          super().__init__()
+
+      def _calculate(self, event):
+          _, η, x_1, x_2 = event
+          rap = np.arctanh((x_1 - x_2) / (x_1 + x_2))
+          return c_xs.second_eta(η, x_1, x_2)
 #+end_src
 
 #+RESULTS:
@@ -231,7 +302,7 @@ We begin by implementing the same sermon for the lab frame.
       xfxQ2 = pdf.xfxQ2
 
       def distribution(event: np.ndarray) -> float:
-          if cut and not cut(e_hadron, *event):
+          if cut and not cut([e_hadron, *event]):
               return 0
 
           angle_arg, x_1, x_2 = event
@@ -289,7 +360,7 @@ We set up a new distribution.
       c_xs.diff_xs_eta,
       c_xs.averaged_tchanel_q2,
       e_proton,
-      cut=cut_pT_from_η(greater_than=200),
+      cut=(CutpT() > min_pT) & (interval_η[0] < CutOtherEta() < interval_η[1]),
   )
 #+end_src
 
@@ -304,13 +375,13 @@ Plotting it, we can see that the variance is reduced.
   pts = np.linspace(*interval_η, 1000)
 
   ax.plot(pts, [dist_η(np.array([η, 0.04, 0.04])) for η in pts])
-  ax2.plot(pts, [dist_η(np.array([η, 1, .1])) for η in pts])
+  ax2.plot(pts, [dist_η(np.array([η, 1, .5])) for η in pts])
 #+end_src
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7f4a197dd820> |
-[[file:./.ob-jupyter/038d0c37fa9e7b737cccb79451c0b89c0a6ea14c.png]]
+| <matplotlib.lines.Line2D | at | 0x7f5b8a2f5730> |
+[[file:./.ob-jupyter/a789e914a8b6fcb110b211ea596d41c763ea15e7.png]]
 :END:
 
 Lets plot how the pdf looks.
@@ -324,7 +395,7 @@ Lets plot how the pdf looks.
 
 #+RESULTS:
 :RESULTS:
-| <matplotlib.lines.Line2D | at | 0x7f4a195ddbb0> |
+| <matplotlib.lines.Line2D | at | 0x7f5b98e0e700> |
 [[file:./.ob-jupyter/b92f0c4b2c9f2195ae14444748fcdb7708d81c19.png]]
 : LHAPDF 6.2.3 loading /usr/share/lhapdf/LHAPDF/NNPDF31_lo_as_0118/NNPDF31_lo_as_0118_0000.dat
 : NNPDF31_lo_as_0118 PDF set, member #0, version 1; LHAPDF ID = 315000
@@ -395,8 +466,8 @@ gev_to_pb(eff * (intervals_η[:, 1] - intervals_η[:, 0]).prod() * 5.5e-10) * 2*
 
 #+RESULTS:
 :RESULTS:
-| <Figure | size | 432x288 | with | 2 | Axes> | (<matplotlib.axes._subplots.AxesSubplot at 0x7f4a19670e80> <matplotlib.axes._subplots.AxesSubplot at 0x7f4a193e42e0>) |
-[[file:./.ob-jupyter/7cb88ad6d47c5257514683853b4d63dbe3d9c349.png]]
+| <Figure | size | 432x288 | with | 2 | Axes> | (<matplotlib.axes._subplots.AxesSubplot at 0x7f5b975d3610> <matplotlib.axes._subplots.AxesSubplot at 0x7f5b975791c0>) |
+[[file:./.ob-jupyter/f29962e850ee73c0d1352def5de55d738f7748f1.png]]
 :END:
 
 That looks OK.
@@ -426,7 +497,7 @@ scipy.integrate.quad(lambda x: gev_to_pb(dist_η(x)), intervals_η[:, 0], interv
       c_xs.diff_xs_eta,
       c_xs.averaged_tchanel_q2,
       e_proton,
-      cut=cut_pT_from_η(greater_than=200),
+      cut=(CutpT() > min_pT) & (-1 < CutOtherEta() < 1),
       vectorize=True,
       quarks=np.array([[1, -1/3]])
   )
@@ -434,18 +505,17 @@ scipy.integrate.quad(lambda x: gev_to_pb(dist_η(x)), intervals_η[:, 0], interv
   xs_int_res = monte_carlo.integrate(
       lambda x: gev_to_pb(dist_η_vec(x)),
       np.array([[-1, 1], [pdf.xMin, 1], [pdf.xMin, 1]]),
-      num_points=800000,
+      num_points=8000000,
       adapt=False,
       epsilon=0.01,
   )
-  xs_int_res.result*2*np.pi, xs_int_res.sigma*2*np.pi
-  xs_int_res
+  xs_int_res.result*np.pi, xs_int_res.sigma*np.pi
 #+end_src
 
 #+RESULTS:
 :RESULTS:
-: 0.20826397472383126 [[-0.4523158   0.0201646   0.04736293]]
-: IntegrationResult(result=0.00042109900435229747, sigma=7.494522882354145e-06, N=800000)
+: 0.2113311302255604 [[0.20834856 0.03665743 0.02590634]]
+| 0.0007706235831937771 | 6.407338225905968e-06 |
 : LHAPDF 6.2.3 loading /usr/share/lhapdf/LHAPDF/NNPDF31_lo_as_0118/NNPDF31_lo_as_0118_0000.dat
 : NNPDF31_lo_as_0118 PDF set, member #0, version 1; LHAPDF ID = 315000
 :END:
@@ -460,7 +530,7 @@ scipy.integrate.quad(lambda x: gev_to_pb(dist_η(x)), intervals_η[:, 0], interv
 #+end_src
 
 #+RESULTS:
-: IntegrationResult(result=0.00042109900435229747, sigma=7.494522882354145e-06, N=800000)
+: IntegrationResult(result=0.0004284811348615043, sigma=7.608520813519879e-06, N=800000)
 
 
 #+begin_src jupyter-python :exports both :results raw drawer
@@ -468,7 +538,7 @@ scipy.integrate.quad(lambda x: gev_to_pb(dist_η(x)), intervals_η[:, 0], interv
 #+end_src
 
 #+RESULTS:
-: 3.9 µs ± 75.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
+: 3.95 µs ± 9.81 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
 
 #+begin_src jupyter-python :exports both :results raw drawer
 dist_η([-0.54149063, 0.01765461, 0.05391543])
@@ -637,8 +707,8 @@ I checked back with the d, dbar process. The ME there has a weird factor ~1/20**
 
 #+RESULTS:
 :RESULTS:
-: 0.9997155172485124 [[0.99999095 0.99972456]]
-: IntegrationResult(result=0.2499698460592311, sigma=9.536202506031499e-05, N=5346792)
+: 0.9996193836913989 [[0.9999086  0.99971076]]
+: IntegrationResult(result=0.24991322928674975, sigma=8.781624783123888e-05, N=6300453)
 :END:
 
 
@@ -653,3 +723,69 @@ Let's see how the pts are distributed:
 : 200.00006712310318
 
 Looks ok.
+
+**** Check out the partonic xs.
+Let's set up a cut for the η of the other photon.
+#+begin_src jupyter-python :exports both :results raw drawer
+  other_eta_cut = -2.5 < CutOtherEta() < 2.5
+#+end_src
+
+#+RESULTS:
+
+#+begin_src jupyter-python :exports both :results raw drawer
+  def part_dist(eta):
+      if isinstance(eta, np.ndarray):
+          return np.array([part_dist(s_η) for s_η in eta])
+
+      if not other_eta_cut([0, eta, .5, 1]):
+          return 0
+      return diff_xs_η(e_proton, -1 / 3, eta, 0.5, 1)
+
+  part_samples = monte_carlo.sample_unweighted_array(
+      100000,
+      part_dist,
+      interval=[-2.5, 2.5],
+      proc="auto",
+  )
+  part_samples.min()
+#+end_src
+
+#+RESULTS:
+: -2.499977473915307
+
+#+begin_src jupyter-python :exports both :results raw drawer
+part_hist = np.histogram(part_samples, bins=50, range=[-2.5, 2.5])
+fig, ax = set_up_plot()
+draw_histogram(ax, part_hist)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+: <matplotlib.axes._subplots.AxesSubplot at 0x7f5b9721d070>
+[[file:./.ob-jupyter/ab7ecd9b677441ca163faf6e2c51e0eebaa90e2f.png]]
+:END:
+
+#+begin_src jupyter-python :exports both :results raw drawer
+  yoda_sherpa_part = yoda.read("../../runcards/pp_partonic/analysis/Analysis.yoda")
+  sherpa_part_hist = yoda_to_numpy(yoda_sherpa_part["/MC_DIPHOTON_PARTONIC/eta"])
+  draw_ratio_plot(
+      [
+          dict(hist=sherpa_part_hist),
+          dict(hist=part_hist),
+      ]
+  )
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+: /home/hiro/Documents/Projects/UNI/Bachelor/prog/python/qqgg/tangled/plot_utils.py:108: RuntimeWarning: invalid value encountered in true_divide
+:   [heights / reference, edges],
+: /home/hiro/Documents/Projects/UNI/Bachelor/prog/python/qqgg/tangled/plot_utils.py:109: RuntimeWarning: invalid value encountered in true_divide
+:   errorbars=errors / reference,
+: /home/hiro/Documents/Projects/UNI/Bachelor/prog/python/qqgg/tangled/plot_utils.py:108: RuntimeWarning: invalid value encountered in true_divide
+:   [heights / reference, edges],
+: /home/hiro/Documents/Projects/UNI/Bachelor/prog/python/qqgg/tangled/plot_utils.py:109: RuntimeWarning: invalid value encountered in true_divide
+:   errorbars=errors / reference,
+| <Figure | size | 432x288 | with | 2 | Axes> | (<matplotlib.axes._subplots.AxesSubplot at 0x7f5b9542bb20> <matplotlib.axes._subplots.AxesSubplot at 0x7f5b953cd160>) |
+[[file:./.ob-jupyter/a28edcc7a43fa591d70bc343d429de6a73e8107b.png]]
+:END:

prog/python/qqgg/tangled/pdf.py

@@ -11,6 +11,7 @@ import lhapdf
 from numba import jit, vectorize, float64, boolean
 import lab_xs.lab_xs as c_xs
 
+
 @vectorize([float64(float64, float64, float64, float64)], nopython=True)
 def energy_factor(e_proton, charge, x_1, x_2):
     """Calculates the factor common to all other values in this module.
@@ -98,15 +99,84 @@ def diff_xs_η(e_proton, charge, η, x_1, x_2):
     return f * ((np.tanh(η - rap)) ** 2 + 1)
 
 
+class Cut:
+    def __init__(self):
+        self._other = None
+        self._current_comb = self._call
+
+        self._greater_than = 0
+        self._lower_than = np.inf
+
+    def __gt__(self, greater_than):
+        self._greater_than = greater_than
+
+        return self
+
+    def __lt__(self, lower_than):
+        self._lower_than = lower_than
+
+        return self
+
+    def _or_comb(self, event):
+        return self._call(event) or self._other(event)
+
+    def _and_comb(self, event):
+        return self._call(event) and self._other(event)
+
+    def _call(self, event):
+        return self._greater_than < self._calculate(event) < self._lower_than
+
+    def _calculate(self, event):
+        raise NotImplementedError('"_calulate" must be implemented.')
+
+    def __call__(self, event):
+        return self._current_comb(event)
+
+    def __and__(self, other):
+        self._other = other
+        self._current_comb = self._and_comb
+
+        return self
+
+    def __or__(self, other):
+        self._other = other
+        self._current_comb = self._or_comb
+
+        return self
+
+    def apply(self, function):
+        @wraps(function)
+        def wrapper(event):
+            if self(event):
+                return function(event)
+
+            return 0
+
+        return wrapper
+
+
 @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_η(greater_than=0):
-    def cut(e_proton, η, x_1, x_2):
-        return c_xs.pT(e_proton, η, x_1, x_2) > greater_than
 
-    return cut
+class CutpT(Cut):
+    def __init__(self):
+        super().__init__()
+
+    def _calculate(self, event):
+        e_hadron, eta, x_1, x_2 = event
+        return c_xs.pT(e_hadron, eta, x_1, x_2)
+
+
+class CutOtherEta(Cut):
+    def __init__(self):
+        super().__init__()
+
+    def _calculate(self, event):
+        _, η, x_1, x_2 = event
+        rap = np.arctanh((x_1 - x_2) / (x_1 + x_2))
+        return c_xs.second_eta(η, x_1, x_2)
 
 def cached_pdf(pdf, q, points, e_hadron):
     x_min = pdf.xMin
@@ -184,7 +254,7 @@ def get_xs_distribution_with_pdf(
     xfxQ2 = pdf.xfxQ2
 
     def distribution(event: np.ndarray) -> float:
-        if cut and not cut(e_hadron, *event):
+        if cut and not cut([e_hadron, *event]):
             return 0
 
         angle_arg, x_1, x_2 = event

prog/runcards/pp/Sherpa.yaml

@@ -2,7 +2,7 @@
 OUTPUT: 3
 
 # event count
-EVENTS: 1000000/8  # / num_cpus
+EVENTS: 10000/8  # / num_cpus
 
 # save results
 GENERATE_RESULT_DIRECTORY: true

prog/runcards/pp_partonic/Sherpa.yaml

@@ -9,7 +9,7 @@ GENERATE_RESULT_DIRECTORY: true
 
 # set up $p$ $\bar{p}$ beams, each at 100 GeV
 BEAMS: [1, -1]
-BEAM_ENERGIES: [6500, -
+BEAM_ENERGIES: [6500/2, 6500]
 
 # matrix-element calculation
 ME_GENERATORS:
@@ -34,8 +34,8 @@ MI_HANDLER: None
 
 # cut to $\abs{\eta} \leq 2.5$, and $p_T \geq 20$ GeV
 SELECTORS:
-- [Eta, 22, -1, 1]
-- [PT, 22, 200, 200000000]
+- [Eta, 22, -2.5, 2.5]
+# - [PT, 22, 200, 200000000]
 
 # no transverse impulses
 BEAM_REMNANTS: false
@@ -44,4 +44,4 @@ ANALYSIS: Rivet
 RIVET:
   IGNOREBEAMS: 1
   ANALYSES:
-    - MC_DIPHOTON_PROTON
+    - MC_DIPHOTON_PARTONIC

prog/runcards/pp_partonic/makefile

@@ -3,7 +3,7 @@ sherpa_xs: Sherpa.yaml
 
 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 "OUTPUT: 0"
+	RIVET_ANALYSIS_PATH="$(ROOT_DIR)/../../analysis/qqgg_partonic" mpirun --use-hwthread-cpus --use-hwthread-cpus Sherpa "OUTPUT: 0"
 	yodamerge *.yoda -o Analysis_all.yoda
 	rivet-mkhtml Analysis_all.yoda -o analysis
 	mv Analysis_all.yoda analysis/Analysis.yoda