diff --git a/prog/python/qqgg/lab_xs/lab_xs.html b/prog/python/qqgg/lab_xs/lab_xs.html deleted file mode 100644 index a157518..0000000 --- a/prog/python/qqgg/lab_xs/lab_xs.html +++ /dev/null @@ -1,861 +0,0 @@ - - - - - - Cython: lab_xs.pyx - - - -

Generated by Cython 0.29.17

-

- Yellow lines hint at Python interaction.
- Click on a line that starts with a "+" to see the C code that Cython generated for it. -

-

Raw output: lab_xs.c

-
 01: from libc.math cimport tanh, atanh, sqrt
-
 02: 
-
+03: def diff_xs_eta(double e_proton, double charge, double eta, double x_1, double x_2):
-
/* Python wrapper */
-static PyObject *__pyx_pw_6lab_xs_1diff_xs_eta(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
-static char __pyx_doc_6lab_xs_diff_xs_eta[] = "Calculates the differential cross section as a function of the\n    cosine of the pseudo rapidity eta of one photon in units of 1/GeV\302\262.\n    Here d\316\251=detad\317\206\n    :param e_proton: proton energy per beam [GeV]\n    :param charge: charge of the quark\n    :param x_1: momentum fraction of the first quark\n    :param x_2: momentum fraction of the second quark\n    :param eta: pseudo rapidity\n    :return: the differential cross section [GeV^{-2}]\n    ";
-static PyMethodDef __pyx_mdef_6lab_xs_1diff_xs_eta = {"diff_xs_eta", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_6lab_xs_1diff_xs_eta, METH_VARARGS|METH_KEYWORDS, __pyx_doc_6lab_xs_diff_xs_eta};
-static PyObject *__pyx_pw_6lab_xs_1diff_xs_eta(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
-  double __pyx_v_e_proton;
-  double __pyx_v_charge;
-  double __pyx_v_eta;
-  double __pyx_v_x_1;
-  double __pyx_v_x_2;
-  PyObject *__pyx_r = 0;
-  __Pyx_RefNannyDeclarations
-  __Pyx_RefNannySetupContext("diff_xs_eta (wrapper)", 0);
-  {
-    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_e_proton,&__pyx_n_s_charge,&__pyx_n_s_eta,&__pyx_n_s_x_1,&__pyx_n_s_x_2,0};
-    PyObject* values[5] = {0,0,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  5: values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
-        CYTHON_FALLTHROUGH;
-        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
-        CYTHON_FALLTHROUGH;
-        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_e_proton)) != 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_charge)) != 0)) kw_args--;
-        else {
-          __Pyx_RaiseArgtupleInvalid("diff_xs_eta", 1, 5, 5, 1); __PYX_ERR(0, 3, __pyx_L3_error)
-        }
-        CYTHON_FALLTHROUGH;
-        case  2:
-        if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_eta)) != 0)) kw_args--;
-        else {
-          __Pyx_RaiseArgtupleInvalid("diff_xs_eta", 1, 5, 5, 2); __PYX_ERR(0, 3, __pyx_L3_error)
-        }
-        CYTHON_FALLTHROUGH;
-        case  3:
-        if (likely((values[3] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_x_1)) != 0)) kw_args--;
-        else {
-          __Pyx_RaiseArgtupleInvalid("diff_xs_eta", 1, 5, 5, 3); __PYX_ERR(0, 3, __pyx_L3_error)
-        }
-        CYTHON_FALLTHROUGH;
-        case  4:
-        if (likely((values[4] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_x_2)) != 0)) kw_args--;
-        else {
-          __Pyx_RaiseArgtupleInvalid("diff_xs_eta", 1, 5, 5, 4); __PYX_ERR(0, 3, __pyx_L3_error)
-        }
-      }
-      if (unlikely(kw_args > 0)) {
-        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "diff_xs_eta") < 0)) __PYX_ERR(0, 3, __pyx_L3_error)
-      }
-    } else if (PyTuple_GET_SIZE(__pyx_args) != 5) {
-      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);
-      values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
-      values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
-    }
-    __pyx_v_e_proton = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_e_proton == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)
-    __pyx_v_charge = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_charge == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)
-    __pyx_v_eta = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_eta == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)
-    __pyx_v_x_1 = __pyx_PyFloat_AsDouble(values[3]); if (unlikely((__pyx_v_x_1 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)
-    __pyx_v_x_2 = __pyx_PyFloat_AsDouble(values[4]); if (unlikely((__pyx_v_x_2 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)
-  }
-  goto __pyx_L4_argument_unpacking_done;
-  __pyx_L5_argtuple_error:;
-  __Pyx_RaiseArgtupleInvalid("diff_xs_eta", 1, 5, 5, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 3, __pyx_L3_error)
-  __pyx_L3_error:;
-  __Pyx_AddTraceback("lab_xs.diff_xs_eta", __pyx_clineno, __pyx_lineno, __pyx_filename);
-  __Pyx_RefNannyFinishContext();
-  return NULL;
-  __pyx_L4_argument_unpacking_done:;
-  __pyx_r = __pyx_pf_6lab_xs_diff_xs_eta(__pyx_self, __pyx_v_e_proton, __pyx_v_charge, __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_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) {
-  double __pyx_v_rap;
-  PyObject *__pyx_r = NULL;
-  __Pyx_RefNannyDeclarations
-  __Pyx_RefNannySetupContext("diff_xs_eta", 0);
-/* … */
-  /* function exit code */
-  __pyx_L1_error:;
-  __Pyx_XDECREF(__pyx_t_1);
-  __Pyx_AddTraceback("lab_xs.diff_xs_eta", __pyx_clineno, __pyx_lineno, __pyx_filename);
-  __pyx_r = NULL;
-  __pyx_L0:;
-  __Pyx_XGIVEREF(__pyx_r);
-  __Pyx_RefNannyFinishContext();
-  return __pyx_r;
-}
-/* … */
-  __pyx_tuple_ = PyTuple_Pack(6, __pyx_n_s_e_proton, __pyx_n_s_charge, __pyx_n_s_eta, __pyx_n_s_x_1, __pyx_n_s_x_2, __pyx_n_s_rap); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 3, __pyx_L1_error)
-  __Pyx_GOTREF(__pyx_tuple_);
-  __Pyx_GIVEREF(__pyx_tuple_);
-/* … */
-  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6lab_xs_1diff_xs_eta, NULL, __pyx_n_s_lab_xs); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 3, __pyx_L1_error)
-  __Pyx_GOTREF(__pyx_t_1);
-  if (PyDict_SetItem(__pyx_d, __pyx_n_s_diff_xs_eta, __pyx_t_1) < 0) __PYX_ERR(0, 3, __pyx_L1_error)
-  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-  __pyx_codeobj__2 = (PyObject*)__Pyx_PyCode_New(5, 0, 6, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple_, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_lab_xs_pyx, __pyx_n_s_diff_xs_eta, 3, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__2)) __PYX_ERR(0, 3, __pyx_L1_error)
-
 04:     """Calculates the differential cross section as a function of the
-
 05:     cosine of the pseudo rapidity eta of one photon in units of 1/GeV².
-
 06:     Here dΩ=detadφ
-
 07:     :param e_proton: proton energy per beam [GeV]
-
 08:     :param charge: charge of the quark
-
 09:     :param x_1: momentum fraction of the first quark
-
 10:     :param x_2: momentum fraction of the second quark
-
 11:     :param eta: pseudo rapidity
-
 12:     :return: the differential cross section [GeV^{-2}]
-
 13:     """
