Mix bands and replicas

validphys2/examples/mixbandreplicas.yaml

@@ -0,0 +1,23 @@
+meta:
+  title: Example for mixing bands and replicas
+  author: The Creator
+  keywords: [replicas, bands, mix]
+
+Q: 1.7
+pdfs:
+- {id: NNPDF40_nnlo_as_01180, label: NNLO}
+- {id: NNPDF40_nlo_as_01180, label: NLO}
+
+mixband_as_replicas: [2,]
+basis: flavour
+xscale: log
+npoints: 20
+template_text: |
+  Plots
+  --------------
+
+  {@plot_pdfs_mixed@}
+  {@plot_pdfs_mixed_kinetic_energy@}
+actions_:
+- report(main=True)
+

validphys2/src/validphys/checks.py

@@ -230,6 +230,45 @@ def check_pdfs_noband(pdfs, pdfs_noband):
     return {'pdfs_noband': pdfs_noband_combined}
 
 
+@make_argcheck
+def check_mixband_as_replicas(pdfs, mixband_as_replicas):
+    """Same as check_pdfs_noband, but for the mixband_as_replicas key.
+    Allows mixband_as_replicas to be specified as a list of PDF IDs or a list of
+    PDF indexes (starting from one)."""
+
+    msg = ("mixband_as_replicas should be a list of PDF IDs (strings) or a list of "
+           "PDF indexes (integers, starting from one)")
+    msg_range = ("At least one of the choices in mixband_as_replicas indexes is out of range. "
+                 "Note that pdf_noband indexing starts at 1, not 0.")
+
+    if mixband_as_replicas is None:
+        return {'mixband_as_replicas': []}
+
+    names = [pdf.name for pdf in pdfs]
+    # A list to which PDF IDs can be added, when the PDF is specified by either
+    # its PDF ID (i.e. a string) or an index (i.e. an int)
+    mixband_as_replicas_combined = []
+
+    for pdf_noband in mixband_as_replicas:
+        if isinstance(pdf_noband, int):
+            if not pdf_noband <= len(names) or pdf_noband < 0:
+                raise CheckError(msg_range)
+            # Convert PDF index to list index (i.e. starting from zero)
+            pdf_noband -= 1
+            mixband_as_replicas_combined.append(pdfs[pdf_noband])
+
+        elif isinstance(pdf_noband, str):
+            try:
+                pdf_index = names.index(pdf_noband)
+                mixband_as_replicas_combined.append(pdfs[pdf_index])
+            except ValueError:
+                raise CheckError(msg, pdf_noband, alternatives=names)
+
+        else:
+            raise CheckError(msg)
+
+    return {'mixband_as_replicas': mixband_as_replicas_combined}
+
 def _check_list_different(l, name):
     strs = [str(item) for item in l]
     if not len(set(strs))==len(l):

validphys2/src/validphys/pdfplots.py

@@ -25,7 +25,7 @@
 from validphys.gridvalues import LUMI_CHANNELS
 from validphys.utils import scale_from_grid
 from validphys.checks import check_pdf_normalize_to, check_scale, check_have_two_pdfs
-from validphys.checks import check_pdfs_noband
+from validphys.checks import check_pdfs_noband, check_mixband_as_replicas
 
 log = logging.getLogger(__name__)
 
@@ -55,7 +55,6 @@ def __init__(self, pdfs, xplotting_grids, xscale, normalize_to, ymin, ymax):
     def setup_flavour(self, flstate):
         pass
 
-
     def normalize(self):
         normalize_to = self.normalize_to
         if normalize_to is not None:
