From 5be840069a7f129bd75aa27d57e260f6199f4445 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Tue, 29 Apr 2025 18:28:40 +0200 Subject: [PATCH 1/6] Add preliminary version of 3d plot --- pdfplotter/pdf_set_nuclear.py | 308 +++++++++++++++++++++++++++++++++- 1 file changed, 307 insertions(+), 1 deletion(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index b042d26..32590e3 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -1,13 +1,23 @@ from __future__ import annotations from itertools import zip_longest -from typing import Sequence +from math import log +from typing import Sequence, Literal, Any +#from mpl_toolkits.mplot3d import Axes3D +from mpl_toolkits.mplot3d.art3d import Poly3DCollection +import matplotlib.ticker as mticker +import matplotlib.cm as cm import numpy as np import numpy.typing as npt import pandas as pd +import matplotlib.pyplot as plt +import sympy as sp from pdfplotter.pdf_set import PDFSet +from pdfplotter.util import update_kwargs +from pdfplotter.util import log_tick_formatter +from pdfplotter import util class NuclearPDFSet(PDFSet): @@ -153,3 +163,299 @@ def get( raise ValueError(f"Multiple PDFSets found for Z = {Z}") else: return pdf_set.iloc[0]["pdf_set"] + + def plot_A_dep_3d( + self, + ax: plt.Axes | npt.NDArray[plt.Axes], # pyright: ignore[reportInvalidTypeForm] + A: int | float | list[int | float], + observables: ( + sp.Basic + | npt.NDArray[sp.Basic] # pyright: ignore[reportInvalidTypeForm] + | list[sp.Basic] + ), + Q: float | None = None, + Q2: float | None = None, + colors: list = [], + logA: bool = True, + plot_uncertainty: bool = True, + plot_ratio: bool = False, + pdf_label: Literal["ylabel", "annotate"] = "annotate", + A_label: Literal["legend", "ticks"] = "ticks", + proj_type: Literal["ortho", "persp"] = "ortho", + view_init: tuple[float, float] | list[tuple[float, float]] = (15, -75), + kwargs_theory: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_uncertainty: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_uncertainty_edges: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_title: dict[str, Any] = {}, + kwargs_notation: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_ylabel: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_xlabel: dict[str, Any] = {}, + kwargs_zlabel: dict[str, Any] = {}, + kwargs_legend: dict[str, Any] = {}, + ): + + my_sets = self.pdf_sets + x = self.get(A=my_sets["A"][0]).x_values + + + if Q is None and Q2 is None: + raise ValueError("Please pass either `Q` or `Q2`") + + elif Q is not None and Q2 is not None: + raise ValueError("Please pass either `Q` or `Q2`, not both") + + elif Q is not None: + if Q not in self.get(A=my_sets["A"][0]).Q_values and Q not in np.sqrt( + np.array(self.get(A=my_sets["A"][0]).Q2_values) + ): + raise ValueError( + f"Chosen Q value {Q} was not used for defining nuclear pdf set. \n Please choose Q that was used in initialization" + ) + else: + if ( + Q2 not in self.get(A=my_sets["A"][0]).Q2_values + and Q2 not in np.array(self.get(A=my_sets["A"][0]).Q_values) ** 2 + ): + raise ValueError( + f"Chosen Q2 value {Q2} was not used for defining nuclear pdf set. \n Please choose Q2 that was used in initialization" + ) + + if not isinstance(A, list): + A = [A] + + if isinstance(observables, np.ndarray): + observables = list(observables.flatten()) + + if not isinstance(observables, list): + observables = [observables] + + if not isinstance(ax, np.ndarray): + ax = np.array([ax]) + + if colors == []: + cmap = cm.get_cmap("tab10", lut=len(A)) + colors = [cmap(i) for i in range(len(A))] + + for i, (obs_i, ax_i) in enumerate(zip(observables, ax.flat)): + + ax_i.set_proj_type(proj_type) + ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) + for j, (A_j, col_j) in enumerate(zip(A, colors)): + z_lower, z_upper = self.get(A=A_j).get_uncertainties( + observable=obs_i, x=x, Q=Q, Q2=Q2 + ) + kwargs_default = { + "color": col_j, + "label": f"A={A_j}", + "linewidth": 1.5, + } + if not isinstance(kwargs_theory, list): + kwargs = update_kwargs( + kwargs_default, + kwargs_theory, + ) + else: + kwargs = update_kwargs( + kwargs_default, + kwargs_theory, + i=j, + ) + if logA: + ax_i.plot( + np.log10(x), + np.log10(len(x) * [A_j]), + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) + else: + ax_i.plot( + np.log10(x), + len(x) * [A_j], + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) + if plot_uncertainty: + kwargs_uncertainty_default = { + "color": col_j, + "alpha": 0.3, + } + if not isinstance(kwargs_uncertainty, list): + kwargs = update_kwargs( + kwargs_uncertainty_default, + kwargs_uncertainty, + ) + else: + kwargs = update_kwargs( + kwargs_uncertainty_default, + kwargs_uncertainty, + i=j, + ) + + vertices = [] + z_lower = np.array(z_lower) + z_upper = np.array(z_upper) + if not logA: + for xi, ai, zl, zu in zip( + np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper + ): + vertices.append([xi, ai, zl]) + + for xi, ai, zl, zu in reversed( + list(zip(np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper)) + ): + vertices.append([xi, ai, zu]) + else: + for xi, ai, zl, zu in zip( + np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper + ): + vertices.append([xi, ai, zl]) + + for xi, ai, zl, zu in reversed( + list(zip(np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper)) + ): + vertices.append([xi, ai, zu]) + poly = Poly3DCollection([vertices], **kwargs) + ax_i.add_collection3d(poly) + + kwargs_uncertainty_edges_default = { + "color": col_j, + "alpha": 1, + } + if not isinstance(kwargs_uncertainty_edges, list): + kwargs = update_kwargs( + kwargs_uncertainty_edges_default, + kwargs_uncertainty_edges, + ) + else: + kwargs = update_kwargs( + kwargs_uncertainty_edges_default, + kwargs_uncertainty_edges, + i=j, + ) + if not logA: + ax_i.plot(np.log10(x), len(x) * [A_j], z_upper, **kwargs) + ax_i.plot(np.log10(x), len(x) * [A_j], z_lower, **kwargs) + else: + ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_upper, **kwargs) + ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_lower, **kwargs) + if pdf_label == "annotate": + kwargs_notation_default = { + "fontsize": 12, + "xy": (0.97, 0.96), + "xycoords": "axes fraction", + "va": "top", + "ha": "right", + "bbox": dict( + facecolor=(1, 1, 1), + edgecolor=(0.8, 0.8, 0.8), + lw=0.9, + boxstyle="round, pad=0.2", + ), + } + if not isinstance(kwargs_notation, list): + kwargs_n = update_kwargs( + kwargs_notation_default, + kwargs_notation, + ) + else: + kwargs_n = update_kwargs( + kwargs_notation_default, + kwargs_notation, + i=i, + ) + ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) + + if pdf_label == "ylabel": + kwargs_ylabel_default = { + "fontsize": 14, + "zlabel": f"${util.to_str(obs_i,Q=Q,Q2=Q2)}$", + #"labelpad":-200 + } + if not isinstance(kwargs_ylabel, list): + kwargs = update_kwargs( + kwargs_ylabel_default, + kwargs_ylabel, + ) + else: + kwargs = update_kwargs( + kwargs_ylabel_default, + kwargs_ylabel, + i=i, + ) + ax_i.set_zlabel(**kwargs) + + else: + kwargs_notation_default = { + "fontsize": 12, + "xy": (0.47, 0.96), + "xycoords": "axes fraction", + "va": "top", + "ha": "right", + "bbox": dict( + facecolor=(1, 1, 1), + edgecolor=(0.8, 0.8, 0.8), + lw=0.9, + boxstyle="round, pad=0.2", + ), + } + kwargs_n = update_kwargs(kwargs_notation_default, kwargs_notation, i=i) + + ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) + ax_i.xaxis.set_major_formatter(mticker.FuncFormatter(log_tick_formatter)) + if A_label == "ticks": + if logA: + ax_i.set_yticks(np.log10(A),A) + else: + ax_i.set_yticks(A,A) + kwargs_zlabel_default = { + "fontsize": 14, + "ylabel": f"$A$", + + } + kwargs = update_kwargs( + kwargs_zlabel_default, + kwargs_zlabel, + ) + + ax_i.set_ylabel(**kwargs) + else: + ax_i.set_yticks([]) + kwargs_legend_default = { + "fontsize": 12, + "bbox_to_anchor": (0.95, 0.95), + "frameon": False, + } + kwargs = update_kwargs( + kwargs_legend_default, + kwargs_legend, + ) + ax_i.legend() + kwargs_zlabel_default = { + "fontsize": 14, + "ylabel": f"$A$", + "labelpad":-10 + } + kwargs = update_kwargs( + kwargs_zlabel_default, + kwargs_zlabel, + ) + ax_i.set_ylabel(**kwargs) + ax_i.xaxis.pane.fill=False + ax_i.yaxis.pane.fill=False + ax_i.zaxis.pane.fill=False + ax_i.xaxis.pane.set_edgecolor("white") + ax_i.yaxis.pane.set_edgecolor("white") + ax_i.zaxis.pane.set_edgecolor("white") + + ax_i.zaxis._axinfo["juggled"]=(1,2,0) + + kwargs_xlabel_default = { + "fontsize": 14, + "xlabel": f"$x$", + + } + kwargs = update_kwargs( + kwargs_xlabel_default, + kwargs_xlabel, + ) + ax_i.set_xlabel(**kwargs) \ No newline at end of file From 4f620f59e3dadfcac9f0f4390a599d278a36c0ed Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Wed, 30 Apr 2025 10:12:05 +0200 Subject: [PATCH 2/6] Add log_tick_formatter to util.py. Add possibility to plot ratios. --- pdfplotter/pdf_set_nuclear.py | 129 +++++++++++++++++++++++++++------- pdfplotter/util.py | 3 + 2 files changed, 107 insertions(+), 25 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 32590e3..7011a95 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -175,13 +175,14 @@ def plot_A_dep_3d( ), Q: float | None = None, Q2: float | None = None, + x_lines: float | list[float] | None = None, colors: list = [], logA: bool = True, plot_uncertainty: bool = True, plot_ratio: bool = False, pdf_label: Literal["ylabel", "annotate"] = "annotate", A_label: Literal["legend", "ticks"] = "ticks", - proj_type: Literal["ortho", "persp"] = "ortho", + proj_type: Literal["ortho", "persp"] = "persp", view_init: tuple[float, float] | list[tuple[float, float]] = (15, -75), kwargs_theory: dict[str, Any] | list[dict[str, Any] | None] = {}, kwargs_uncertainty: dict[str, Any] | list[dict[str, Any] | None] = {}, @@ -192,6 +193,7 @@ def plot_A_dep_3d( kwargs_xlabel: dict[str, Any] = {}, kwargs_zlabel: dict[str, Any] = {}, kwargs_legend: dict[str, Any] = {}, + kwargs_xlines: dict[str, Any] | list[dict[str, Any] | None] = {}, ): my_sets = self.pdf_sets @@ -223,6 +225,11 @@ def plot_A_dep_3d( if not isinstance(A, list): A = [A] + if 1 not in A and plot_ratio: + raise ValueError( + "Please pass A=1 if you want to plot the ratio to Proton." + ) + if isinstance(observables, np.ndarray): observables = list(observables.flatten()) @@ -238,12 +245,17 @@ def plot_A_dep_3d( for i, (obs_i, ax_i) in enumerate(zip(observables, ax.flat)): - ax_i.set_proj_type(proj_type) - ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) for j, (A_j, col_j) in enumerate(zip(A, colors)): - z_lower, z_upper = self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2 - ) + if not plot_ratio: + z_lower, z_upper = self.get(A=A_j).get_uncertainties( + observable=obs_i, x=x, Q=Q, Q2=Q2 + ) + else: + z_lower, z_upper = self.get(A=A_j).get_uncertainties( + observable=obs_i, x=x, Q=Q, Q2=Q2 + ) + z_lower = z_lower / self.get(A=1).get_central(observable=obs_i,x=x, Q=Q, Q2=Q2) + z_upper = z_upper / self.get(A=1).get_central(observable=obs_i,x=x, Q=Q, Q2=Q2) kwargs_default = { "color": col_j, "label": f"A={A_j}", @@ -261,19 +273,37 @@ def plot_A_dep_3d( i=j, ) if logA: - ax_i.plot( - np.log10(x), - np.log10(len(x) * [A_j]), - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), - **kwargs, - ) + if plot_ratio: + ax_i.plot( + np.log10(x), + np.log10(len(x) * [A_j]), + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i) + / self.get(A=1).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) + else: + ax_i.plot( + np.log10(x), + np.log10(len(x) * [A_j]), + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) else: - ax_i.plot( - np.log10(x), - len(x) * [A_j], - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), - **kwargs, - ) + if plot_ratio: + ax_i.plot( + np.log10(x), + len(x) * [A_j], + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i) + / self.get(A=1).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) + else: + ax_i.plot( + np.log10(x), + len(x) * [A_j], + self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), + **kwargs, + ) if plot_uncertainty: kwargs_uncertainty_default = { "color": col_j, @@ -295,15 +325,16 @@ def plot_A_dep_3d( z_lower = np.array(z_lower) z_upper = np.array(z_upper) if not logA: + for xi, ai, zl, zu in zip( np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper ): vertices.append([xi, ai, zl]) - for xi, ai, zl, zu in reversed( list(zip(np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper)) ): vertices.append([xi, ai, zu]) + else: for xi, ai, zl, zu in zip( np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper @@ -338,6 +369,50 @@ def plot_A_dep_3d( else: ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_upper, **kwargs) ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_lower, **kwargs) + + centrals={} + if x_lines is not None: + if not isinstance(x_lines, list): + x_lines = [x_lines] + for k,x_line in enumerate(x_lines): + if x_line not in x: + raise ValueError( + f"Chosen x value {x_line} was not used for defining nuclear pdf set. \n Please choose x that was used in initialization" + ) + kwargs_xlines_default = { + "color": "black", + "linestyle": "--", + "linewidth": 1.5, + } + if not isinstance(kwargs_xlines, list): + kwargs = update_kwargs( + kwargs_xlines_default, + kwargs_xlines, + ) + else: + kwargs = update_kwargs( + kwargs_xlines_default, + kwargs_xlines, + i=k, + ) + for a in A: + if not plot_ratio: + if x_line not in centrals.keys(): + centrals[x_line] = [self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)] + else: + centrals[x_line].append(self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)) + else: + if x_line not in centrals.keys(): + centrals[x_line] = [self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)/self.get(A=1).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)] + else: + centrals[x_line].append(self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)/self.get(A=1).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)) + if logA: + ax_i.plot( + np.ones(len(A))*np.log10(x_line), np.log10(A), centrals[x_line],**kwargs) + else: + ax_i.plot( + np.ones(len(A))*np.log10(x_line), A, centrals[x_line],**kwargs) + if pdf_label == "annotate": kwargs_notation_default = { "fontsize": 12, @@ -369,7 +444,6 @@ def plot_A_dep_3d( kwargs_ylabel_default = { "fontsize": 14, "zlabel": f"${util.to_str(obs_i,Q=Q,Q2=Q2)}$", - #"labelpad":-200 } if not isinstance(kwargs_ylabel, list): kwargs = update_kwargs( @@ -433,7 +507,8 @@ def plot_A_dep_3d( kwargs_zlabel_default = { "fontsize": 14, "ylabel": f"$A$", - "labelpad":-10 + "labelpad":-10, + "linespacing": -4 } kwargs = update_kwargs( kwargs_zlabel_default, @@ -443,9 +518,9 @@ def plot_A_dep_3d( ax_i.xaxis.pane.fill=False ax_i.yaxis.pane.fill=False ax_i.zaxis.pane.fill=False - ax_i.xaxis.pane.set_edgecolor("white") - ax_i.yaxis.pane.set_edgecolor("white") - ax_i.zaxis.pane.set_edgecolor("white") + ax_i.xaxis.pane.set_edgecolor("w") + ax_i.yaxis.pane.set_edgecolor("w") + ax_i.zaxis.pane.set_edgecolor("w") ax_i.zaxis._axinfo["juggled"]=(1,2,0) @@ -458,4 +533,8 @@ def plot_A_dep_3d( kwargs_xlabel_default, kwargs_xlabel, ) - ax_i.set_xlabel(**kwargs) \ No newline at end of file + ax_i.set_xlabel(**kwargs) + ax_i.set_xlim(np.log10(x[0])-np.log10(x[0])/100, np.log10(x[-1])) + ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 + ax_i.set_proj_type(proj_type) + ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) \ No newline at end of file diff --git a/pdfplotter/util.py b/pdfplotter/util.py index f0478b1..9f6ad64 100644 --- a/pdfplotter/util.py +++ b/pdfplotter/util.py @@ -198,3 +198,6 @@ def update_kwargs( return kwargs else: raise ValueError("kwargs_user must be a dict or a list") + +def log_tick_formatter(val,pos=None): + return r"$10^{{{:.0f}}}$".format(val) \ No newline at end of file From 4825af7bdf0e3764c43d9f967069c9f720c2e109 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Tue, 6 May 2025 10:03:06 +0200 Subject: [PATCH 3/6] Remove y and x axis offset. Add vertical grid lines for case: A_label:yticks --- pdfplotter/pdf_set_nuclear.py | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 7011a95..4c7c021 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -492,18 +492,20 @@ def plot_A_dep_3d( ) ax_i.set_ylabel(**kwargs) + ax_i.set_xlim(np.log10(x[0])*0.98) else: ax_i.set_yticks([]) - kwargs_legend_default = { - "fontsize": 12, - "bbox_to_anchor": (0.95, 0.95), - "frameon": False, - } - kwargs = update_kwargs( - kwargs_legend_default, - kwargs_legend, - ) - ax_i.legend() + if i==len(observables)-1: + kwargs_legend_default = { + "fontsize": 12, + "bbox_to_anchor": (0.95, 0.95), + "frameon": False, + } + kwargs = update_kwargs( + kwargs_legend_default, + kwargs_legend, + ) + ax_i.legend() kwargs_zlabel_default = { "fontsize": 14, "ylabel": f"$A$", @@ -534,7 +536,8 @@ def plot_A_dep_3d( kwargs_xlabel, ) ax_i.set_xlabel(**kwargs) - ax_i.set_xlim(np.log10(x[0])-np.log10(x[0])/100, np.log10(x[-1])) - ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 + #, np.log10(x[-1])) + ax_i.set_zlim(ax_i.get_zlim()[1]*0.02) + #ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 ax_i.set_proj_type(proj_type) ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) \ No newline at end of file From 97e22df840c6853f799f3e98f347eed101681255 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Tue, 6 May 2025 16:12:29 +0200 Subject: [PATCH 4/6] Add plotting function to plot PDFs for different x as scatter --- pdfplotter/pdf_set_nuclear.py | 580 ++++++++++++++++------------------ 1 file changed, 273 insertions(+), 307 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 4c7c021..2e77b80 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -1,6 +1,6 @@ from __future__ import annotations -from itertools import zip_longest +from itertools import zip_longest, cycle from math import log from typing import Sequence, Literal, Any @@ -8,15 +8,19 @@ from mpl_toolkits.mplot3d.art3d import Poly3DCollection import matplotlib.ticker as mticker import matplotlib.cm as cm +import matplotlib.colors as cls +import math import numpy as np import numpy.typing as npt import pandas as pd import matplotlib.pyplot as plt import sympy as sp + from pdfplotter.pdf_set import PDFSet from pdfplotter.util import update_kwargs from pdfplotter.util import log_tick_formatter +from pdfplotter import elements from pdfplotter import util @@ -164,10 +168,10 @@ def get( else: return pdf_set.iloc[0]["pdf_set"] - def plot_A_dep_3d( + def plot_A_dep_2d_scatter( self, ax: plt.Axes | npt.NDArray[plt.Axes], # pyright: ignore[reportInvalidTypeForm] - A: int | float | list[int | float], + x: float | list[float], observables: ( sp.Basic | npt.NDArray[sp.Basic] # pyright: ignore[reportInvalidTypeForm] @@ -175,30 +179,76 @@ def plot_A_dep_3d( ), Q: float | None = None, Q2: float | None = None, - x_lines: float | list[float] | None = None, - colors: list = [], - logA: bool = True, - plot_uncertainty: bool = True, + A_lines: float | list[float] | None = None, + colors: list[str] | str | cycle = [], + offset: float = 0, + labels_Bjx: Literal["lines", "legend", "none"] = "legend", + name: str = "", plot_ratio: bool = False, - pdf_label: Literal["ylabel", "annotate"] = "annotate", - A_label: Literal["legend", "ticks"] = "ticks", - proj_type: Literal["ortho", "persp"] = "persp", - view_init: tuple[float, float] | list[tuple[float, float]] = (15, -75), + pdf_label: Literal["title", "annotate", "none"] | None = "annotate", + plot_legend: bool = True, kwargs_theory: dict[str, Any] | list[dict[str, Any] | None] = {}, - kwargs_uncertainty: dict[str, Any] | list[dict[str, Any] | None] = {}, - kwargs_uncertainty_edges: dict[str, Any] | list[dict[str, Any] | None] = {}, - kwargs_title: dict[str, Any] = {}, - kwargs_notation: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_legend: dict[str, Any] = {}, kwargs_ylabel: dict[str, Any] | list[dict[str, Any] | None] = {}, - kwargs_xlabel: dict[str, Any] = {}, - kwargs_zlabel: dict[str, Any] = {}, - kwargs_legend: dict[str, Any] = {}, - kwargs_xlines: dict[str, Any] | list[dict[str, Any] | None] = {}, - ): + kwargs_title: dict[str, Any] | list[dict[str, Any] | None]= {}, + kwargs_annotate: dict[str, Any] | list[dict[str, Any] | None] = {}, + kwargs_xticks: list[dict[str, Any] | None] = {}, + ) -> None: + """Plot nuclear PDFs in the A-f plane for different values of x. + + Parameters + ---------- + ax : matplotlib.axes.Axes | numpy.ndarray[matplotlib.axes.Axes] + The axes to plot on. + x : float | list[float] + The x values to plot for. + observables : sympy.Basic | numpy.ndarray[sympy.Basic] | list[sympy.Basic] + The observables to plot. + Q : float, optional + The scale at which to plot the PDFs + Q2 : float, optional + The Q^2 scale at which to plot the PDFs. Either Q or Q2 has to be passed. + colors : list, optional + The colors to use for the different x values, by default [], tab color palette is used if == []. + logx : bool, optional + If True, use a logarithmic scale for the x axis, by default True. + title : str | list[str] | None, optional + The title of the plot, by default None. If a list is passed, the titles are set for each subplot. If a single string is passed, it is set for the first subplot. + plot_unc : bool, optional + If True, plot the uncertainties, by default False. + plot_ratio : bool, optional + If True, plot the ratio of the PDFs to the Proon PDF, by default False. + pdf_label : str, optional + The label for the PDF, by default "annotate". If "ylabel", the label is set as the y-axis label. If "annotate", the label is set as an anannotate in the top right corner of the plot. + kwargs_theory : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the plot function for the central PDF, by default {}. + kwargs_legend : dict[str, Any], optional + The keyword arguments to pass to the legend function, by default {}. + kwargs_xlabel : dict[str, Any], optional + The keyword arguments to pass to the xlabel function, by default {}. + kwargs_ylabel : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the ylabel function, by default {}. + kwargs_title : dict[str, Any], optional + The keyword arguments to pass to the title function, by default {}. + kwargs_annotate : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the anannotate function, by default {}. + kwargs_uncertainty : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the fill_between function for the uncertainties, by default {}. + kwargs_uncertainty_edges : dict[str, Any] | list[dict[str, Any] | None], optional + The keyword arguments to pass to the plot function for the uncertainty edges, by default {}. + """ my_sets = self.pdf_sets - x = self.get(A=my_sets["A"][0]).x_values + my_data = {} + + if not isinstance(x, list): + x = [x] + for x_i in x: + if x_i not in self.get(A=my_sets["A"][0]).x_values: + raise ValueError( + f"Chosen x value {x_i} was not used for defining nuclear pdf set. \n Pleas choose x that was used in initialization" + ) if Q is None and Q2 is None: raise ValueError("Please pass either `Q` or `Q2`") @@ -221,200 +271,216 @@ def plot_A_dep_3d( raise ValueError( f"Chosen Q2 value {Q2} was not used for defining nuclear pdf set. \n Please choose Q2 that was used in initialization" ) - - if not isinstance(A, list): - A = [A] - - if 1 not in A and plot_ratio: - raise ValueError( - "Please pass A=1 if you want to plot the ratio to Proton." - ) - if isinstance(observables, np.ndarray): observables = list(observables.flatten()) if not isinstance(observables, list): observables = [observables] + if isinstance(colors, str): + colors = len(x) * [colors] + + elif isinstance(colors, list) and colors != []: + if len(colors) != len(x): + raise ValueError("No. of colors must match no. of x-values") + + for obs in observables: + data_obs = {} + list_x = [] + list_central = [] + list_unc1 = [] + list_unc2 = [] + list_A = [] + for A in list(my_sets["A"]): + + for x_i in self.get(A=A).x_values: + list_x.append(x_i) + list_central.append( + self.get(A=A).get_central(x=x_i, Q=Q, Q2=Q2, observable=obs) + ) + unc1 = self.get(A=A).get_uncertainties( + x=x_i, Q=Q, Q2=Q2, observable=obs + )[0] + unc2 = self.get(A=A).get_uncertainties( + x=x_i, Q=Q, Q2=Q2, observable=obs + )[1] + if math.isnan(unc1): + list_unc1.append( + self.get(A=A).get_central(x=x_i, Q=Q, Q2=Q2, observable=obs) + ) + else: + list_unc1.append(abs(self.get(A=A).get_central(x=x_i, Q=Q, Q2=Q2, observable=obs)-unc1)) + if math.isnan(unc2): + list_unc2.append( + self.get(A=A).get_central(x=x_i, Q=Q, Q2=Q2, observable=obs) + ) + else: + list_unc2.append(abs(-self.get(A=A).get_central(x=x_i, Q=Q, Q2=Q2, observable=obs)+unc2)) + i = 0 + while i < len(self.get(A=A).x_values): + list_A.append(A) + i += 1 + + data_obs["A"] = list_A + data_obs["x"] = list_x + data_obs["central"] = list_central + data_obs["err_mi"] = list_unc1 + data_obs["err_pl"] = list_unc2 + + dataframe_obs = pd.DataFrame(data_obs) + my_data[obs] = dataframe_obs + + # fig, ax = plt.subplots(1, len(observables), figsize=(9 * len(observables), 5)) + if not isinstance(ax, np.ndarray): ax = np.array([ax]) - if colors == []: - cmap = cm.get_cmap("tab10", lut=len(A)) - colors = [cmap(i) for i in range(len(A))] + for m, (obs_m, ax_m) in enumerate(zip(observables, ax.flat)): + ax_m: plt.Axes + - for i, (obs_i, ax_i) in enumerate(zip(observables, ax.flat)): + if colors == []: + colors = cycle(plt.rcParams["axes.prop_cycle"].by_key()["color"]) - for j, (A_j, col_j) in enumerate(zip(A, colors)): + markers=cycle(["x","*","o","^"]) + + + for j, x_j in enumerate(x): + if isinstance(colors, str): + col = colors + elif isinstance(colors, list): + col = colors[j] + else: + col = next(colors) + mk=next(markers) + if j==0 and labels_Bjx=="legend": + kwargs_default = { + "color": col, + "label": f"x={x_j}, {name}", + "linestyle":"", + "capsize":5, + "fmt":mk, + "ecolor":col + } + elif j!=0 and labels_Bjx=="legend": + kwargs_default = { + "color": col, + "label": f"x={x_j}", + "linestyle":"", + "capsize":5, + "fmt":mk, + "ecolor":col + } + elif j==0 and labels_Bjx!="legend": + kwargs_default = { + "color": col, + "label": f"{name}", + "linestyle":"", + "capsize":5, + "fmt":mk, + "ecolor":col + } + else: + kwargs_default = { + "color": col, + "linestyle":"", + "capsize":5, + "fmt":mk, + "ecolor":col + } + + kwargs = update_kwargs( + kwargs_default, + kwargs_theory, + i=j, + ) if not plot_ratio: - z_lower, z_upper = self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2 - ) - else: - z_lower, z_upper = self.get(A=A_j).get_uncertainties( - observable=obs_i, x=x, Q=Q, Q2=Q2 - ) - z_lower = z_lower / self.get(A=1).get_central(observable=obs_i,x=x, Q=Q, Q2=Q2) - z_upper = z_upper / self.get(A=1).get_central(observable=obs_i,x=x, Q=Q, Q2=Q2) - kwargs_default = { - "color": col_j, - "label": f"A={A_j}", - "linewidth": 1.5, - } - if not isinstance(kwargs_theory, list): - kwargs = update_kwargs( - kwargs_default, - kwargs_theory, + ax_m.errorbar( + np.array(range(len(A_lines)))+offset*np.ones(shape=len(A_lines)), + my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["central"], + yerr=( + my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["err_mi"], + my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["err_pl"]), + **kwargs, ) else: - kwargs = update_kwargs( - kwargs_default, - kwargs_theory, - i=j, + ax_m.errorbar( + np.array(range(len(my_data[obs_m].query(f"x=={x_j}and A in {A_lines}")["A"])))+offset*np.ones(shape=len(A_lines)), + np.array(list(my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["central"]))/np.array(my_data[obs_m].query(f"A=={1} & x=={x_j}")["central"]), + yerr=( + np.array(list(my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["err_mi"]))/np.array(my_data[obs_m].query(f"A=={1} & x=={x_j}")["central"]), + np.array(list(my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["err_pl"]))/np.array(my_data[obs_m].query(f"A=={1} & x=={x_j}")["central"])), + **kwargs, ) - if logA: - if plot_ratio: - ax_i.plot( - np.log10(x), - np.log10(len(x) * [A_j]), - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i) - / self.get(A=1).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), - **kwargs, - ) - else: - ax_i.plot( - np.log10(x), - np.log10(len(x) * [A_j]), - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), - **kwargs, - ) - else: - if plot_ratio: - ax_i.plot( - np.log10(x), - len(x) * [A_j], - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i) - / self.get(A=1).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), - **kwargs, - ) - else: - ax_i.plot( - np.log10(x), - len(x) * [A_j], - self.get(A=A_j).get_central(x=x, Q=Q, Q2=Q2, observable=obs_i), - **kwargs, - ) - if plot_uncertainty: - kwargs_uncertainty_default = { - "color": col_j, - "alpha": 0.3, - } - if not isinstance(kwargs_uncertainty, list): - kwargs = update_kwargs( - kwargs_uncertainty_default, - kwargs_uncertainty, - ) - else: - kwargs = update_kwargs( - kwargs_uncertainty_default, - kwargs_uncertainty, - i=j, - ) - vertices = [] - z_lower = np.array(z_lower) - z_upper = np.array(z_upper) - if not logA: - - for xi, ai, zl, zu in zip( - np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper - ): - vertices.append([xi, ai, zl]) - for xi, ai, zl, zu in reversed( - list(zip(np.log10(x), np.ones(len(x)) * A_j, z_lower, z_upper)) - ): - vertices.append([xi, ai, zu]) - - else: - for xi, ai, zl, zu in zip( - np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper - ): - vertices.append([xi, ai, zl]) - - for xi, ai, zl, zu in reversed( - list(zip(np.log10(x), np.ones(len(x)) * np.log10(A_j), z_lower, z_upper)) - ): - vertices.append([xi, ai, zu]) - poly = Poly3DCollection([vertices], **kwargs) - ax_i.add_collection3d(poly) - - kwargs_uncertainty_edges_default = { - "color": col_j, - "alpha": 1, - } - if not isinstance(kwargs_uncertainty_edges, list): - kwargs = update_kwargs( - kwargs_uncertainty_edges_default, - kwargs_uncertainty_edges, - ) - else: - kwargs = update_kwargs( - kwargs_uncertainty_edges_default, - kwargs_uncertainty_edges, - i=j, - ) - if not logA: - ax_i.plot(np.log10(x), len(x) * [A_j], z_upper, **kwargs) - ax_i.plot(np.log10(x), len(x) * [A_j], z_lower, **kwargs) - else: - ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_upper, **kwargs) - ax_i.plot(np.log10(x), len(x) * [np.log10(A_j)], z_lower, **kwargs) - - centrals={} - if x_lines is not None: - if not isinstance(x_lines, list): - x_lines = [x_lines] - for k,x_line in enumerate(x_lines): - if x_line not in x: - raise ValueError( - f"Chosen x value {x_line} was not used for defining nuclear pdf set. \n Please choose x that was used in initialization" - ) - kwargs_xlines_default = { - "color": "black", - "linestyle": "--", - "linewidth": 1.5, + + if A_lines is not None: + if not isinstance(A_lines, list): + A_lines = [A_lines] + + ax_m.set_xticks( + range(len(A_lines)), + labels=[ + f"{elements.element_to_str(A=A_line,long=True)} {A_line}" + for A_line in A_lines + ], + ha="right", + rotation=30, + ) + ax_m.xaxis.set_tick_params(which="minor", size=0) + + ax_m.grid(visible=True) + + + kwargs_ylabel_default = { + "ylabel": f"${util.to_str(obs_m)}$", + } + + if isinstance(kwargs_ylabel, list): + kwargs_y = update_kwargs(kwargs_ylabel_default, kwargs_ylabel, i=m) + else: + kwargs_y = update_kwargs( + kwargs_ylabel_default, + kwargs_ylabel, + ) + + ax_m.set_ylabel(**kwargs_y) + + if labels_Bjx == "lines": + if m == len(ax.flat) - 1: + + kwargs_legend_default = { + "loc": "upper left", + "bbox_to_anchor": (1, 1), + "frameon": False, } - if not isinstance(kwargs_xlines, list): - kwargs = update_kwargs( - kwargs_xlines_default, - kwargs_xlines, - ) - else: - kwargs = update_kwargs( - kwargs_xlines_default, - kwargs_xlines, - i=k, + kwargs_legend = update_kwargs( + kwargs_legend_default, + kwargs_legend, + ) + + ax_m.legend(**kwargs_legend) + + for x_j in x: + if not plot_ratio: + ax_m.annotate(f"x={x_j}",fontsize= 10, + xy= (0, my_data[obs_m].query(f"x=={x_j} and A ==1 ")["central"].iloc[0])) + if labels_Bjx == "legend": + if plot_legend: + kwargs_legend_default = { + "loc": "upper left", + "bbox_to_anchor": (1, 1), + "frameon": False, + } + kwargs_legend = update_kwargs( + kwargs_legend_default, + kwargs_legend, ) - for a in A: - if not plot_ratio: - if x_line not in centrals.keys(): - centrals[x_line] = [self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)] - else: - centrals[x_line].append(self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)) - else: - if x_line not in centrals.keys(): - centrals[x_line] = [self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)/self.get(A=1).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)] - else: - centrals[x_line].append(self.get(A=a).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)/self.get(A=1).get_central(x=x_line, Q=Q, Q2=Q2, observable=obs_i)) - if logA: - ax_i.plot( - np.ones(len(A))*np.log10(x_line), np.log10(A), centrals[x_line],**kwargs) - else: - ax_i.plot( - np.ones(len(A))*np.log10(x_line), A, centrals[x_line],**kwargs) + + ax[-1].legend(**kwargs_legend) if pdf_label == "annotate": - kwargs_notation_default = { + kwargs_annotate_default = { "fontsize": 12, "xy": (0.97, 0.96), "xycoords": "axes fraction", @@ -427,117 +493,17 @@ def plot_A_dep_3d( boxstyle="round, pad=0.2", ), } - if not isinstance(kwargs_notation, list): - kwargs_n = update_kwargs( - kwargs_notation_default, - kwargs_notation, - ) - else: - kwargs_n = update_kwargs( - kwargs_notation_default, - kwargs_notation, - i=i, - ) - ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) - - if pdf_label == "ylabel": - kwargs_ylabel_default = { - "fontsize": 14, - "zlabel": f"${util.to_str(obs_i,Q=Q,Q2=Q2)}$", - } - if not isinstance(kwargs_ylabel, list): - kwargs = update_kwargs( - kwargs_ylabel_default, - kwargs_ylabel, - ) - else: - kwargs = update_kwargs( - kwargs_ylabel_default, - kwargs_ylabel, - i=i, - ) - ax_i.set_zlabel(**kwargs) - - else: - kwargs_notation_default = { - "fontsize": 12, - "xy": (0.47, 0.96), - "xycoords": "axes fraction", - "va": "top", - "ha": "right", - "bbox": dict( - facecolor=(1, 1, 1), - edgecolor=(0.8, 0.8, 0.8), - lw=0.9, - boxstyle="round, pad=0.2", - ), + kwargs_n = update_kwargs(kwargs_annotate_default, kwargs_annotate, i=m) + ax_m.annotate(f"${util.to_str(obs_m, Q=Q,Q2=Q2)}$", **kwargs_n) + + if pdf_label == "title": + kwargs_title_default = { + "y": 1.05, + "loc": "center", + "label": f"${util.to_str(obs_m,Q=Q,Q2=Q2)}$", } - kwargs_n = update_kwargs(kwargs_notation_default, kwargs_notation, i=i) - - ax_i.annotate(f"${util.to_str(obs_i, Q=Q,Q2=Q2)}$", **kwargs_n) - ax_i.xaxis.set_major_formatter(mticker.FuncFormatter(log_tick_formatter)) - if A_label == "ticks": - if logA: - ax_i.set_yticks(np.log10(A),A) + if isinstance(kwargs_title, list): + kwargs_t = update_kwargs(kwargs_title_default, kwargs_title, i=m) else: - ax_i.set_yticks(A,A) - kwargs_zlabel_default = { - "fontsize": 14, - "ylabel": f"$A$", - - } - kwargs = update_kwargs( - kwargs_zlabel_default, - kwargs_zlabel, - ) - - ax_i.set_ylabel(**kwargs) - ax_i.set_xlim(np.log10(x[0])*0.98) - else: - ax_i.set_yticks([]) - if i==len(observables)-1: - kwargs_legend_default = { - "fontsize": 12, - "bbox_to_anchor": (0.95, 0.95), - "frameon": False, - } - kwargs = update_kwargs( - kwargs_legend_default, - kwargs_legend, - ) - ax_i.legend() - kwargs_zlabel_default = { - "fontsize": 14, - "ylabel": f"$A$", - "labelpad":-10, - "linespacing": -4 - } - kwargs = update_kwargs( - kwargs_zlabel_default, - kwargs_zlabel, - ) - ax_i.set_ylabel(**kwargs) - ax_i.xaxis.pane.fill=False - ax_i.yaxis.pane.fill=False - ax_i.zaxis.pane.fill=False - ax_i.xaxis.pane.set_edgecolor("w") - ax_i.yaxis.pane.set_edgecolor("w") - ax_i.zaxis.pane.set_edgecolor("w") - - ax_i.zaxis._axinfo["juggled"]=(1,2,0) - - kwargs_xlabel_default = { - "fontsize": 14, - "xlabel": f"$x$", - - } - kwargs = update_kwargs( - kwargs_xlabel_default, - kwargs_xlabel, - ) - ax_i.set_xlabel(**kwargs) - #, np.log10(x[-1])) - ax_i.set_zlim(ax_i.get_zlim()[1]*0.02) - #ax_i.yaxis._axinfo["grid"]["linewidth"] = 0 - ax_i.set_proj_type(proj_type) - ax_i.view_init(*view_init[i] if isinstance(view_init, list) else view_init) \ No newline at end of file + kwargs_t = update_kwargs(kwargs_title_default, kwargs_title) + ax_m.set_title(**kwargs_t) From c59c080b39f55f17d3cd073556998c4408fec64f Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Wed, 7 May 2025 09:35:19 +0200 Subject: [PATCH 5/6] Add possibility to label different x with an annotation --- pdfplotter/pdf_set_nuclear.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index 2e77b80..f9243d6 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -182,7 +182,7 @@ def plot_A_dep_2d_scatter( A_lines: float | list[float] | None = None, colors: list[str] | str | cycle = [], offset: float = 0, - labels_Bjx: Literal["lines", "legend", "none"] = "legend", + labels_Bjx: Literal["annotate", "legend", "none"] = "legend", name: str = "", plot_ratio: bool = False, pdf_label: Literal["title", "annotate", "none"] | None = "annotate", @@ -446,7 +446,7 @@ def plot_A_dep_2d_scatter( ax_m.set_ylabel(**kwargs_y) - if labels_Bjx == "lines": + if labels_Bjx == "annotate": if m == len(ax.flat) - 1: kwargs_legend_default = { @@ -464,7 +464,8 @@ def plot_A_dep_2d_scatter( for x_j in x: if not plot_ratio: ax_m.annotate(f"x={x_j}",fontsize= 10, - xy= (0, my_data[obs_m].query(f"x=={x_j} and A ==1 ")["central"].iloc[0])) + xy= (0, my_data[obs_m].query(f"x=={x_j} and A ==1 ")["central"].iloc[0]+max(my_data[obs_m].query(f"A ==1 ")["central"])*0.025)) + if labels_Bjx == "legend": if plot_legend: kwargs_legend_default = { From 381b28d4a5e3dc1d6399e06f5b333d62598ddd27 Mon Sep 17 00:00:00 2001 From: Katrin-Greve Date: Mon, 12 May 2025 10:46:46 +0200 Subject: [PATCH 6/6] Format code, update description --- pdfplotter/pdf_set_nuclear.py | 213 ++++++++++++++++++++-------------- 1 file changed, 129 insertions(+), 84 deletions(-) diff --git a/pdfplotter/pdf_set_nuclear.py b/pdfplotter/pdf_set_nuclear.py index f9243d6..8b5a20b 100644 --- a/pdfplotter/pdf_set_nuclear.py +++ b/pdfplotter/pdf_set_nuclear.py @@ -4,7 +4,7 @@ from math import log from typing import Sequence, Literal, Any -#from mpl_toolkits.mplot3d import Axes3D +# from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d.art3d import Poly3DCollection import matplotlib.ticker as mticker import matplotlib.cm as cm @@ -190,9 +190,8 @@ def plot_A_dep_2d_scatter( kwargs_theory: dict[str, Any] | list[dict[str, Any] | None] = {}, kwargs_legend: dict[str, Any] = {}, kwargs_ylabel: dict[str, Any] | list[dict[str, Any] | None] = {}, - kwargs_title: dict[str, Any] | list[dict[str, Any] | None]= {}, + kwargs_title: dict[str, Any] | list[dict[str, Any] | None] = {}, kwargs_annotate: dict[str, Any] | list[dict[str, Any] | None] = {}, - kwargs_xticks: list[dict[str, Any] | None] = {}, ) -> None: """Plot nuclear PDFs in the A-f plane for different values of x. @@ -208,34 +207,28 @@ def plot_A_dep_2d_scatter( The scale at which to plot the PDFs Q2 : float, optional The Q^2 scale at which to plot the PDFs. Either Q or Q2 has to be passed. + A_lines: + A values to evaluate the PDFs at. x-ticks are set to those A values colors : list, optional The colors to use for the different x values, by default [], tab color palette is used if == []. - logx : bool, optional - If True, use a logarithmic scale for the x axis, by default True. - title : str | list[str] | None, optional - The title of the plot, by default None. If a list is passed, the titles are set for each subplot. If a single string is passed, it is set for the first subplot. - plot_unc : bool, optional - If True, plot the uncertainties, by default False. + offset: + The offset in A, where the errorbar is plottet (useful, when several PDF errorbars are plotted in the same panel) + labels_BjX: + Whether the different values for x are shown via a legend or an annotation plot_ratio : bool, optional - If True, plot the ratio of the PDFs to the Proon PDF, by default False. + If True, plot the ratio of the PDFs to the Proton PDF, by default False. (only recomended if only one PDF set is plotted in the panel) pdf_label : str, optional - The label for the PDF, by default "annotate". If "ylabel", the label is set as the y-axis label. If "annotate", the label is set as an anannotate in the top right corner of the plot. + The label for the PDF, by default "annotate". If "title", it is used in ax.set_title(). If "annotate", the label is set as an anannotate. kwargs_theory : dict[str, Any] | list[dict[str, Any] | None], optional The keyword arguments to pass to the plot function for the central PDF, by default {}. kwargs_legend : dict[str, Any], optional The keyword arguments to pass to the legend function, by default {}. - kwargs_xlabel : dict[str, Any], optional - The keyword arguments to pass to the xlabel function, by default {}. kwargs_ylabel : dict[str, Any] | list[dict[str, Any] | None], optional The keyword arguments to pass to the ylabel function, by default {}. kwargs_title : dict[str, Any], optional The keyword arguments to pass to the title function, by default {}. kwargs_annotate : dict[str, Any] | list[dict[str, Any] | None], optional The keyword arguments to pass to the anannotate function, by default {}. - kwargs_uncertainty : dict[str, Any] | list[dict[str, Any] | None], optional - The keyword arguments to pass to the fill_between function for the uncertainties, by default {}. - kwargs_uncertainty_edges : dict[str, Any] | list[dict[str, Any] | None], optional - The keyword arguments to pass to the plot function for the uncertainty edges, by default {}. """ my_sets = self.pdf_sets @@ -256,21 +249,14 @@ def plot_A_dep_2d_scatter( elif Q is not None and Q2 is not None: raise ValueError("Please pass either `Q` or `Q2`, not both") - elif Q is not None: + elif Q is not None and Q2 is None: if Q not in self.get(A=my_sets["A"][0]).Q_values and Q not in np.sqrt( np.array(self.get(A=my_sets["A"][0]).Q2_values) ): raise ValueError( f"Chosen Q value {Q} was not used for defining nuclear pdf set. \n Please choose Q that was used in initialization" ) - else: - if ( - Q2 not in self.get(A=my_sets["A"][0]).Q2_values - and Q2 not in np.array(self.get(A=my_sets["A"][0]).Q_values) ** 2 - ): - raise ValueError( - f"Chosen Q2 value {Q2} was not used for defining nuclear pdf set. \n Please choose Q2 that was used in initialization" - ) + if isinstance(observables, np.ndarray): observables = list(observables.flatten()) @@ -309,13 +295,27 @@ def plot_A_dep_2d_scatter( self.get(A=A).get_central(x=x_i, Q=Q, Q2=Q2, observable=obs) ) else: - list_unc1.append(abs(self.get(A=A).get_central(x=x_i, Q=Q, Q2=Q2, observable=obs)-unc1)) + list_unc1.append( + abs( + self.get(A=A).get_central( + x=x_i, Q=Q, Q2=Q2, observable=obs + ) + - unc1 + ) + ) if math.isnan(unc2): list_unc2.append( self.get(A=A).get_central(x=x_i, Q=Q, Q2=Q2, observable=obs) ) else: - list_unc2.append(abs(-self.get(A=A).get_central(x=x_i, Q=Q, Q2=Q2, observable=obs)+unc2)) + list_unc2.append( + abs( + -self.get(A=A).get_central( + x=x_i, Q=Q, Q2=Q2, observable=obs + ) + + unc2 + ) + ) i = 0 while i < len(self.get(A=A).x_values): list_A.append(A) @@ -330,21 +330,17 @@ def plot_A_dep_2d_scatter( dataframe_obs = pd.DataFrame(data_obs) my_data[obs] = dataframe_obs - # fig, ax = plt.subplots(1, len(observables), figsize=(9 * len(observables), 5)) - if not isinstance(ax, np.ndarray): ax = np.array([ax]) for m, (obs_m, ax_m) in enumerate(zip(observables, ax.flat)): ax_m: plt.Axes - if colors == []: colors = cycle(plt.rcParams["axes.prop_cycle"].by_key()["color"]) - markers=cycle(["x","*","o","^"]) + markers = cycle(["x", "*", "o", "^"]) - for j, x_j in enumerate(x): if isinstance(colors, str): col = colors @@ -352,43 +348,43 @@ def plot_A_dep_2d_scatter( col = colors[j] else: col = next(colors) - mk=next(markers) - if j==0 and labels_Bjx=="legend": + mk = next(markers) + if j == 0 and labels_Bjx == "legend": kwargs_default = { "color": col, "label": f"x={x_j}, {name}", - "linestyle":"", - "capsize":5, - "fmt":mk, - "ecolor":col + "linestyle": "", + "capsize": 5, + "fmt": mk, + "ecolor": col, } - elif j!=0 and labels_Bjx=="legend": + elif j != 0 and labels_Bjx == "legend": kwargs_default = { "color": col, "label": f"x={x_j}", - "linestyle":"", - "capsize":5, - "fmt":mk, - "ecolor":col + "linestyle": "", + "capsize": 5, + "fmt": mk, + "ecolor": col, } - elif j==0 and labels_Bjx!="legend": + elif j == 0 and labels_Bjx != "legend": kwargs_default = { "color": col, "label": f"{name}", - "linestyle":"", - "capsize":5, - "fmt":mk, - "ecolor":col + "linestyle": "", + "capsize": 5, + "fmt": mk, + "ecolor": col, } else: kwargs_default = { "color": col, - "linestyle":"", - "capsize":5, - "fmt":mk, - "ecolor":col - } - + "linestyle": "", + "capsize": 5, + "fmt": mk, + "ecolor": col, + } + kwargs = update_kwargs( kwargs_default, kwargs_theory, @@ -396,24 +392,66 @@ def plot_A_dep_2d_scatter( ) if not plot_ratio: ax_m.errorbar( - np.array(range(len(A_lines)))+offset*np.ones(shape=len(A_lines)), + np.array(range(len(A_lines))) + + offset * np.ones(shape=len(A_lines)), my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["central"], yerr=( - my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["err_mi"], - my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["err_pl"]), + my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")[ + "err_mi" + ], + my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")[ + "err_pl" + ], + ), **kwargs, ) else: ax_m.errorbar( - np.array(range(len(my_data[obs_m].query(f"x=={x_j}and A in {A_lines}")["A"])))+offset*np.ones(shape=len(A_lines)), - np.array(list(my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["central"]))/np.array(my_data[obs_m].query(f"A=={1} & x=={x_j}")["central"]), + np.array( + range( + len( + my_data[obs_m].query(f"x=={x_j}and A in {A_lines}")[ + "A" + ] + ) + ) + ) + + offset * np.ones(shape=len(A_lines)), + np.array( + list( + my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")[ + "central" + ] + ) + ) + / np.array( + my_data[obs_m].query(f"A=={1} & x=={x_j}")["central"] + ), yerr=( - np.array(list(my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["err_mi"]))/np.array(my_data[obs_m].query(f"A=={1} & x=={x_j}")["central"]), - np.array(list(my_data[obs_m].query(f"x=={x_j} and A in {A_lines}")["err_pl"]))/np.array(my_data[obs_m].query(f"A=={1} & x=={x_j}")["central"])), + np.array( + list( + my_data[obs_m].query( + f"x=={x_j} and A in {A_lines}" + )["err_mi"] + ) + ) + / np.array( + my_data[obs_m].query(f"A=={1} & x=={x_j}")["central"] + ), + np.array( + list( + my_data[obs_m].query( + f"x=={x_j} and A in {A_lines}" + )["err_pl"] + ) + ) + / np.array( + my_data[obs_m].query(f"A=={1} & x=={x_j}")["central"] + ), + ), **kwargs, ) - if A_lines is not None: if not isinstance(A_lines, list): A_lines = [A_lines] @@ -431,11 +469,10 @@ def plot_A_dep_2d_scatter( ax_m.grid(visible=True) - kwargs_ylabel_default = { "ylabel": f"${util.to_str(obs_m)}$", } - + if isinstance(kwargs_ylabel, list): kwargs_y = update_kwargs(kwargs_ylabel_default, kwargs_ylabel, i=m) else: @@ -450,33 +487,41 @@ def plot_A_dep_2d_scatter( if m == len(ax.flat) - 1: kwargs_legend_default = { - "loc": "upper left", - "bbox_to_anchor": (1, 1), - "frameon": False, + "loc": "upper left", + "bbox_to_anchor": (1, 1), + "frameon": False, } kwargs_legend = update_kwargs( - kwargs_legend_default, - kwargs_legend, + kwargs_legend_default, + kwargs_legend, ) ax_m.legend(**kwargs_legend) - - for x_j in x: + + for x_j in x: if not plot_ratio: - ax_m.annotate(f"x={x_j}",fontsize= 10, - xy= (0, my_data[obs_m].query(f"x=={x_j} and A ==1 ")["central"].iloc[0]+max(my_data[obs_m].query(f"A ==1 ")["central"])*0.025)) + ax_m.annotate( + f"x={x_j}", + fontsize=10, + xy=( + 0 + offset, + my_data[obs_m] + .query(f"x=={x_j} and A ==1 ")["central"] + .iloc[0], + ), + ) if labels_Bjx == "legend": if plot_legend: kwargs_legend_default = { - "loc": "upper left", - "bbox_to_anchor": (1, 1), - "frameon": False, - } + "loc": "upper left", + "bbox_to_anchor": (1, 1), + "frameon": False, + } kwargs_legend = update_kwargs( - kwargs_legend_default, - kwargs_legend, - ) + kwargs_legend_default, + kwargs_legend, + ) ax[-1].legend(**kwargs_legend) @@ -499,9 +544,9 @@ def plot_A_dep_2d_scatter( if pdf_label == "title": kwargs_title_default = { - "y": 1.05, - "loc": "center", - "label": f"${util.to_str(obs_m,Q=Q,Q2=Q2)}$", + "y": 1.05, + "loc": "center", + "label": f"${util.to_str(obs_m,Q=Q,Q2=Q2)}$", } if isinstance(kwargs_title, list): kwargs_t = update_kwargs(kwargs_title_default, kwargs_title, i=m)