-
+14:     cdef double rap = atanh((x_1 - x_2) / (x_1 + x_2))
-
  __pyx_v_rap = atanh(((__pyx_v_x_1 - __pyx_v_x_2) / (__pyx_v_x_1 + __pyx_v_x_2)));
-
+15:     return (
-
  __Pyx_XDECREF(__pyx_r);
-
 16:         charge ** 4
-
 17:         / (137.036 * e_proton) ** 2
-
 18:         / (24 * x_1 * x_2)
-
+19:         * ((tanh(eta - rap)) ** 2 + 1)
-
  __pyx_t_1 = PyFloat_FromDouble((((pow(__pyx_v_charge, 4.0) / pow((137.036 * __pyx_v_e_proton), 2.0)) / ((24.0 * __pyx_v_x_1) * __pyx_v_x_2)) * (pow(tanh((__pyx_v_eta - __pyx_v_rap)), 2.0) + 1.0))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 19, __pyx_L1_error)
-  __Pyx_GOTREF(__pyx_t_1);
-  __pyx_r = __pyx_t_1;
-  __pyx_t_1 = 0;
-  goto __pyx_L0;
-
 20:     )
-
 21: 
-
+22: def pT(double e_proton, double eta, double x_1, double x_2):
-
/* Python wrapper */
-static PyObject *__pyx_pw_6lab_xs_3pT(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
-static PyMethodDef __pyx_mdef_6lab_xs_3pT = {"pT", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_6lab_xs_3pT, METH_VARARGS|METH_KEYWORDS, 0};
-static PyObject *__pyx_pw_6lab_xs_3pT(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
-  double __pyx_v_e_proton;
-  double __pyx_v_eta;
-  double __pyx_v_x_1;
-  double __pyx_v_x_2;
-  PyObject *__pyx_r = 0;
-  __Pyx_RefNannyDeclarations
-  __Pyx_RefNannySetupContext("pT (wrapper)", 0);
-  {
-    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_e_proton,&__pyx_n_s_eta,&__pyx_n_s_x_1,&__pyx_n_s_x_2,0};
-    PyObject* values[4] = {0,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  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
-        CYTHON_FALLTHROUGH;
-        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_e_proton)) != 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_eta)) != 0)) kw_args--;
-        else {
-          __Pyx_RaiseArgtupleInvalid("pT", 1, 4, 4, 1); __PYX_ERR(0, 22, __pyx_L3_error)
-        }
-        CYTHON_FALLTHROUGH;
-        case  2:
-        if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_x_1)) != 0)) kw_args--;
-        else {
-          __Pyx_RaiseArgtupleInvalid("pT", 1, 4, 4, 2); __PYX_ERR(0, 22, __pyx_L3_error)
-        }
-        CYTHON_FALLTHROUGH;
-        case  3:
-        if (likely((values[3] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_x_2)) != 0)) kw_args--;
-        else {
-          __Pyx_RaiseArgtupleInvalid("pT", 1, 4, 4, 3); __PYX_ERR(0, 22, __pyx_L3_error)
-        }
-      }
-      if (unlikely(kw_args > 0)) {
-        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "pT") < 0)) __PYX_ERR(0, 22, __pyx_L3_error)
-      }
-    } else if (PyTuple_GET_SIZE(__pyx_args) != 4) {
-      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);
-      values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
-    }
-    __pyx_v_e_proton = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_e_proton == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 22, __pyx_L3_error)
-    __pyx_v_eta = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_eta == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 22, __pyx_L3_error)
-    __pyx_v_x_1 = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_x_1 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 22, __pyx_L3_error)
-    __pyx_v_x_2 = __pyx_PyFloat_AsDouble(values[3]); if (unlikely((__pyx_v_x_2 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 22, __pyx_L3_error)
-  }
-  goto __pyx_L4_argument_unpacking_done;
-  __pyx_L5_argtuple_error:;
-  __Pyx_RaiseArgtupleInvalid("pT", 1, 4, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 22, __pyx_L3_error)
-  __pyx_L3_error:;
-  __Pyx_AddTraceback("lab_xs.pT", __pyx_clineno, __pyx_lineno, __pyx_filename);
-  __Pyx_RefNannyFinishContext();
-  return NULL;
-  __pyx_L4_argument_unpacking_done:;
-  __pyx_r = __pyx_pf_6lab_xs_2pT(__pyx_self, __pyx_v_e_proton, __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_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) {
-  double __pyx_v_tanh_eta;
-  PyObject *__pyx_r = NULL;
-  __Pyx_RefNannyDeclarations
-  __Pyx_RefNannySetupContext("pT", 0);
-/* … */
-  /* function exit code */
-  __pyx_L1_error:;
-  __Pyx_XDECREF(__pyx_t_1);
-  __Pyx_AddTraceback("lab_xs.pT", __pyx_clineno, __pyx_lineno, __pyx_filename);
-  __pyx_r = NULL;
-  __pyx_L0:;
-  __Pyx_XGIVEREF(__pyx_r);
-  __Pyx_RefNannyFinishContext();
-  return __pyx_r;
-}
-/* … */
-  __pyx_tuple__3 = PyTuple_Pack(5, __pyx_n_s_e_proton, __pyx_n_s_eta, __pyx_n_s_x_1, __pyx_n_s_x_2, __pyx_n_s_tanh_eta); if (unlikely(!__pyx_tuple__3)) __PYX_ERR(0, 22, __pyx_L1_error)
-  __Pyx_GOTREF(__pyx_tuple__3);
-  __Pyx_GIVEREF(__pyx_tuple__3);
-/* … */
-  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_6lab_xs_3pT, NULL, __pyx_n_s_lab_xs); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 22, __pyx_L1_error)
-  __Pyx_GOTREF(__pyx_t_1);
-  if (PyDict_SetItem(__pyx_d, __pyx_n_s_pT, __pyx_t_1) < 0) __PYX_ERR(0, 22, __pyx_L1_error)
-  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-  __pyx_codeobj__4 = (PyObject*)__Pyx_PyCode_New(4, 0, 5, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__3, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_lab_xs_pyx, __pyx_n_s_pT, 22, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__4)) __PYX_ERR(0, 22, __pyx_L1_error)
-
+23:     cdef double tanh_eta = tanh(eta)
-
  __pyx_v_tanh_eta = tanh(__pyx_v_eta);