@@ -218,32 +217,6 @@ def plot_pdfreplicas(pdfs, xplotting_grids, xscale:(str,type(None))=None,
     yield from ReplicaPDFPlotter(pdfs=pdfs, xplotting_grids=xplotting_grids,
                                  xscale=xscale, normalize_to=normalize_to, ymin=ymin, ymax=ymax)
 
-
-@figuregen
-@check_pdf_normalize_to
-@check_scale('xscale', allow_none=True)
-@_warn_any_pdf_not_montecarlo
-def plot_pdfreplicas_kinetic_energy(
-    pdfs,
-    kinetic_xplotting_grids,
-    xscale: (str, type(None)) = None,
-    normalize_to: (int, str, type(None)) = None,
-    ymin=None,
-    ymax=None,
-):
-    """Plot the kinetic energy of the replicas of the specified PDFs.
-    Otherise it works the same as ``plot_pdfs_kinetic_energy``.
-    """
-    return plot_pdfreplicas(
-        pdfs,
-        kinetic_xplotting_grids,
-        xscale=xscale,
-        normalize_to=normalize_to,
-        ymin=ymin,
-        ymax=ymax,
-    )
-
-
 class UncertaintyPDFPlotter(PDFPlotter):
 
     def get_ylabel(self, parton_name):
@@ -444,10 +417,8 @@ def draw(self, pdf, grid, flstate):
         handles = flstate.handles
         # Take only the flavours we are interested in
         stats = grid.select_flavour(flstate.flindex).grid_values
-        pcycler = ax._get_lines.prop_cycler
-        #This is ugly but can't think of anything better
 
-        next_prop = next(pcycler)
+        next_prop = next(ax._get_lines.prop_cycler)
         cv = stats.central_value()
         xgrid = grid.xgrid
         #Ignore spurious normalization warnings
@@ -556,6 +527,50 @@ def plot_pdfs(
     )
 
 
+@figuregen
+@check_pdf_normalize_to
+@check_mixband_as_replicas
+@check_pdfs_noband
+@check_scale("xscale", allow_none=True)
+def plot_pdfs_mixed(
+    pdfs,
+    xplotting_grids,
+    xscale: (str, type(None)) = None,
+    normalize_to: (int, str, type(None)) = None,
+    ymin=None,
+    ymax=None,
+    pdfs_noband: (list, type(None)) = None,
+    show_mc_errors: bool = True,
+    legend_stat_labels: bool = True,
+    mixband_as_replicas: (list, type(None)) = None,
+):
+    """This function is similar to plot_pdfs, except instead of only plotting
+    the central value and the uncertainty of the PDFs, those PDFs indicated by
+    mixband_as_replicas will be plotted as replicas without the central value.
+
+    Inputs are the same as plot_pdfs, with the exeption of mixband_as_replicas,
+    which only exists here.
+
+    mixband_as_replicas: A list of PDFs to plot as replicas, i.e. the
+    central values and replicas of these PDFs will be plotted. The list can be
+    formed of strings, corresponding to PDF IDs, integers (starting from one),
+    corresponding to the index of the PDF in the list of PDFs, or a mixture
+    of both.
+    """
+    yield from MixBandPDFPlotter(
+        pdfs,
+        xplotting_grids,
+        xscale,
+        normalize_to,
+        ymin,
+        ymax,
+        pdfs_noband=pdfs_noband,
+        show_mc_errors=show_mc_errors,
+        legend_stat_labels=legend_stat_labels,
+        mixband_as_replicas=mixband_as_replicas,
+    )
+
+
 @figuregen
 @check_pdf_normalize_to
 @check_pdfs_noband
@@ -587,6 +602,65 @@ def plot_pdfs_kinetic_energy(
     )
 
 
+@figuregen
+@check_pdf_normalize_to
+@check_scale('xscale', allow_none=True)
+@_warn_any_pdf_not_montecarlo
+def plot_pdfreplicas_kinetic_energy(
+    pdfs,
+    kinetic_xplotting_grids,
+    xscale: (str, type(None)) = None,
+    normalize_to: (int, str, type(None)) = None,
+    ymin=None,
+    ymax=None,
+):
+    """Plot the kinetic energy of the replicas of the specified PDFs.
+    Otherise it works the same as ``plot_pdfs_kinetic_energy``.
+    """
+    return plot_pdfreplicas(
+        pdfs,
+        kinetic_xplotting_grids,
+        xscale=xscale,
+        normalize_to=normalize_to,
+        ymin=ymin,
+        ymax=ymax,
+    )
+
+
+@figuregen
+@check_pdf_normalize_to
+@check_pdfs_noband
+@check_mixband_as_replicas
+@check_scale("xscale", allow_none=True)
+def plot_pdfs_mixed_kinetic_energy(
+    pdfs,
+    kinetic_xplotting_grids,
+    xscale: (str, type(None)) = None,
+    normalize_to: (int, str, type(None)) = None,
+    ymin=None,
+    ymax=None,
+    pdfs_noband: (list, type(None)) = None,
+    show_mc_errors: bool = True,
+    legend_stat_labels: bool = True,
+    mixband_as_replicas: (list, type(None)) = None,
+):
+    """Mixed band and replica plotting of the "kinetic energy" of the PDF as a
+    function of x for a given value of Q.
+    """
+    return plot_pdfs_mixed(
+        pdfs,
+        kinetic_xplotting_grids,
+        xscale=xscale,
+        normalize_to=normalize_to,
+        ymin=ymin,
+        ymax=ymax,
+        pdfs_noband=pdfs_noband,
+        show_mc_errors=show_mc_errors,
+        legend_stat_labels=legend_stat_labels,
+        mixband_as_replicas=mixband_as_replicas,
+    )
+
+
 class FlavoursPlotter(AllFlavoursPlotter, BandPDFPlotter):
     def get_title(self, parton_name):
         return f'{self.pdfs[0]} Q={self.Q : .1f} GeV'
@@ -722,7 +796,7 @@ def plot_lumi1d(
             f"${LUMI_CHANNELS[lumi_channel]}$ luminosity\n"
             f"$\\sqrt{{s}}={format_number(sqrts/1000)}$ TeV   "
             f"$\\|y|<{format_number(y_cut)}$"
-        )      
+        )
 
     return fig
 
