Merge branch 'master' of github.com:APN-Pucky/HEPi

distribute/tmp.tmp

@@ -0,0 +1,55 @@
+   1        -1.0635e-12   7.63e-13          -1.0635e-12      7.63e-13    0.00e+00
+   2        -6.9063e-13   7.66e-13          -8.7766e-13      5.41e-13    1.19e-01
+   3        -8.9693e-13   7.65e-13          -8.8407e-13      4.42e-13    5.97e-02
+   3        -9.8977e-13   7.66e-13          -9.1043e-13      3.83e-13    4.45e-02
+   3        -1.9127e-12   7.67e-13          -1.1100e-12      3.42e-13    3.75e-01
+   3        -2.0852e-12   7.65e-13          -1.2725e-12      3.13e-13    5.71e-01
+   3        -1.3260e-12   7.69e-13          -1.2801e-12      2.90e-13    4.76e-01
+   3        -4.8079e-13   7.69e-13          -1.1810e-12      2.71e-13    5.43e-01
+   3        -8.5426e-13   7.75e-13          -1.1454e-12      2.56e-13    4.95e-01
+   3        -1.4276e-12   7.73e-13          -1.1732e-12      2.43e-13    4.53e-01
+   3        1.3890e-13   7.73e-13          -1.0554e-12      2.32e-13    6.70e-01
+   3        -2.2356e-12   7.77e-13          -1.1517e-12      2.22e-13    8.02e-01
+   3        -2.1872e-12   7.78e-13          -1.2297e-12      2.14e-13    8.71e-01
+   3        -3.4192e-12   7.78e-13          -1.3832e-12      2.06e-13    1.37e+00
+   3        -1.7759e-12   7.85e-13          -1.4085e-12      1.99e-13    1.29e+00
+   3        4.9687e-13   7.86e-13          -1.2935e-12      1.93e-13    1.57e+00
+   3        -2.1280e-12   7.83e-13          -1.3414e-12      1.87e-13    1.54e+00
+   3        -6.2883e-13   7.82e-13          -1.3027e-12      1.82e-13    1.50e+00
+   3        -1.8292e-12   7.85e-13          -1.3296e-12      1.78e-13    1.44e+00
+   3        -3.4403e-14   7.91e-13          -1.2675e-12      1.73e-13    1.50e+00
+   3        -2.9205e-12   7.90e-13          -1.3433e-12      1.69e-13    1.63e+00
+   3        -6.3022e-13   7.95e-13          -1.3124e-12      1.66e-13    1.59e+00
+   3        -2.4205e-12   8.01e-13          -1.3578e-12      1.62e-13    1.60e+00
+   3        -1.5846e-12   8.04e-13          -1.3666e-12      1.59e-13    1.53e+00
+   3        -1.2509e-12   8.06e-13          -1.3623e-12      1.56e-13    1.47e+00
+   3        -1.8451e-12   8.00e-13          -1.3800e-12      1.53e-13    1.43e+00
+   3        -1.3302e-12   8.03e-13          -1.3782e-12      1.50e-13    1.37e+00
+   3        -1.2982e-12   8.12e-13          -1.3756e-12      1.48e-13    1.32e+00
+   3        -2.4247e-12   8.13e-13          -1.4091e-12      1.45e-13    1.33e+00
+   3        -6.7669e-13   8.27e-13          -1.3872e-12      1.43e-13    1.31e+00
+   3        -1.6549e-12   8.18e-13          -1.3951e-12      1.41e-13    1.27e+00
+   3        -5.0768e-13   8.27e-13          -1.3701e-12      1.39e-13    1.27e+00
+   3        -2.8081e-12   8.33e-13          -1.4091e-12      1.37e-13    1.32e+00
+   3        -3.1958e-12   8.21e-13          -1.4576e-12      1.35e-13    1.42e+00
+   3        -4.6863e-13   8.22e-13          -1.4316e-12      1.33e-13    1.42e+00
+   3        -1.4441e-12   8.26e-13          -1.4319e-12      1.32e-13    1.38e+00
+   3        -1.1787e-12   8.29e-13          -1.4256e-12      1.30e-13    1.34e+00
+   3        -2.6902e-12   8.32e-13          -1.4558e-12      1.29e-13    1.37e+00
+   3        -9.9002e-13   8.34e-13          -1.4450e-12      1.27e-13    1.34e+00
+   3        -1.0315e-12   8.36e-13          -1.4357e-12      1.26e-13    1.31e+00
+   3        -1.6542e-12   8.38e-13          -1.4405e-12      1.24e-13    1.28e+00
+   3        5.1905e-13   8.42e-13          -1.3987e-12      1.23e-13    1.38e+00
+   3        -1.8640e-12   8.41e-13          -1.4085e-12      1.22e-13    1.35e+00
+   3        -8.8271e-13   8.41e-13          -1.3977e-12      1.20e-13    1.33e+00
+   3        -1.4344e-12   8.49e-13          -1.3984e-12      1.19e-13    1.30e+00
+   3        -1.9147e-12   8.44e-13          -1.4085e-12      1.18e-13    1.28e+00
+   3        -2.0361e-12   8.49e-13          -1.4204e-12      1.17e-13    1.26e+00
+   3        -2.1018e-12   8.52e-13          -1.4330e-12      1.16e-13    1.25e+00
+   3        -1.5627e-12   8.45e-13          -1.4354e-12      1.15e-13    1.22e+00
+   3        -1.6809e-12   8.53e-13          -1.4397e-12      1.14e-13    1.20e+00
+   3        -1.9045e-12   8.48e-13          -1.4479e-12      1.13e-13    1.18e+00
+   3        -1.6423e-12   8.59e-13          -1.4512e-12      1.12e-13    1.16e+00
+   3        -1.6388e-12   8.56e-13          -1.4544e-12      1.11e-13    1.14e+00
+NLO+NLL = (1.2043830e-11 +- 1.1167217e-13) pb
+