-
+24:     return (
-
  __Pyx_XDECREF(__pyx_r);
-
 25:         2
-
 26:         * e_proton
-
 27:         * x_1
-
 28:         * x_2
-
 29:         / (x_1 + x_2 - (x_1 - x_2) * tanh_eta)
-
+30:         * sqrt(1 - tanh_eta ** 2)
-
  __pyx_t_1 = PyFloat_FromDouble((((((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))))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 30, __pyx_L1_error)
-  __Pyx_GOTREF(__pyx_t_1);
-  __pyx_r = __pyx_t_1;
-  __pyx_t_1 = 0;
-  goto __pyx_L0;
-
 31:     )
-
 32: 
-
 33: 
-
+34: def averaged_tchanel_q2(double e_proton, double x_1, double x_2):
-
/* 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) {
-  double __pyx_v_e_proton;
-  double __pyx_v_x_1;
-  double __pyx_v_x_2;
-  PyObject *__pyx_r = 0;
-  __Pyx_RefNannyDeclarations
-  __Pyx_RefNannySetupContext("averaged_tchanel_q2 (wrapper)", 0);
-  {
-    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_e_proton,&__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_e_proton)) != 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("averaged_tchanel_q2", 1, 3, 3, 1); __PYX_ERR(0, 34, __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("averaged_tchanel_q2", 1, 3, 3, 2); __PYX_ERR(0, 34, __pyx_L3_error)
-        }
-      }
-      if (unlikely(kw_args > 0)) {
-        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "averaged_tchanel_q2") < 0)) __PYX_ERR(0, 34, __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_e_proton = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_e_proton == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 34, __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, 34, __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, 34, __pyx_L3_error)
-  }
-  goto __pyx_L4_argument_unpacking_done;
-  __pyx_L5_argtuple_error:;
-  __Pyx_RaiseArgtupleInvalid("averaged_tchanel_q2", 1, 3, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 34, __pyx_L3_error)
-  __pyx_L3_error:;
-  __Pyx_AddTraceback("lab_xs.averaged_tchanel_q2", __pyx_clineno, __pyx_lineno, __pyx_filename);
-  __Pyx_RefNannyFinishContext();
-  return NULL;
-  __pyx_L4_argument_unpacking_done:;
-  __pyx_r = __pyx_pf_6lab_xs_4averaged_tchanel_q2(__pyx_self, __pyx_v_e_proton, __pyx_v_x_1, __pyx_v_x_2);
-
-  /* function exit code */
-  __Pyx_RefNannyFinishContext();
-  return __pyx_r;
-}
-
-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) {
-  PyObject *__pyx_r = NULL;
-  __Pyx_RefNannyDeclarations
-  __Pyx_RefNannySetupContext("averaged_tchanel_q2", 0);
-/* … */
-  /* function exit code */
-  __pyx_L1_error:;
-  __Pyx_XDECREF(__pyx_t_1);
-  __Pyx_AddTraceback("lab_xs.averaged_tchanel_q2", __pyx_clineno, __pyx_lineno, __pyx_filename);
-  __pyx_r = NULL;
-  __pyx_L0:;
-  __Pyx_XGIVEREF(__pyx_r);
-  __Pyx_RefNannyFinishContext();
-  return __pyx_r;
-}
-/* … */
-  __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_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;
-  __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)
-
+35:     return 2 * x_1 * x_2 * e_proton ** 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)
-  __Pyx_GOTREF(__pyx_t_1);
-  __pyx_r = __pyx_t_1;
-  __pyx_t_1 = 0;
-  goto __pyx_L0;
-
 36: 