@@ -734,7 +808,7 @@ def plot_lumi1d_uncertainties(
     pdfs_lumis,
     lumi_channel,
     sqrts: numbers.Real,
-    y_cut: (numbers.Real, type(None)) = None,    
+    y_cut: (numbers.Real, type(None)) = None,
     normalize_to=None,
     ymin: (numbers.Real, type(None)) = None,
     ymax: (numbers.Real, type(None)) = None,
@@ -781,7 +855,7 @@ def plot_lumi1d_uncertainties(
     else:
         ax.set_title(
             f"${LUMI_CHANNELS[lumi_channel]}$ luminosity uncertainty\n"
-            f"$\\sqrt{{s}}={format_number(sqrts/1000)}$ TeV   "    
+            f"$\\sqrt{{s}}={format_number(sqrts/1000)}$ TeV   "
             f"$\\|y|<{format_number(y_cut)}$"
         )
     ax.set_ylim(ymin, ymax)
@@ -807,7 +881,7 @@ def plot_lumi1d_replicas(
     scale="log",
 ):
     """This function is similar to `plot_lumi1d`, but instead of plotting
-    the standard deviation and 68% c.i. it plots the luminosities for 
+    the standard deviation and 68% c.i. it plots the luminosities for
     individual replicas.
 
     Plot PDF replica luminosities at a given center of mass energy.
@@ -869,7 +943,7 @@ def plot_lumi1d_replicas(
             f"${LUMI_CHANNELS[lumi_channel]}$ luminosity\n"
             f"$\\sqrt{{s}}={format_number(sqrts/1000)}$ TeV   "
             f"$\\|y|<{format_number(y_cut)}$"
-        )      
+        )
 
     return fig
 
@@ -1044,3 +1118,32 @@ def plot_lumi2d_uncertainty(pdf, lumi_channel, lumigrid2d, sqrts:numbers.Real):
     ax.grid(False)
 
     return fig
+
+
+class MixBandPDFPlotter(BandPDFPlotter):
+    """Special wrapper class to plot, in the same figure, PDF bands and PDF replicas
+    depending on the type of PDF.
+    Practical use: plot together the PDF central values with the NNPDF bands
+    """
+    def __init__(self, *args, mixband_as_replicas, **kwargs):
+        self.mixband_as_replicas = mixband_as_replicas
+        super().__init__(*args, **kwargs)
+
+    def draw(self, pdf, grid, flstate):
+        if pdf in self.mixband_as_replicas:
+            labels = flstate.labels
+            handles = flstate.handles
+            ax = flstate.ax
+            next_prop = next(ax._get_lines.prop_cycler)
+            color = next_prop['color']
+            stats = grid.select_flavour(flstate.flindex).grid_values
+            gv = stats.data
+            ax.plot(grid.xgrid, gv.T, alpha=0.2, linewidth=0.5,
+                    color=color, zorder=1)
+            cv_line = ax.plot(grid.xgrid[0:1], stats.central_value()[0:1], 
+                    color=color, linewidth=2)
+            handle = cv_line[0]
+            labels.append(pdf.label)
+            handles.append(handle)
+            return gv
+        return super().draw(pdf, grid, flstate)