docs/source/examples/test_interpolate.ipynb

@@ -1,5 +1,13 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "41876a00",
+   "metadata": {},
+   "source": [
+    "# Interpolation"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 1,
@@ -150,7 +158,7 @@
   {
    "cell_type": "code",
    "execution_count": 6,
-   "id": "8d6dd918",
+   "id": "b14f324c",
    "metadata": {},
    "outputs": [
     {
@@ -177,7 +185,7 @@
   {
    "cell_type": "code",
    "execution_count": 7,
-   "id": "c7a2d5ba",
+   "id": "ba1d705f",
    "metadata": {},
    "outputs": [
     {
@@ -215,7 +223,7 @@
   {
    "cell_type": "code",
    "execution_count": 10,
-   "id": "6ed7f08e",
+   "id": "28d178a5",
    "metadata": {},
    "outputs": [
     {
@@ -249,7 +257,7 @@
   {
    "cell_type": "code",
    "execution_count": 11,
-   "id": "d4c5cf40",
+   "id": "efc9b48b",
    "metadata": {},
    "outputs": [
     {
@@ -525,7 +533,7 @@
   {
    "cell_type": "code",
    "execution_count": 26,
-   "id": "ce6b8b44",
+   "id": "91517bf5",
    "metadata": {
     "scrolled": false
    },
@@ -573,7 +581,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "3bfbc278",
+   "id": "5f5c40b5",
    "metadata": {},
    "outputs": [],
    "source": []

docs/source/examples/test_madgraph.ipynb

@@ -729,7 +729,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },

hepi/plot.py

@@ -18,8 +18,8 @@
 from matplotlib.ticker import NullFormatter
 
 
-def title(axe,
-          i: Input,
+def title(i: Input,
+          axe=None,
           scenario="",
           diff_L_R=None,
           extra="",
@@ -28,6 +28,8 @@ def title(axe,
           id=False,
           **kwargs):
     """Sets the title on axis `axe`."""
+    if axe is None:
+        axe = plt.gca()
     axe.set_title("$pp\\to" + get_name(i.particle1) + get_name(i.particle2) +
                   "$" + (" at $\\sqrt{S} = " + str(i.energy / 1000) +
                          "$ TeV" if cms_energy else "") + " for " +

hepi/results.py

@@ -72,8 +72,94 @@ def __init__(self,
         else:
             self.aNNLO_PLUS_NNLL_OVER_NLO = None
 
+def pdf_errors(li, dl,ordernames=["LO","NLO","aNNLO_PLUS_NNLL"], confidence_level=90):
+    """
+    Just like `pdf_error` but over a list of ordernames.
+    """
+    r_dl = dl
+    for o in ordernames:
+        r_dl = pdf_error(li,r_dl,o,confidence_level=confidence_level)
+
+    return r_dl
+
+def pdf_error(li, dl,ordername="LO", confidence_level=90):
+    """
+    Computes Parton Density Function (PDF) uncertainties through :func:`lhapdf.set.uncertainty`.
+
+    Args:
+        li (:obj:`list` of :class:`Input`): Input list.
+        dl (:obj:`dict`): :class:`Result` dictionary with lists per entry.
+        ordername (`str`): Name of the order.
+        confidence_level (:obj:`double`): Confidence Level for PDF uncertainty
+
+    Returns:
+        :obj:`dict`: Modified `dl` with new `LO`/`NLO`/`NLO_PLUS_NLL` _ `PDF`/`PDF_CENTRAL`/`PDF_ERRPLUS`/`PDF_ERRMINUS`/`PDF_ERRSYM` entries.
+            - `LO`/`NLO`/`NLO_PLUS_NLL` _ `PDF` contains a symmetrized :mod:`uncertainties` object.
+    """
+    global required_numerical_uncertainty_factor
+    example = li[0]
+    members = [
+        attr for attr in dir(example)
+        if not callable(getattr(example, attr)) and not attr.startswith("__")
+    ]
+
+    dl[ordername + "_PDF"] = np.array([None] * len(dl["pdfset_nlo"]))
+    dl[ordername+"_PDF_CENTRAL"] = np.array([None] * len(dl["pdfset_nlo"]))
+    dl[ordername+"_PDF_ERRPLUS"] = np.array([None] * len(dl["pdfset_nlo"]))
+    dl[ordername+"_PDF_ERRMINUS"] = np.array([None] * len(dl["pdfset_nlo"]))
+    dl[ordername+"_PDF_ERRSYM"] = np.array([None] * len(dl["pdfset_nlo"]))
+
+    for i in range(len(dl["pdfset_nlo"])):
+        if dl["pdfset_nlo"][i] == 0 and dl["mu_f"][i] == 1.0 and dl["mu_r"][
+                i] == 1.0:
+            pdfset = lhapdf.getPDFSet(dl["pdf_nlo"][i])
+            pdfs = [0.0] * pdfset.size
+            for j in range(len(dl["pdfset_nlo"])):
+                same = True
+                for s in members:
+                    if not (
+                            dl[s][i] == dl[s][j]
+                    ) and s != "pdfset_nlo" and s != "precision" and s != "max_iters":
+                        same = False
+                if same:
+                    pdfs[dl["pdfset_nlo"][j]] = j
+
+           # lo_unc = pdfset.uncertainty(
+            #    [plot.unv(dl["LO"][k]) for k in pdfs], -1)
+            if ordername=="LO":
+                dl.loc[i,"LO_PDF_CENTRAL"] = plot.unv(dl["LO"][i])
+                dl.loc[i,"LO_PDF_ERRPLUS"] = 0.0
+                dl.loc[i,"LO_PDF_ERRMINUS"] = 0.0
+                dl.loc[i,"LO_PDF_ERRSYM"] = 0.0
+            else:
+                nlo_unc = pdfset.uncertainty(
+                    [plot.unv(dl["NLO"][k]) for k in pdfs], confidence_level)
+                dl.loc[i,ordername + "_PDF_CENTRAL"] = nlo_unc.central
+                dl.loc[i,ordername + "_PDF_ERRPLUS"] = nlo_unc.errplus
+                dl.loc[i,ordername + "_PDF_ERRMINUS"] = -nlo_unc.errminus
+                dl.loc[i,ordername + "_PDF_ERRSYM"] = nlo_unc.errsymm
+            #TODO error sym to minus and plus
+            if ordername!= "LO":
+                if (plot.usd(dl[ordername][i]) * required_numerical_uncertainty_factor
+                    > dl[ordername + "_PDF_ERRPLUS"][i] or plot.usd(dl[ordername][i]) *
+                    required_numerical_uncertainty_factor >
+                    -dl[ordername + "_PDF_ERRMINUS"][i]):
+                    rel = plot.unv(dl[ordername][i])
+                    warnings.warn("too bad numerical precision vs pdf @ "+ordername+ " num: "+
+                                    str(plot.usd(dl[ordername][i])/rel*100.)+ "% vs "+ str(dl[ordername + "_PDF_ERRPLUS"][i]/rel*100.)+ "% to pdf: "+ str(dl[ordername + "_PDF_ERRMINUS"][i]/rel*100. + "%"),
+                            RuntimeWarning)
+
+    mask = dl[ordername+"_PDF_CENTRAL"].notnull() 
+    dl.loc[mask,ordername+"_PDF"] = unumpy.uarray(
+        plot.unv(dl[ordername+"_PDF_CENTRAL"][mask]) +
+        dl[ordername+"_PDF_ERRPLUS"][mask] / 2. + dl[ordername+"_PDF_ERRMINUS"][mask] / 2.,
+        (+dl[ordername+"_PDF_ERRPLUS"][mask] - dl[ordername+"_PDF_ERRMINUS"][mask]) / 2.)
+
+    return dl
 
-def pdf_error(li, dl, confidence_level=90):
+
+
+def pdf_error_old(li, dl, confidence_level=90):
     """
     Computes Parton Density Function (PDF) uncertainties through :func:`lhapdf.set.uncertainty`.
 
@@ -180,8 +266,80 @@ def pdf_error(li, dl, confidence_level=90):
          dl["NLO_PLUS_NLL_PDF_ERRMINUS"][mask]) / 2.)
     return dl
 
+def scale_errors(li, dl,ordernames=["LO","NLO","aNNLO_PLUS_NNLL"]):
+    """
+    Just like `scale_error` but over a list of ordernames.
+    """
+    r_dl = dl
+    for o in ordernames:
+        r_dl = scale_error(li,r_dl,o)
+
+    return r_dl
+def scale_error(li, dl,ordername="LO"):
+    """
+    Computes seven-point scale uncertainties from the results where the renormalization and factorization scales are varied by factors of 2 and  relative factors of four are excluded (cf. :meth:`seven_point_scan`).
+
+    Args:
+        li (:obj:`list` of :class:`Input`): Input list.
+        dl (:obj:`dict`): :class:`Result` dictionary with lists per entry.
+
+    Returns:
+        :obj:`dict`: Modified `dl` with new `LO`/`NLO`/`NLO_PLUS_NLL` _ `SCALE`/`SCALE_ERRPLUS`/`SCALE_ERRMINUS`/`SCALE_ERRSYM` entries.
+            - `LO`/`NLO`/`NLO_PLUS_NLL` _ `SCALE` contains a symmetrized :mod:`uncertainties` object.
+    """
+    global required_numerical_uncertainty_factor
+    example = li[0]
+    members = [
+        attr for attr in dir(example)
+        if not callable(getattr(example, attr)) and not attr.startswith("__")
+    ]
+    dl[ordername+"_SCALE"] = np.array([None] * len(dl["pdfset_nlo"]))
+    dl[ordername+"_SCALE_ERRPLUS"] = np.array([None] * len(dl["pdfset_nlo"]))
+    dl[ordername+"_SCALE_ERRMINUS"] = np.array([None] * len(dl["pdfset_nlo"]))
+
+    for i in range(len(dl["pdfset_nlo"])):
+        if dl["pdfset_nlo"][i] == 0 and dl["mu_f"][i] == 1.0 and dl["mu_r"][
+                i] == 1.0:
+            scales = []
+            for j in range(len(dl["pdfset_nlo"])):
+                same = True
+                for s in members:
+                    if not (
+                            dl[s][i] == dl[s][j]
+                    ) and s != "mu_f" and s != "mu_r" and s != "precision" and s != "max_iters":
+                        same = False
+                if same:
+                    scales.append(j)
+
+            # lo_unc = pdfset.uncertainty(
+            #    [plot.unv(dl["LO"][k]) for k in pdfs], -1)
+            dl.loc[i,ordername+"_SCALE_ERRPLUS"] = np.max(
+                [plot.unv(dl[ordername][k])
+                 for k in scales]) - plot.unv(dl[ordername][i])
+            dl.loc[i,ordername+"_SCALE_ERRMINUS"] = np.min(
+                [plot.unv(dl[ordername][k])
+                 for k in scales]) - plot.unv(dl[ordername][i])
+            if (plot.usd(dl[ordername][i]) * required_numerical_uncertainty_factor >
+                    dl[ordername+"_SCALE_ERRPLUS"][i] or plot.usd(dl[ordername][i]) *
+                    required_numerical_uncertainty_factor >
+                    -dl[ordername+"_SCALE_ERRMINUS"][i]):
+                rel = plot.unv(dl[ordername][i])
+                warnings.warn("too bad numerical precision vs scale @ num:"+ordername + " " + str(plot.usd(dl[ordername][i])/rel*100.) + "% vs scale:" +str(dl[ordername+"_SCALE_ERRPLUS"][i]/rel*100.) + "% to " + str(dl[ordername+"_SCALE_ERRMINUS"][i]/rel*100.) + "%"  ,  
+                              RuntimeWarning)
+
+
+    mask = dl[ordername+"_SCALE_ERRPLUS"].notnull()
+    dl.loc[mask ,ordername+"_SCALE"] = unumpy.uarray(
+        plot.unv(dl[ordername][mask]) + dl[ordername+"_SCALE_ERRPLUS"][mask] / 2. +
+        dl[ordername+"_SCALE_ERRMINUS"][mask] / 2.,
+        (+dl[ordername+"_SCALE_ERRPLUS"][mask] - dl[ordername+"_SCALE_ERRMINUS"][mask]) / 2.)
+
+
+    return dl
+
 
-def scale_error(li, dl):
+
+def scale_error_old(li, dl):
     """
     Computes seven-point scale uncertainties from the results where the renormalization and factorization scales are varied by factors of 2 and  relative factors of four are excluded (cf. :meth:`seven_point_scan`).
 
@@ -290,7 +448,53 @@ def scale_error(li, dl):
     return dl
 
 
-def combine_errors(dl: dict):
+def combine_errors( dl,ordernames=["LO","NLO","aNNLO_PLUS_NNLL"]):
+    """
+    Just like `combine_error` but over a list of ordernames.
+    """
+    r_dl = dl
+    for o in ordernames:
+        r_dl = combine_error(r_dl,o)
+
+    return r_dl
+
+def combine_error(dl: dict,ordername="LO"):
+    """
+    Combines seven-point scale uncertainties and pdf uncertainties from the results by Pythagorean addition.
+
+    Note:
+        Running :func:`scale_errors` and :func:`pdf_errors` before is necessary.
+
+    Args:
+        dl (:obj:`dict`): :class:`Result` dictionary with lists per entry.
+
+    Returns:
+        :obj:`dict`: Modified `dl` with new `LO`/`NLO`/`NLO_PLUS_NLL` _ `COMBINED`/`ERRPLUS`/`ERRMINUS` entries.
+            - `LO`/`NLO`/`NLO_PLUS_NLL` _ `COMBINED` contains a symmetrized :mod:`uncertainties` object.
+    """
+    dl[ordername+"_NOERR"] = np.array([None] * len(dl["pdfset_nlo"]))
+    dl[ordername+"_ERRPLUS"] = np.array([None] * len(dl["pdfset_nlo"]))
+    dl[ordername+"_ERRMINUS"] = np.array([None] * len(dl["pdfset_nlo"]))
+    dl[ordername+"_COMBINED"] = np.array([None] * len(dl["pdfset_nlo"]))
+
+
+    mask = dl[ordername+"_PDF_CENTRAL"].notnull() 
+    dl.loc[mask,ordername+"_NOERR"] = plot.unv(dl[ordername+""][mask]).astype(float)
+    dl.loc[mask,ordername+"_ERRPLUS"] = np.sqrt(
+        dl[ordername+"_PDF_ERRPLUS"][mask].astype(float)**2 +
+        dl[ordername+"_SCALE_ERRPLUS"][mask].astype(float)**2)
+    dl.loc[mask,ordername+"_ERRMINUS"] = -np.sqrt(dl[ordername+"_PDF_ERRMINUS"][mask].astype(
+        float)**2 + dl[ordername+"_SCALE_ERRMINUS"][mask].astype(float)**2)
+    dl.loc[mask,ordername+"_COMBINED"] = unumpy.uarray(
+        plot.unv(dl[ordername+""][mask]) + dl[ordername+"_ERRPLUS"][mask] / 2. +
+        dl[ordername+"_ERRMINUS"][mask] / 2.,
+        +dl[ordername+"_ERRPLUS"][mask] - dl[ordername+"_ERRMINUS"][mask])
+
+
+    return dl
+
+
+def combine_errors_old(dl: dict):
     """
     Combines seven-point scale uncertainties and pdf uncertainties from the results by Pythagorean addition.