-
+37: def second_eta(double eta, double x_1, double x_2):
-
/* 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
-  __Pyx_RefNannySetupContext("second_eta", 0);
-/* … */
-  /* 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;
-}
-/* … */
-  __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_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;
-
+38:     cdef double rap = atanh((x_1 - x_2) / (x_1 + x_2))
-
  __pyx_v_rap = atanh(((__pyx_v_x_1 - __pyx_v_x_2) / (__pyx_v_x_1 + __pyx_v_x_2)));
-
+39:     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;
-
diff --git a/prog/python/qqgg/monte_carlo.py b/prog/python/qqgg/monte_carlo.py index 6aeea6f..fbd76ba 100644 --- a/prog/python/qqgg/monte_carlo.py +++ b/prog/python/qqgg/monte_carlo.py @@ -741,7 +741,7 @@ def sample_unweighted_array( @utility.numpy_cache("cache") def integrate_vegas_nd( f, - interval, + intervals, seed=None, num_increments=5, epsilon=1e-3, @@ -779,12 +779,12 @@ def integrate_vegas_nd( :rtype: tuple """ - intervals = np.asarray(_process_interval(interval)) - ndim = len(interval) + intervals = np.asarray(_process_interval(intervals)) + ndim = len(intervals) integration_volume = (intervals[:, 1] - intervals[:, 0]).prod() if not isinstance(num_increments, collections.Iterable): - num_increments = np.ones(ndim) * num_increments + num_increments = np.ones(ndim).astype(int) * num_increments else: num_increments = np.asarray(num_increments) @@ -841,9 +841,16 @@ def integrate_vegas_nd( while True: vegas_iterations += 1 new_increment_borders = [] + remainder = 0 for dim in range(ndim): increment_weights = generate_integral_steps(increment_borders.copy(), dim) + nonzero_increments = increment_weights[increment_weights > 0] + + remainder += (nonzero_increments.max() - nonzero_increments.min()) * len( + nonzero_increments + ) + new_borders = reshuffle_increments_nd( increment_weights, num_increments[dim], @@ -853,18 +860,7 @@ def integrate_vegas_nd( ) new_increment_borders.append(new_borders) - remainder = np.array( - [ - np.abs( - (borders[1:-1] - old_borders[1:-1]) / (borders[0] - borders[-1]) - ).max() - for borders, old_borders in zip( - new_increment_borders, increment_borders - ) - ] - ).max() - - print(remainder) + remainder /= ndim increment_borders = new_increment_borders if abs(remainder) < increment_epsilon: break diff --git a/prog/python/qqgg/vegas_diagnostics.py b/prog/python/qqgg/vegas_diagnostics.py new file mode 100644 index 0000000..8051ba7 --- /dev/null +++ b/prog/python/qqgg/vegas_diagnostics.py @@ -0,0 +1,31 @@ +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from matplotlib import cm +from monte_carlo import * + + +def plot_cubes(f, increments): + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + cubes = generate_cubes(increments) + for cube in cubes: + x_range, y_range = cube + xs, ys = np.mgrid[ + x_range[0] : x_range[1] : 0.01, y_range[0] : y_range[1] : 0.01, + ] + + vol = get_integration_volume(cube) + fun = f(xs, ys) + + ax.plot_surface(xs, ys, fun * vol * len(cubes), cmap=cm.coolwarm, linewidth=0) + ax.plot_surface(xs, ys, fun, cmap=cm.coolwarm, linewidth=0) + # ax.plot_surface(xs, ys, np.ones_like(fun) / vol, cmap=cm.coolwarm, linewidth=0) + + return fig, ax + + +def test_f(f, *args, **kwargs): + res = integrate_vegas_nd(f, *args, **kwargs) + fig, ax = plot_cubes(f, res.increment_borders) + plt.show() + return res