DQMServices/Components/scripts/dqm-plot - master - githubmirror/cmssw

DQMServices/Components/scripts/dqm-plot
2 029 строк80 KB
Bruno Alves
Remove duplicated warning message.
20 июл 2026, 21:17
20 июл 2026, 21:173ad94ac
#!/usr/bin/env python3

"""
Simple DQM comparison plotter using mplhep with CMS styling.

Usage: dqm-plot -s "DQMData/Run 1/HLT/Run summary/Muon/*" [options] $DQM_FILES
See all available options with `dqm-plot -h`.
"""

import ROOT
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.legend as mlegend
import mplhep as hep
import argparse
import os
import re
import random
import string
import multiprocessing
import atexit
from functools import partial
import warnings
import urllib.request
import urllib.error
from pathlib import Path
from matplotlib import colors as _mcolors

# Suppress numpy warnings about subnormal float values
warnings.filterwarnings("ignore", category=UserWarning, module="numpy")

# Set CMS style globally
plt.style.use(hep.style.CMS)

# Use matplotlib's mathtext
plt.rcParams["mathtext.default"] = "regular"


# Per-process cache of open ROOT files. A worker process handles many plot
# tasks, so caching the open TFile keeps each input file opened at most once
# per process instead of once per histogram. Histograms handed back to callers
# are detached clones (SetDirectory(0)), so the cached files can stay open.
_OPEN_ROOT_FILES = {}


def _get_cached_root_file(file_path):
    """Return an open TFile for file_path, reusing a cached handle if present."""
    root_file = _OPEN_ROOT_FILES.get(file_path)
    if not root_file or root_file.IsZombie() or not root_file.IsOpen():
        root_file = ROOT.TFile.Open(file_path, "READ")
        _OPEN_ROOT_FILES[file_path] = root_file
    return root_file


@atexit.register
def _close_cached_root_files():
    """Close any ROOT files cached in this process on interpreter shutdown."""
    for root_file in _OPEN_ROOT_FILES.values():
        try:
            if root_file and root_file.IsOpen():
                root_file.Close()
        except Exception:
            pass
    _OPEN_ROOT_FILES.clear()


class DQMPlotter:
    def __init__(self, figsize=(12, 9), ratio_height=0.3, processors=None, debug=False):
        """
        Initialize the plotter.

        Args:
            figsize: Figure size (width, height) in inches
            ratio_height: Fraction of figure height for ratio plot
            processors: Number of processors to use for multiprocessing (None for auto-detect)
        """
        self.figsize = figsize
        self.ratio_height = ratio_height
        # Color palette from https://cms-analysis.docs.cern.ch/guidelines/plotting/colors/
        self.colors = [
            "#3f8fda",
            "#ffa90e",
            "#bd1f01",
            "#94a4a2",
            "#832db6",
            "#a96b59",
            "#e76300",
            "#b9ac70",
            "#717581",
            "#92dadd",
        ]
        self.markers = [
            "o",  # circle
            "s",  # square
            "^",  # triangle up
            "D",  # diamond
            "v",  # triangle down
            "p",  # pentagon
            "*",  # star
            "h",  # hexagon
            "<",  # triangle left
            ">",  # triangle right
        ]
        self.processors = (
            processors if processors is not None else len(os.sched_getaffinity(0))
        )
        self.debug = debug

    def _extend_color_palette(self, needed: int):
        """Ensure self.colors has at least 'needed' distinct entries."""
        if needed <= len(self.colors):
            return

        extra_needed = needed - len(self.colors)
        new_colors = []

        cmap = plt.colormaps["hsv"]
        for i in range(max(extra_needed, 1)):
            rgba = cmap(i / max(extra_needed, 1))
            hexcol = _mcolors.to_hex(rgba, keep_alpha=False)
            if hexcol not in self.colors and hexcol not in new_colors:
                new_colors.append(hexcol)
            if len(new_colors) >= extra_needed:
                break

        self.colors.extend(new_colors)

    def _clamp_yerrors(self, contents, yerrors):
        """
        Clamp yerrors so content +/- yerr stays within [0, 1].
        Returns a (2, N) numpy array of asymmetric error *distances* (lower, upper)
        as expected by matplotlib's `yerr` argument for asymmetric error bars.
        """
        assert len(contents) == len(yerrors)
        lower = [cont - max(cont - yerr, 0.0) for cont, yerr in zip(contents, yerrors)]
        upper = [min(cont + yerr, 1.0) - cont for cont, yerr in zip(contents, yerrors)]
        return np.array([lower, upper])
    
    def divide_histograms(self, histos1, histos2):
        results = []
        for h1, h2 in zip(histos1, histos2):
            if h1.GetNbinsX() != h2.GetNbinsX():
                raise ValueError("Histogram binning mismatch.")

            rndstr = ''.join(random.choices(string.ascii_letters + string.digits, k=10))
            if isinstance(h1, ROOT.TH2) and isinstance(h2, ROOT.TH2):
                result = h1.ProjectionX('ProjectionX1_' + rndstr)
                tmp = h2.ProjectionX('ProjectionX2_' + rndstr)
            elif isinstance(h1, ROOT.TH1) and isinstance(h2, ROOT.TH1):
                result = h1
                tmp = h2
            tmp.Sumw2()
            result.Sumw2()
            result.Divide(tmp)
            results.append(result)
     
        return results

    def operate_histograms(self, histos1, histos2, operation_type):
        if operation_type == "division":
            return self.divide_histograms(histos1, histos2)
        raise ValueError(f"Unsupported operation '{operation_type}'")

    def _rebin_histogram(self, hist, rebin):
        """
        Rebin a ROOT histogram, leaving the original untouched.

        Args:
            hist: ROOT histogram
            rebin: Either a single int/float, interpreted as the number of
                existing bins to merge together (ROOT's "ngroup" rebinning,
                e.g. rebin=2 merges every 2 bins into 1), or an array-like of
                new bin edges to rebin into variable-width bins.

        Returns:
            A new, independent ROOT histogram with the requested binning.
        """
        rndstr = ''.join(random.choices(string.ascii_letters + string.digits, k=10))
        hist_clone = hist.Clone(f"{hist.GetName()}_rebin_{rndstr}")
        hist_clone.SetDirectory(0)

        if isinstance(rebin, (int, np.integer)):
            if rebin < 1:
                raise ValueError(f"rebin group size must be >= 1, got {rebin}")
            hist_clone.Rebin(int(rebin))
            return hist_clone

        # Otherwise treat as an array-like of new bin edges (variable rebinning)
        edges = np.asarray(rebin, dtype=float)
        if edges.ndim != 1 or len(edges) < 2:
            raise ValueError(
                "rebin must be either a single int (number of bins to merge) "
                "or an array-like of at least 2 bin edges."
            )
        if np.any(np.diff(edges) <= 0):
            raise ValueError("rebin bin edges must be strictly increasing.")

        n_new_bins = len(edges) - 1
        rebinned = hist_clone.Rebin(n_new_bins, f"{hist.GetName()}_rebin_var_{rndstr}", edges)
        if not rebinned:
            raise ValueError(
                "Rebin failed - check that the provided bin edges are compatible "
                "with the histogram's original binning."
            )
        rebinned.SetDirectory(0)
        return rebinned

    def root_to_numpy(self, hist, rebin=None):
        """
        Convert ROOT histogram to numpy arrays.

        Args:
            hist: ROOT histogram
            rebin: Optional rebinning specification. Either a single int
                (number of existing bins to merge together) or an array-like
                of new bin edges for variable-width rebinning.

        Returns:
            tuple: (bin_centers, bin_contents, bin_errors, bin_edges, bin_labels, has_labels)
        """
        if rebin is not None:
            hist = self._rebin_histogram(hist, rebin)

        n_bins = hist.GetNbinsX()
        bin_edges = np.array([hist.GetBinLowEdge(i) for i in range(1, n_bins + 2)])
        bin_centers = np.array([hist.GetBinCenter(i) for i in range(1, n_bins + 1)])
        
        bin_contents = np.array([hist.GetBinContent(i) for i in range(1, n_bins + 1)])
        bin_errors = np.array([hist.GetBinError(i) for i in range(1, n_bins + 1)])

        bin_labels = []
        for i in range(1, n_bins + 1):
            label = hist.GetXaxis().GetBinLabel(i)
            bin_labels.append(self._clean_bin_label(label) if label else "")

        # Only use extracted labels if meaningful
        has_labels = any(label and not label.isdigit() for label in bin_labels)

        return bin_centers, bin_contents, bin_errors, bin_edges, bin_labels, has_labels

    def _clean_bin_label(self, label):
        """Clean up bin label for better readability."""
        return label.replace("TrackSelectionHighPurity", "HP")

    def _is_tex_math(self, label):
        """Whether the string is to be interpreted as a latex string."""
        if not label:
            return False
        return any(x in label for x in ('\\','~'))
        
    def _apply_axis_scaling(self, ax, xlabel, ylabel, title, logy=False, xlim=None, has_negative_values=False):
        """Apply log scaling to axes based on content and user settings."""
        if logy:
            if has_negative_values:
                output_file = getattr(self, '_current_output_path', 'unknown')
                print(f"\nWarning: Cannot set log scale for plot '{title}' (output: {output_file}) - plot contains negative y-values")
            else:
                ax.set_yscale("log")

        # Skip auto x-log-scale if the user gave an explicit xlim that includes non-positive values
        xlim_blocks_log = xlim is not None and xlim[0] <= 0

        # Automatic log scale for specific variable types
        if ( any(x in xlabel for x in (r'p_{\text{T}}', r'$p_{\text{T}}$', 'p_{T}'))
             and "pull" not in xlabel.lower() and not xlim_blocks_log ):
            ax.set_xscale("log")

        if ( (any(x in ylabel for x in (r'p_{\text{T}}', r'$p_{\text{T}}$', 'p_{T}'))
             and "pull" not in ylabel.lower()) or (title is not None and "Sigma" in title) ):
            if not has_negative_values:
                ax.set_yscale("log")

        if any(x in xlabel for x in ('vertpos', 'Track~r')) and not xlim_blocks_log:
            ax.set_xscale("log")

        return xlabel, ylabel

    def _plot_histogram_data(self, ax, bin_centers, bin_contents, bin_errors, bin_edges, 
                           label, color_idx, plot_histograms):
        """Plot histogram data either as histogram or error bars."""
        if plot_histograms:
            hep.histplot(
                bin_contents,
                bins=bin_edges,
                yerr=bin_errors,
                label=label,
                color=self.colors[color_idx % len(self.colors)],
                histtype="step",
                linewidth=2,
                ax=ax,
            )
        else:
            ax.errorbar(
                bin_centers,
                bin_contents,
                yerr=bin_errors,
                label=label,
                color=self.colors[color_idx % len(self.colors)],
                fmt=self.markers[color_idx % len(self.markers)],
                markersize=5,
                capsize=2,
                linewidth=1.5,
            )

    def _calculate_and_plot_ratio(self, ax_ratio, bin_centers, bin_contents, bin_errors, bin_edges,
                                ref_centers, ref_contents, ref_errors, color_idx, plot_histograms):
        """Calculate and plot ratio between current and reference histogram."""
        if ax_ratio is None:
            return []

        tolerance = 1e-6
        matching_indices = []

        for idx, center in enumerate(bin_centers):
            ref_idx = np.argmin(np.abs(ref_centers - center))
            if np.abs(ref_centers[ref_idx] - center) < tolerance:
                matching_indices.append((idx, ref_idx))

        if not matching_indices:
            return []

        curr_idxs, ref_idxs = zip(*matching_indices)
        matching_centers = bin_centers[list(curr_idxs)]

        matching_ref_contents = ref_contents[list(ref_idxs)]
        matching_contents = bin_contents[list(curr_idxs)]
        matching_ref_errors = ref_errors[list(ref_idxs)]
        matching_errors = bin_errors[list(curr_idxs)]

        ratio = np.divide(
            matching_contents,
            matching_ref_contents,
            out=np.zeros_like(matching_contents),
            where=matching_ref_contents != 0,
        )

        ratio_errors = np.zeros_like(ratio)
        mask = (matching_ref_contents != 0) & (matching_contents != 0)
        ratio_errors[mask] = np.abs(ratio[mask]) * np.sqrt(
            np.power(matching_errors[mask] / np.maximum(matching_contents[mask], 1e-10), 2)
            + np.power(matching_ref_errors[mask] / np.maximum(matching_ref_contents[mask], 1e-10), 2)
        )
        ratio_errors = np.nan_to_num(ratio_errors, nan=0.0, posinf=0.0, neginf=0.0)

        if plot_histograms and len(matching_centers) > 1:
            # Use the histogram's real bin edges for the matched bins instead
            # of reconstructing fake uniform-width edges from the average bin
            # spacing - the latter silently diverges from the main plot's
            # (possibly variable-width, e.g. after a variable-edge rebin)
            # bin_edges, so the two panels ended up drawing different bin
            # boundaries.
            curr_idx_arr = np.array(curr_idxs)
            is_contiguous = np.all(np.diff(curr_idx_arr) == 1)
            if is_contiguous:
                matching_edges = np.concatenate([
                    bin_edges[curr_idx_arr],
                    [bin_edges[curr_idx_arr[-1] + 1]],
                ])
            else:
                # Fallback for non-contiguous matches (gaps between matched
                # bins): fall back to the previous approximation rather than
                # constructing a misleading edges array.
                bin_widths = np.diff(np.sort(matching_centers))
                avg_width = np.mean(bin_widths)
                matching_edges = np.concatenate([
                    [matching_centers[0] - avg_width / 2],
                    matching_centers + avg_width / 2,
                ])
            hep.histplot(
                ratio, bins=matching_edges, yerr=ratio_errors,
                color=self.colors[color_idx % len(self.colors)],
                histtype="step", linewidth=2, ax=ax_ratio,
            )
        else:
            ax_ratio.errorbar(
                matching_centers, ratio, yerr=ratio_errors,
                color=self.colors[color_idx % len(self.colors)],
                fmt=self.markers[color_idx % len(self.markers)],
                markersize=5, capsize=2, linewidth=1.5,
            )

        return ratio[(ratio > 0) & np.isfinite(ratio)]

    def _wrap_legend_labels(self, labels, width=25):
        """Soft-wrap legend labels at natural break points to reduce horizontal size."""
        wrapped = []
        for lab in labels:
            # Explicit new lines
            if "\\n" in lab:
                explicit_lines = lab.split("\\n")
                wrapped.append("\n".join(explicit_lines))
                continue

            # Preserve " - " separator (for overlay labels like "File - Collection")
            if " - " in lab:
                parts = lab.split(" - ", 1)  # Split only on first occurrence
                file_part = parts[0]
                collection_part = parts[1] if len(parts) > 1 else ""

                # If the combined length is too long, put on separate lines
                if len(lab) > width:
                    wrapped.append(f"{file_part}\n- {collection_part}")
                else:
                    wrapped.append(lab)
                continue

            # Automatic split using / or _
            parts = re.split(r'(/|_)+', lab)
            tokens = []
            buffer = ""
            for p in parts:
                if not p:
                    continue
                candidate = (buffer + p) if buffer else p
                if len(candidate) > width and buffer:
                    tokens.append(buffer.rstrip("_/"))
                    buffer = p
                else:
                    buffer = candidate
            if buffer:
                tokens.append(buffer.rstrip("_/"))

            # CamelCase and digit splitting
            final_tokens = []
            for tok in tokens:
                if len(tok) > width:
                    subtoks = re.findall(r'[A-Z]?[a-z]+|[A-Z]+(?![a-z])|\d+', tok)
                    line = ""
                    for st in subtoks:
                        if len(line) + len(st) + 1 > width and line:
                            final_tokens.append(line)
                            line = st
                        else:
                            line = (line + st) if not line else (line + st)
                    if line:
                        final_tokens.append(line)
                else:
                    final_tokens.append(tok)
            wrapped_label = "\n".join(final_tokens) if final_tokens else lab
            wrapped.append(wrapped_label)
        return wrapped

    def _configure_legend(self, ax, labels, legend_title, legend_location,
                          legend_title_fontsize, legend_fontsize,
                          logy=False, force_outside=False):
        """Configure legend; wrap long entries and move outside if needed."""
        if not labels:
            return

        wrapped_labels = self._wrap_legend_labels(labels)
        max_label_length = max(len(l.replace("\n", "")) for l in wrapped_labels)

        place_outside = (
            force_outside
            or max_label_length > 30
            or (len(wrapped_labels) > 8 and max_label_length > 20)
        )

        if legend_fontsize is None:
            if len(wrapped_labels) > 10:
                legend_fontsize = 18
            elif len(wrapped_labels) > 6:
                legend_fontsize = 20
            else:
                legend_fontsize = 22

        if legend_title_fontsize is None:
            if legend_title and len(legend_title) > 30:
                legend_title_fontsize = 20
            else:
                legend_title_fontsize = 24

        if self._is_tex_math(legend_title):
            legend_title = fr'${legend_title}$'

        if place_outside:
            y_min, y_max = ax.get_ylim()
            y_range = y_max - y_min
            ax.set_ylim(y_min, y_max + y_range * 0.1)
            ax.figure.subplots_adjust(right=0.9)
            ax.legend(
                wrapped_labels,
                loc="upper left",
                bbox_to_anchor=(1.01, 1.0),
                borderaxespad=0.0,
                title=legend_title,
                fontsize=legend_fontsize,
                title_fontsize=legend_title_fontsize,
                frameon=False,
            )
            return

        legend_padding = 0.5 if len(wrapped_labels) <= 2 else np.ceil(len(wrapped_labels) / 2) * 0.25
        legend_padding = min(legend_padding, 2.0)
        if logy or ax.get_yscale() == "log":
            legend_padding *= 50
        if any("\n" in l for l in wrapped_labels):
            legend_padding *= 2
        y_min, y_max = ax.get_ylim()
        y_range = y_max - y_min
        ax.set_ylim(y_min, y_max + y_range * legend_padding)

        legend_columns = 1 if len(wrapped_labels) <= 2 else 2
        ax.legend(
            wrapped_labels,
            loc=legend_location,
            ncols=legend_columns,
            title=legend_title,
            fontsize=legend_fontsize,
            title_fontsize=legend_title_fontsize,
            columnspacing=1.0,
        )

    def _apply_custom_formatter(self, ax):
        """Apply custom scientific notation formatter to y-axis."""
        if ax.get_yscale() != "log":
            from matplotlib.ticker import ScalarFormatter

            class CustomScalarFormatter(ScalarFormatter):
                def format_data_short(self, value):
                    if self.orderOfMagnitude != 0:
                        return f"×10$^{{{self.orderOfMagnitude}}}$"
                    return ""

            formatter = CustomScalarFormatter(useOffset=True, useMathText=True)
            ax.yaxis.set_major_formatter(formatter)

            # Move y-axis scientific notation to avoid overlap with CMS label
            ax.yaxis.get_offset_text().set_position((-0.01, 1.02))
            ax.yaxis.get_offset_text().set_horizontalalignment("right")
            ax.yaxis.get_offset_text().set_verticalalignment("bottom")

    def _load_single_histogram(self, file_path, hist_path):
        """
        Load a single histogram from a ROOT file.
        
        Args:
        file_path: ROOT file path
        hist_path: Histogram path inside the ROOT file
        
        Returns:
        ROOT histogram clone, or None if not found.
        """
        try:
            root_file = _get_cached_root_file(file_path)
            if not root_file or root_file.IsZombie():
                print(f'Warning: Could not open {file_path}.')
                return None

            hist = root_file.Get(hist_path)

            if not hist:
                # A histogram can legitimately be absent from a file (e.g. a
                # collection present in only some inputs). Missing histograms are
                # summarised once per file in compare_files, so only note each one
                # individually in debug mode instead of spamming a line per plot.
                if self.debug:
                    print(f'Warning: Could not find {hist_path}.')
                return None

            hist_clone = hist.Clone()
            hist_clone.SetDirectory(0)

            return hist_clone

        except Exception as e:
            print(f"Error processing {file_path}: {e}")
            return None
        
    def _load_histograms(self, file_paths, hist_path, labels):
        """
        Load histograms from files, returning list of histograms and valid labels.
        
        Args:
        file_paths: list of ROOT files
        hist_path: histogram path (same for every file)
        labels: legend labels
        
        Returns:
        tuple: (histograms, valid_labels)
        """
        histograms = []
        valid_labels = []

        for file_path, label in zip(file_paths, labels):
            hist = self._load_single_histogram(file_path, hist_path)

            if hist is None: # a warning is already present in _load_single_histogram()
                continue

            histograms.append(hist)
            valid_labels.append(label)

        return histograms, valid_labels

    def extract_labels_from_hist(self, hist, xtitle=None, ytitle=None):
        """
        Extract title and axis labels from ROOT histogram.
        Overwrite result with user optional arguments.
        Args:
            hist: ROOT histogram

        Returns:
            tuple: (title, xlabel, ylabel)
        """
        title = hist.GetTitle()
        xlabel = title
        ylabel = "Occurrences"

        # Match "vs", "vs.", and flexible spacing/periods between v and s; allow underscores or spaces as delimiters
        vs_regex = re.compile(r'[_\s]+v\s*\.?\s*s\s*\.?[_\s]+', re.IGNORECASE)
        if vs_regex.search(title):
            # Detect explicit underscore-delimited form even with optional dots/spaces
            used_underscore_delim = bool(re.search(r'_v\s*\.?\s*s\s*\.?_', title.lower()))
            parts = vs_regex.split(title, maxsplit=1)
            left, right = parts[0].strip(), parts[1].strip()

            if used_underscore_delim:
                left = left.replace("_", " ")
                right = right.replace("_", " ")

            if "#sigma(" in title.lower():
                core = left[left.find("(") + 1 : left.rfind(")")]
                ylabel = r"$\delta$" + core
                if "p_{T}" in ylabel:
                    ylabel = ylabel + "/" + core
                right_clean = right
                if "Mean" in title:
                    ylabel = "<" + ylabel + ">"
                    right_clean = right_clean.replace("Mean", "")
                elif "Sigma" in title:
                    ylabel = r"$\sigma$(" + ylabel + ")"
                    right_clean = right_clean.replace("Sigma", "")
                xlabel = right_clean.strip()

            elif "Mean" in title:
                right_clean = right.replace("Mean", "").strip()
                ylabel = "<" + left + ">"
                xlabel = right_clean
            elif "Sigma" in title:
                right_clean = right.replace("Sigma", "").strip()
                ylabel = r"$\sigma$<" + left + ">"
                xlabel = right_clean
            else:
                # Default: "ylabel vs xlabel"
                ylabel = left
                xlabel = right
                if hist.InheritsFrom("TProfile") and "mean " in ylabel:
                    ylabel = ylabel.replace("mean ", "<") + ">"
        else:
            # Pull plots
            if "pull" not in title.lower():
                if "eta" in title.lower():
                    xlabel = r"$\eta$"
                elif "pt2" in title.lower():
                    xlabel = r"$p_{\text{T}}^2$"
                elif "pt" in title.lower():
                    xlabel = r"$p_{\text{T}}$"
                elif "phi" in title.lower():
                    xlabel = r"$\phi$"
            # Efficiency and turn-on plots
            if "eff" in title.lower():
                ylabel = "Efficiency"
            elif "turn-on" in title.lower():
                ylabel = "Turn-On"

        # Always convert ROOT (#eta, #phi, ...) notation in the automatic labels to
        # mathtext. These labels can mix hand-written mathtext ($\sigma$) with ROOT
        # tokens taken from the histogram title, so the conversion must always run;
        # skipping it (as when a label merely contains a backslash) leaves raw ROOT
        # tokens that matplotlib's mathtext parser rejects.
        title = self.convert_root_to_latex(title)
        xlabel = self.convert_root_to_latex(xlabel)
        ylabel = self.convert_root_to_latex(ylabel)

        # A user-provided --xtitle/--ytitle overrides the automatic label: wrap it as
        # math if it already looks like LaTeX, otherwise convert its ROOT notation.
        if xtitle is not None:
            xlabel = fr'${xtitle}$' if self._is_tex_math(xtitle) else self.convert_root_to_latex(xtitle)
        if ytitle is not None:
            ylabel = fr'${ytitle}$' if self._is_tex_math(ytitle) else self.convert_root_to_latex(ytitle)

        return title, xlabel, ylabel

    def convert_root_to_latex(self, text):
        """
        Convert ROOT-style notation to matplotlib mathtext format.

        Args:
            text: String with ROOT notation

        Returns:
            String with mathtext notation
        """
        conversions = [
            ("#eta", r"$\eta$"),
            ("#phi", r"$\phi$"),
            ("#theta", r"$\theta$"),
            ("#alpha", r"$\alpha$"),
            ("#beta", r"$\beta$"),
            ("#gamma", r"$\gamma$"),
            ("#delta", r"$\delta$"),
            ("#Delta", r"$\Delta$"),
            ("#sigma", r"$\sigma$"),
            ("d_{xy}", r"$d_{xy}$"),
            ("d_{z}", r"$d_{z}$"),
            ("#chi", r"$\chi$"),
            ("#nu", r"$\nu$"),
            ("#tau", r"$\tau$"),
            ("pt2", r"$p_{\text{T}}^2$"),
            ("Pt2", r"$p_{\text{T}}^2$"),
            ("pT2", r"$p_{\text{T}}^2$"),
            ("p_{T}", r"$p_{\text{T}}$"),
            ("p_{t}", r"$p_{\text{T}}$"),
            ("pT", r"$p_{\text{T}}$"),
            ("pt", r"$p_{\text{T}}$"),
            ("GeV/c^{2}", r"GeV/$c^2$"),
            ("GeV/c", r"GeV/$c$"),
            ("^{2}", r"$^2$"),
            ("number of", "#"),
            ("Number of", "#"),
        ]

        result = text
        for root_notation, mathtext_notation in conversions:
            result = result.replace(root_notation, mathtext_notation)

        # Replace standalone tokens without #
        token_map = {
            "eta": r"$\eta$",
            "phi": r"$\phi$",
            "theta": r"$\theta$",
            "alpha": r"$\alpha$",
            "beta": r"$\beta$",
            "gamma": r"$\gamma$",
            "delta": r"$\delta$",
            "sigma": r"$\sigma$",
            "chi": r"$\chi$"
        }
        # Tokens not preceded by letter/digit/_ or backslash and not followed by letter/digit/_.
        pattern = re.compile(
            r'(?<![A-Za-z0-9_\\])(' + "|".join(token_map.keys()) + r')(?![A-Za-z0-9_])',
            re.IGNORECASE,
        )

        def _token_repl(m):
            return token_map[m.group(1).lower()]

        result = pattern.sub(_token_repl, result)

        return result

    def find_histograms_in_file(self, root_file, patterns):
        """
        Find all histogram paths matching regex patterns in a ROOT file.

        Traverses the directory tree once and checks each histogram against
        all patterns, avoiding redundant I/O when multiple patterns are used.

        Args:
            root_file: Open ROOT file
            patterns: Single pattern string or list of regex pattern strings

        Returns:
            list: List of (histogram_path, matched_pattern) tuples
        """
        if isinstance(patterns, str):
            patterns = [patterns]
        compiled = [re.compile(p) for p in patterns]

        def traverse_directory(directory, current_path=""):
            """Recursively traverse ROOT file directory structure."""
            hist_paths = []
            for key in directory.GetListOfKeys():
                obj_name = key.GetName()
                obj_path = f"{current_path}/{obj_name}" if current_path else obj_name

                # Decide what to do from the class name stored in the key, without
                # deserializing the object. Only directories are read (to recurse);
                # matching histogram paths are collected without ever reading the
                # (potentially large) histograms, since only their paths are needed.
                cls = ROOT.TClass.GetClass(key.GetClassName())
                if cls is None:
                    continue
                if cls.InheritsFrom("TDirectory"):
                    hist_paths.extend(traverse_directory(key.ReadObj(), obj_path))
                # Skip 2D histograms
                elif cls.InheritsFrom("TH1") and not cls.InheritsFrom("TH2"):
                    for pat in compiled:
                        if pat.search(obj_path):
                            hist_paths.append((obj_path, pat.pattern))
                            break  # first match wins

            return hist_paths

        return traverse_directory(root_file)

    def _parse_overlay_groups(self, overlay_specs, n_files):
        """
        Parse --overlay specifications into overlay jobs.

        Each --overlay value is one job. A job is a colon-separated list of
        specs; each spec is 'collections[@file]', where collections is a
        comma-separated list of collection patterns and the optional @file
        selector is a 1-based file index or '*' (all files, the default).

        Backward compatible: '--overlay collA:collB' has no @ selectors, so
        both collections are overlaid from every file, exactly as before.
        The new form '--overlay collA,collB@1:collC@2' overlays collA and
        collB from file 1 together with collC from file 2.

        Args:
            overlay_specs: List of --overlay strings (one per --overlay flag)
            n_files: Number of input files (to resolve/validate @file indices)

        Returns:
            List of jobs. Each job is a dict with:
              'patterns':      ordered list of unique collection patterns
              'file_patterns': list (indexed by file) of the patterns to take
                               from that file, in command-line order
        """
        if not overlay_specs:
            return []

        jobs = []
        for spec in overlay_specs:
            patterns = []
            file_patterns = [[] for _ in range(n_files)]
            for chunk in spec.split(":"):
                chunk = chunk.strip()
                if not chunk:
                    continue
                collections, _, file_sel = chunk.partition("@")
                file_sel = file_sel.strip()
                collection_list = [c.strip() for c in collections.split(",") if c.strip()]
                if not collection_list:
                    continue

                if file_sel in ("", "*"):
                    targets = range(n_files)
                else:
                    try:
                        idx = int(file_sel)
                    except ValueError:
                        raise ValueError(
                            f"--overlay: invalid file selector '@{file_sel}' in '{spec}' "
                            "(use a 1-based file index or '*')."
                        )
                    if not 1 <= idx <= n_files:
                        raise ValueError(
                            f"--overlay: file index {idx} in '{spec}' is out of range "
                            f"(1..{n_files})."
                        )
                    targets = [idx - 1]

                for collection in collection_list:
                    if collection not in patterns:
                        patterns.append(collection)
                    for target in targets:
                        file_patterns[target].append(collection)

            if patterns:
                jobs.append({"patterns": patterns, "file_patterns": file_patterns})

        return jobs

    def _normalize_path_for_overlay(self, hist_path, overlay_patterns):
        """
        Normalize histogram path by removing overlay pattern components.

        Args:
            hist_path: Full histogram path
            overlay_patterns: List of patterns that should be removed for grouping

        Returns:
            Tuple (normalized_path, matched_pattern), or (None, None) if the path
            does not match any overlay pattern.
        """
        for pattern in overlay_patterns:
            if re.search(pattern, hist_path):
                # Remove the matching pattern part
                normalized = re.sub(pattern, "OVERLAY_PLACEHOLDER", hist_path)
                return normalized, pattern

        return None, None

    def _group_histograms_by_overlay(self, all_hists, job_patterns):
        """
        Group the histograms of one overlay job by their normalized path.

        Patterns are tried longest-first so that a collection name which is a
        prefix of another (e.g. 'hltFoo' vs 'hltFooBar') does not shadow the
        more specific one, letting both end up in the same group.

        Args:
            all_hists: Set of all matched histogram paths
            job_patterns: Collection patterns of a single overlay job

        Returns:
            (grouped, used) where grouped maps a normalized path to a
            {pattern: hist_path} dict, keeping only groups that span at least
            two collections; used is the set of histogram paths that ended up
            in a kept group.
        """
        ordered = sorted(job_patterns, key=len, reverse=True)
        grouped = {}

        for hist_path in all_hists:
            normalized, matched_pattern = self._normalize_path_for_overlay(
                hist_path, ordered
            )
            if normalized is None:
                continue
            grouped.setdefault(normalized, {}).setdefault(matched_pattern, hist_path)

        # Keep only groups that overlay at least two different collections.
        filtered_groups = {
            norm_path: hists
            for norm_path, hists in grouped.items()
            if len(hists) > 1
        }
        used_hists = {hp for hists in filtered_groups.values() for hp in hists.values()}

        return filtered_groups, used_hists

    def _overlay_job_folder(self, job_patterns):
        """Build a filesystem-safe folder name from an overlay job's patterns."""
        names = []
        for pat in job_patterns:
            name = re.sub(r'[^a-zA-Z0-9_]', '', pat)
            if name:
                names.append(name)
        return "+".join(sorted(set(names)))

    def _generate_overlay_output_path(self, normalized_path, base_pattern,
                                      output_dir, job_folder):
        """
        Generate the output path for an overlaid histogram.

        The overlaid collections live in a per-job folder
        (output_dir/overlay/<job_folder>/...), so the collection segment is
        dropped from the preserved directory tree.

        Args:
            normalized_path: Normalized path with the collection placeholder
            base_pattern: Source pattern used to preserve the folder structure
            output_dir: Base output directory
            job_folder: Name of this overlay job's subfolder

        Returns:
            Output file path
        """
        tree_path = normalized_path.replace("OVERLAY_PLACEHOLDER", "")
        tree_path = re.sub(r"/+", "/", tree_path).strip("/")
        overlay_dir = os.path.join(output_dir, "overlay", job_folder)
        return self._generate_output_path(tree_path, base_pattern, overlay_dir)

    def plot_comparison(
        self,
        histograms,
        labels,
        output_path,
        logy=False,
        normalize=False,
        cms_text="Preliminary",
        cms_text_fontsize=18,
        energy_text=None,
        show_energy=False,
        data=False,
        grid=False,
        plot_histograms=False,
        do_ratio=True,
        save_as_pdf=False,
        legend_outside=False,
        legend_title=None,
        legend_location=None,
        legend_title_fontsize=None,
        legend_fontsize=None,
        title=None,
        xtitle=None,
        ytitle=None,
        xtitle_fontsize=None,
        ytitle_fontsize=None,
        xlim=None,
        ylim=None,
        ylim_ratio=None,
        hline=None,
        ytitle_ratio_fontsize=None,
        ratio_label=None,
        complement=False,
        rebin=None,
        clamp_yerrors=False,
        multiply_x_values=None,
    ):
        """
        Create comparison plot with ratio panel.

        Args:
            histograms: List of ROOT histograms to compare
            labels: List of labels for each histogram
            output_path: Output file path
            logy: Use log scale for y-axis
            normalize: Normalize histograms to unit area
            cms_text: CMS label text
            energy_text: Custom energy text (if None, uses default)
            show_energy: Whether to show energy text
            grid: Whether to show grid on both main and ratio plots
            plot_histograms: Whether to plot histograms instead of individual bin points with errors
            pdf: Whether to save as PDF in addition to PNG
            legend_outside: Force legend to be placed outside the plot area
            hline: Add a horizontal dashed line at the specified y axis value
            rebin: Optional rebinning specification, either a single int or an array-like of new bin edges
        """
        if len(histograms) != len(labels):
            raise ValueError("Number of histograms must match number of labels")

        if title is None:
            title, xlabel, ylabel = self.extract_labels_from_hist(histograms[0], xtitle, ytitle)
        else:
            _, xlabel, ylabel = self.extract_labels_from_hist(histograms[0], xtitle, ytitle)

        if legend_title is None: # otherwise use the user's choice
            if hasattr(self, "_current_output_path"):
                output_parts = self._current_output_path.replace("\\", "/").split("/")
                if len(output_parts) > 1:
                    # Get the parent directory name (second to last part)
                    legend_title = output_parts[-2] if len(output_parts) >= 2 else None
                    # Ignore legend name for summary plots
                    if "global" in output_parts[-1].lower() or "coll" in output_parts[-1].lower():
                        legend_title = None
                    # Split long overlay titles at + marker
                    elif legend_title and "+" in legend_title and len(legend_title) > 30:
                        legend_title = legend_title.replace("+", "\n+")

        fig = plt.figure(figsize=self.figsize)
        gs = gridspec.GridSpec(
            2, 1, height_ratios=[1 - self.ratio_height, self.ratio_height], hspace=0.10
        )

        ax_main = fig.add_subplot(gs[0])

        # CMS styling
        if show_energy:
            if energy_text is None:
                # Default energy text (13 TeV)
                hep.cms.label(cms_text, data=data, ax=ax_main, fontsize=cms_text_fontsize)
            else:
                hep.cms.label(cms_text, data=data, ax=ax_main, fontsize=cms_text_fontsize,
                              rlabel=fr"${energy_text}$" if self._is_tex_math(energy_text) else energy_text)
        else:
            hep.cms.label(cms_text, data=data, ax=ax_main, rlabel="", fontsize=cms_text_fontsize)

        ax_ratio = None
        if do_ratio and len(histograms) > 1:
            ax_ratio = fig.add_subplot(gs[1], sharex=ax_main)

        ref_centers = None
        ref_contents = None
        ref_errors = None
        has_labels = False
        bin_labels = None
        has_negative_values = False

        all_ratios = []
        ref_centers = ref_contents = ref_errors = None

        for i, (hist, label) in enumerate(zip(histograms, labels)):
            bin_centers, bin_contents, bin_errors, bin_edges, bin_labels, has_labels = (
                self.root_to_numpy(hist, rebin=rebin)
            )

            if multiply_x_values is not None:
                bin_centers = bin_centers * multiply_x_values
                bin_edges = bin_edges * multiply_x_values

            if normalize and np.sum(bin_contents) > 0:
                norm_factor = 1.0 / np.sum(bin_contents)
                bin_contents *= norm_factor
                bin_errors *= norm_factor

            if complement:
                bin_contents = np.array([1-i for i in bin_contents])

            # clamp_yerrors only affects how the error bars are drawn on the main plot,
            # so its (2, N) asymmetric output is kept in a separate variable used just for display.
            plot_errors = bin_errors
            if clamp_yerrors:
                plot_errors = self._clamp_yerrors(bin_contents, bin_errors)

            if np.any(bin_contents < 0):
                has_negative_values = True

            # Use first histogram as reference
            if i == 0:
                ref_centers, ref_contents, ref_errors = bin_centers, bin_contents, bin_errors

            self._plot_histogram_data(ax_main, bin_centers, bin_contents, plot_errors, 
                                    bin_edges, label, i, plot_histograms)

            if i > 0:
                valid_ratios = self._calculate_and_plot_ratio(
                    ax_ratio, bin_centers, bin_contents, bin_errors, bin_edges,
                    ref_centers, ref_contents, ref_errors, i, plot_histograms
                )
                all_ratios.extend(valid_ratios)

        # Styling for main plot
        # Calculate visible label length by removing LaTeX notation for sizing
        visible_ylabel = re.sub(r"\$[^$]*\$", "X", ylabel)
        label_size = "medium" if len(visible_ylabel) < 20 else "small"
        if title:
            ax_main.set_title(title, pad=50)

        xlabel, ylabel = self._apply_axis_scaling(ax_main, xlabel, ylabel, title, logy, xlim, has_negative_values)

        ax_main.set_xlabel(xlabel, fontsize=label_size if xtitle_fontsize is None else xtitle_fontsize)
        ax_main.set_ylabel(ylabel, fontsize=label_size if ytitle_fontsize is None else ytitle_fontsize)

        ha = "right"
        # args.source is None in --operation mode; guard before indexing it.
        if args.source and "KindOfSignalPV" in args.source[0][0]:
            if rebin != None:
                bin_labels = ['PV not Reco', 'PV Reco\nas Leading', 'PV Reco\nnot as Leading']
                ha = "center"
                
        if has_labels: # Set custom labels if available
            ax_main.set_xticks(
                ref_centers,
                bin_labels,
                size="small" if len(bin_labels) < 10 else "xx-small",
                rotation=45,
                ha=ha,
                va="top",
            )
            ax_main.tick_params(axis="x", which="minor", bottom=False)

        self._configure_legend(ax_main, labels, legend_title, legend_location,
                               legend_title_fontsize, legend_fontsize,
                               logy, force_outside=legend_outside)

        if grid:
            ax_main.grid(True, alpha=0.75, linestyle="dashdot", linewidth=0.75)           
        if xlim:
            ax_main.set_xlim(*xlim)
        if ylim:
            ax_main.set_ylim(*ylim)

        if hline:
            ax_main.axhline(y=1, color="black", linestyle="--", alpha=0.7)
        
        self._apply_custom_formatter(ax_main)

        # Ratio plot styling
        if ax_ratio is not None:

            # Only show label in ratio and keep the same ticks as main plot
            ax_main.set_xlabel("")
            ax_main.tick_params(axis="x", labelbottom=False)

            # Set ratio plot limits
            if ylim_ratio:
                ax_ratio.set_ylim(*ylim_ratio)
            elif len(all_ratios) > 0:
                # Use percentiles to avoid extreme outliers
                ratio_min = np.percentile(all_ratios, 5)
                ratio_max = np.percentile(all_ratios, 95)

                # Set range around 1 in case of small fluctuations
                ratio_min = min(ratio_min, 0.95)
                ratio_max = max(ratio_max, 1.05)

                # Add at least 15% padding
                ratio_range = ratio_max - ratio_min
                padding = max(0.15, ratio_range * 0.1)

                # Clamp limits between 0.1 and 5.0
                final_min = max(0.1, ratio_min - padding)
                final_max = min(5.0, ratio_max + padding)

                ax_ratio.set_ylim(final_min, final_max)
            else:
                ax_ratio.set_ylim(0.5, 1.5)

            if has_labels:
                ax_ratio.set_xticks(
                    ref_centers,
                    bin_labels,
                    size="small" if len(bin_labels) < 10 else "xx-small",
                    rotation=0,
                    ha="center",
                    va="top",
                )
                ax_ratio.tick_params(axis="x", which="minor", bottom=False)
            else:
                ax_ratio.set_xlabel(xlabel, fontsize=label_size if xtitle_fontsize is None else xtitle_fontsize)

            if ratio_label is None:
                ratio_ylabel = "Ratio"
            elif ratio_label == "auto":
                ratio_ylabel = f"Ratio wrt. {labels[0]}"
            else:
                ratio_ylabel = ratio_label
            ax_ratio.set_ylabel(ratio_ylabel, fontsize=label_size if ytitle_ratio_fontsize is None else ytitle_ratio_fontsize)
            ax_ratio.axhline(y=1, color="black", linestyle="--", alpha=0.7)

            if grid:
                ax_ratio.grid(True, alpha=0.75, linestyle="dashdot", linewidth=0.75)

        os.makedirs(os.path.dirname(output_path), exist_ok=True)
        plt.savefig(output_path, dpi=300, bbox_inches="tight")

        if save_as_pdf:
            pdf_path = output_path.rsplit(".", 1)[0] + ".pdf"
            plt.savefig(pdf_path, dpi=300, bbox_inches="tight")
        plt.close()

    def compare_files(self, file_paths, hist_pattern=None, operation_specs=None, operation_type=None,
                      labels=None, output_dir="plots", more_files=False, create_web_index=False,
                      overlay_groups=None, overlay_individual=True, debug=False, **plot_kwargs):
        """
        Compare histograms matching a pattern from multiple files.

        Args:
            file_paths: List of ROOT file paths
            hist_pattern: Regex pattern to match histogram paths
            labels: Labels for each file (uses filename if None)
            output_dir: Output directory for plots
            more_files: Color palette will be extended to support more than 10 files if True
            create_web_index: Whether to create index.php files for web viewing
            overlay_groups: List of pattern groups to overlay
            overlay_individual: Whether to also create individual plots for overlaid collections
            **plot_kwargs: Additional arguments for plot_comparison
        """
        if more_files:
            self._extend_color_palette(len(file_paths))
        else:
            if len(file_paths) > 10:
                raise ValueError(
                    f"Default color palette supports a maximum of 10 files, got {len(file_paths)}. "
                    "Please reduce the number of input files or use the option --more-files if you are really sure of what you are doing."
                )

        # Normalize patterns to a list
        patterns = (
            list(hist_pattern)
            if isinstance(hist_pattern, (list, tuple, set))
            else [hist_pattern]
        )

        if labels is None:
            labels = []
            for f in file_paths:
                legend_label = Path(f).stem
                if legend_label.startswith("DQM_"):
                    legend_label = legend_label.replace("DQM_", "")
                labels.append(legend_label)

        if operation_specs is not None:
         
            plot_tasks = []
            for numerator, denominator in operation_specs:
                output_path = os.path.join(output_dir, os.path.basename(numerator) + '_' + operation_type.upper() + '_' + os.path.basename(denominator))
                plot_tasks.append({
                    "is_operation": True,
                    "file_paths": file_paths,
                    "numerator_path": numerator,
                    "denominator_path": denominator,
                    "operation_type": operation_type,
                    "labels": labels,
                    "colors": self.colors,
                    "output_path": output_path,
                    "plot_kwargs": plot_kwargs,
                })

            self._process_files(plot_tasks)

            if create_web_index:
                print("Adding php index files for web viewing...")
                self._create_web_index_files(output_dir)
         
            return
        
        all_matching_hists = set()
        hist_to_pattern = {}

        print(f"Scanning all files for histograms matching the pattern(s) using {self.processors} parallel processes...")
        scan_tasks = [{"file_path": fp, "patterns": patterns} for fp in file_paths]
        with multiprocessing.Pool(self.processors) as pool:
            scan_file_task_extra_args = partial(scan_file_task, debug=self.debug)
            scan_results = pool.map(scan_file_task_extra_args, scan_tasks)

        for file_hists, file_hist_to_pattern in scan_results:
            all_matching_hists.update(file_hists)
            for hp, pat in file_hist_to_pattern.items():
                hist_to_pattern.setdefault(hp, pat)

        if not all_matching_hists:
            joined = ", ".join(map(str, patterns))
            print(f"No histograms matching pattern(s) '{joined}' found in any file")
            return

        print(
            f"Total unique histograms found across all files: {len(all_matching_hists)}"
        )

        # Warn once (not once per plot) about matched histograms missing from some
        # input files, e.g. a collection present in only a subset of the files. Those
        # histograms are simply skipped in the affected plots.
        for fp, (file_hists, _) in zip(file_paths, scan_results):
            missing = all_matching_hists - set(file_hists)
            if not missing:
                continue
            collections = sorted({hp.rsplit("/", 1)[0] for hp in missing})
            print(
                f"Warning: {os.path.basename(fp)} is missing {len(missing)} of the matched "
                f"histogram(s) in {len(collections)} collection(s); they are skipped in the "
                f"affected plots:"
            )
            for collection in collections:
                print(f"    {collection}")

        plot_tasks = []

        # Handle overlay groups if specified
        if overlay_groups:
            used_folders = set()
            all_used = set()

            for job_index, job in enumerate(overlay_groups):
                job_patterns = job["patterns"]
                file_patterns = job["file_patterns"]
                job_legend = job.get("legend")
                job_legend_title = job.get("legend_title")

                # Per-job output folder, kept unique across jobs.
                job_folder = self._overlay_job_folder(job_patterns) or f"overlay{job_index + 1}"
                base_folder = job_folder
                dup = 2
                while job_folder in used_folders:
                    job_folder = f"{base_folder}_{dup}"
                    dup += 1
                used_folders.add(job_folder)

                grouped_hists, used_hists = self._group_histograms_by_overlay(
                    all_matching_hists, job_patterns
                )
                all_used |= used_hists

                if not grouped_hists:
                    print(f"No overlappable histogram groups found for overlay 'overlay/{job_folder}'")
                    continue
                print(f"Found {len(grouped_hists)} histogram groups to overlay in 'overlay/{job_folder}'")

                for normalized_path, pattern_to_path in grouped_hists.items():
                    base_pattern = hist_to_pattern.get(
                        next(iter(pattern_to_path.values())), patterns[0]
                    )
                    output_path = self._generate_overlay_output_path(
                        normalized_path, base_pattern, output_dir, job_folder
                    )

                    # Build the overlaid series in file-major order: for each file,
                    # take the collections assigned to it in command-line order. A
                    # custom --overlay-legend maps to these (file, collection) slots
                    # in the same order.
                    series_files, series_paths, series_labels = [], [], []
                    slot = 0
                    for file_idx in range(len(file_paths)):
                        for pattern in file_patterns[file_idx]:
                            hist_path = pattern_to_path.get(pattern)
                            if hist_path is not None:
                                if job_legend is not None and slot < len(job_legend):
                                    label = job_legend[slot]
                                else:
                                    collection = re.search(pattern, hist_path)
                                    collection_name = collection.group(0) if collection else pattern
                                    label = f"{labels[file_idx]} - {collection_name}"
                                series_files.append(file_paths[file_idx])
                                series_paths.append(hist_path)
                                series_labels.append(label)
                            slot += 1

                    if len(series_paths) < 2:
                        continue

                    task_plot_kwargs = plot_kwargs
                    if job_legend_title is not None:
                        task_plot_kwargs = dict(plot_kwargs)
                        task_plot_kwargs["legend_title"] = job_legend_title

                    plot_tasks.append({
                        "file_paths": series_files,
                        "hist_path": series_paths,
                        "labels": series_labels,
                        "colors": self.colors,
                        "output_path": output_path,
                        "plot_kwargs": task_plot_kwargs,
                        "is_overlay": True,
                    })

            # Also make individual plots unless disabled.
            if overlay_individual:
                remaining_hists = all_matching_hists
                if remaining_hists:
                    print(f"Processing {len(remaining_hists)} histograms (including individual overlay collections)")
            else:
                remaining_hists = all_matching_hists - all_used
                if remaining_hists:
                    print(f"Processing {len(remaining_hists)} non-overlaid histograms separately")
        else:
            remaining_hists = all_matching_hists

        # Create tasks for non-overlaid histograms (standard behavior)
        for hist_path in sorted(remaining_hists):
            matched_pat = hist_to_pattern.get(hist_path, patterns[0])
            output_path = self._generate_output_path(hist_path, matched_pat, output_dir)

            plot_tasks.append({
                "file_paths": file_paths,
                "hist_path": hist_path,
                "labels": labels,
                "colors": self.colors,
                "output_path": output_path,
                "plot_kwargs": plot_kwargs,
                "is_overlay": False,
            })

        self._process_files(plot_tasks)

        if create_web_index:
            print("Adding php index files for web viewing...")
            self._create_web_index_files(output_dir)

    def _generate_output_path(self, hist_path, pattern, output_dir):
        """
        Generate output path preserving folder structure from regex match.

        Args:
            hist_path: Full histogram path from ROOT file
            pattern: Regex pattern used to find the histogram
            output_dir: Base output directory

        Returns:
            Output file path with preserved folder structure
        """
        match = re.search(pattern, hist_path)
        if match:
            match_end = match.end()
            remaining_path = hist_path[match_end:].lstrip("/")

            if remaining_path:
                path_parts = remaining_path.split("/")

                if len(path_parts) > 1:
                    # Create subdirectories and filename
                    subdirs = "/".join(path_parts[:-1])
                    filename = path_parts[-1]
                    output_path = os.path.join(output_dir, subdirs, f"{filename}")
                else:
                    # Single file, no subdirectories
                    filename = os.path.basename(hist_path)
                    output_path = os.path.join(output_dir, f"{filename}")
            else:
                # If no remaining path, use the last part of the matched path
                matched_part = match.group(0)
                filename = matched_part.split("/")[-1]
                output_path = os.path.join(output_dir, f"{filename}")
        else:
            hist_name = hist_path.replace("/", "_").replace("\\", "_")
            output_path = os.path.join(output_dir, f"{hist_name}")

        return output_path

    def _create_web_index_files(self, output_dir):
        """Create index.php files in each subdirectory for web viewing."""
        index_php_url = "https://cernbox.cern.ch/remote.php/dav/public-files/XCmC5GCFnF7Aqfd/index.php"

        try:
            with urllib.request.urlopen(index_php_url) as response:
                index_php_content = response.read().decode("utf-8")
        except urllib.error.URLError as e:
            print(f"Warning: Could not download index.php from CERNBox: {e}")
            print("Skipping web index creation.")
            return

        n_index = 0
        for root, dirs, files in os.walk(output_dir):
            index_php_path = os.path.join(root, "index.php")
            with open(index_php_path, "w", encoding="utf-8") as f:
                f.write(index_php_content)
            n_index += 1
        print(f"Created index.php files in {n_index} directories for web viewing")

    def _process_files(self, plot_tasks):
        """Process plotting tasks using multiprocessing or sequentially in debug mode."""
        if self.debug:
            print("Debug mode: processing sequentially")
            completed = 0
            failed = 0
            for i, task in enumerate(plot_tasks):
                task_copy = task.copy()
                task_copy.pop("progress_counter", None)
                task_copy.pop("total_tasks", None)

                try:
                    process_histogram_task(task_copy, self.debug)
                    completed += 1
                    if completed % 10 == 0 or completed == len(plot_tasks):
                        print(
                            f"\rProgress: {completed}/{len(plot_tasks)} plots completed",
                            end="",
                            flush=True,
                        )
                except Exception as e:
                    print(f"\nError in task {i}: {e}")
                    failed += 1
        else:
            print(f"Using {self.processors} parallel processes")
            with multiprocessing.Pool(self.processors) as pool:
                results = []
                for i, task in enumerate(plot_tasks):
                    task_copy = task.copy()
                    task_copy.pop("progress_counter", None)
                    task_copy.pop("total_tasks", None)
                    result = pool.apply_async(process_histogram_task, (task_copy, self.debug))
                    results.append((i, result))
                completed = 0
                failed = 0
                for i, result in results:
                    try:
                        result.get(timeout=300)
                        completed += 1
                        if completed % 10 == 0 or completed == len(plot_tasks):
                            print(
                                f"\rProgress: {completed}/{len(plot_tasks)} plots completed",
                                end="",
                                flush=True,
                            )
                    except multiprocessing.TimeoutError:
                        print(f"\nTask {i} timed out")
                        failed += 1
                    except Exception as e:
                        print(f"\nError in task {i}: {e}")
                        failed += 1
                pool.close()
                pool.join()
     
        if failed > 0:
            print(f"\n{failed} tasks failed out of {len(plot_tasks)} total")
        print()

def scan_file_task(task, debug):
    """Scan a single ROOT file for histograms matching patterns, in a separate process."""
    try:
        file_path = task["file_path"]
        patterns = task["patterns"]

        finder = DQMPlotter(debug=debug)
        root_file = ROOT.TFile.Open(file_path, "READ")
        if not root_file or root_file.IsZombie():
            print(f"Error: Could not open {file_path}")
            return [], {}

        results = finder.find_histograms_in_file(root_file, patterns)
        root_file.Close()

        file_hists = [hp for hp, _ in results]
        file_hist_to_pattern = {}
        for hp, pat in results:
            file_hist_to_pattern.setdefault(hp, pat)

        print(f"Found {len(file_hists)} matching histograms in {file_path}")
        return file_hists, file_hist_to_pattern

    except Exception as e:
        print(f"Error scanning {task.get('file_path', '?')}: {e}")
        return [], {}


def process_histogram_task(task, debug):
    """Process a single histogram task in a separate process."""
    try:
        file_paths = task["file_paths"]
        labels = task["labels"]
        output_path = task["output_path"]
        plot_kwargs = task["plot_kwargs"]
        is_overlay = task.get("is_overlay", False)
        is_operation = task.get("is_operation", False)

        plotter = DQMPlotter(debug=debug)
        plotter._current_output_path = output_path
        plotter.colors = task["colors"]
        
        if is_operation:
            numhistos, denhistos = [], []
            valid_labels = []

            for fp, lab in zip(file_paths, labels):
                num, lab = plotter._load_histograms([fp], task["numerator_path"], [lab])
                den, _ = plotter._load_histograms([fp], task["denominator_path"], [lab])

                # _load_histograms returns ([], []) when a histogram is missing; skip
                # this file so numerator/denominator/labels stay aligned.
                if not num or not den:
                    continue

                numhistos.extend(num)
                denhistos.extend(den)
                valid_labels.extend(lab)

            if len(numhistos) > 0:
                results = plotter.operate_histograms(numhistos, denhistos, task["operation_type"])
                plotter.plot_comparison(results, valid_labels, output_path, **plot_kwargs)
            return

        hist_path = task["hist_path"]
            
        if is_overlay:
            # hist_path is a list of paths for overlay mode
            histograms = []
            valid_labels = []

            if isinstance(hist_path, list):
                # Process in order to maintain consistent color/marker assignment
                for fp, hp, lab in zip(file_paths, hist_path, labels):
                    hists, vlabs = plotter._load_histograms([fp], hp, [lab])
                    histograms.extend(hists)
                    valid_labels.extend(vlabs)
            else:
                histograms, valid_labels = plotter._load_histograms(
                    file_paths, hist_path, labels
                )
        else:
            # Keep user-specified order for files
            histograms, valid_labels = plotter._load_histograms(
                file_paths, hist_path, labels
            )

        if len(histograms) > 0:
            plotter.plot_comparison(histograms, valid_labels, output_path, **plot_kwargs)

    except Exception as e:
        print(f"Error in process_histogram_task: {e}")
        raise


if __name__ == "__main__":
    """Command line interface for DQM plotter."""
    parser = argparse.ArgumentParser(
        description="Compare DQM histograms with CMS styling",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog=(
            "argument order:\n"
            "  Pass the positional ROOT files first, then the options, e.g.\n"
            "    dqm-plot file1.root file2.root -s \"path/*\" --ratio-label\n"
            "  Several options (-s/--source, --operation, --overlay, --overlay-legend,\n"
            "  --legend, --ratio-label) take a variable number of values and would\n"
            "  otherwise swallow trailing positional files. To list the files last,\n"
            "  separate them with '--':\n"
            "    dqm-plot -s \"path/*\" --ratio-label -- file1.root file2.root\n"
        ),
    )
    parser.add_argument("files", nargs="+", help="ROOT files to compare")
    parser.add_argument(
        "-s",
        "--source",
        nargs="+",
        action="append",
        help="Histogram path(s) or regex pattern(s).",
    )

    parser.add_argument(
        "--operation",
        type=str,
        nargs=2,
        help="Operate on two histogram specifications. Each argument uses the same syntax accepted by --source.",
    )
    parser.add_argument(
        "--operation-type",
        choices=["division"],
        default="division",
        help="Operation applied between the two histograms.",
    )

    parser.add_argument(
        '-l', '--legend',
        help="Comma-separated list of legend entries for each file "
        "(e.g. -l \"File 1,File 2\"). Use \\n for explicit line breaks.",
    )
    parser.add_argument("-o", "--output-dir", default="plots", help="Output directory.")
    parser.add_argument("--logy", action="store_true", help="Use log scale for y-axis.")
    parser.add_argument("--normalize", action="store_true", help="Normalize histograms.")
    parser.add_argument(
        "--cms-text",
        default="Preliminary",
        help="CMS label text (e.g. Preliminary, Work in progress).",
    )
    parser.add_argument(
        "--cms-text-fontsize",
        default=20,
        type=int,
        help="CMS label text fontsize.",
    )
    parser.add_argument(
        "--energy-text",
        default=None,
        help="Custom energy text in the top right corner (e.g., '14 TeV', 'Run 3').",
    )
    parser.add_argument(
        "--data", action="store_true", help="Use CMS datastyle for plots."
    )
    parser.add_argument(
        "--no-energy", action="store_true", help="Don't show energy text."
    )
    parser.add_argument("--no-ratio", action="store_true", help="Disable ratio plots.")
    parser.add_argument(
        "--no-grid", action="store_true", help="Remove grid from all plots."
    )
    parser.add_argument(
        '--multiply-x-values', type=float, required=False, default=None, dest='multiply_x_values',
        help="Multiply all x-axis values (bin centers and edges) by this fixed factor."
    )
    parser.add_argument(
        '--hline', type=float, required=False, default=None, dest='hline',
        help="Add a horizontal dashed line at the specified Y axis location."
    )    
    parser.add_argument(
        '--title', required=False, default=None,
        help="Plot title."
    )
    parser.add_argument(
        '--xtitle', required=False, default=None,
        help="Label of the x axis."
    )
    parser.add_argument(
        '--ytitle', required=False, default=None,
        help="Label of the y axis."
    )
    parser.add_argument(
        '--xtitle-fontsize', type=int, required=False, default=None,
        help="X axis label font size."
    )
    parser.add_argument(
        '--xlim', nargs=2, type=float, required=False,
        help="Numerical limits of the x axis of the main and ratio histograms."
    )
    parser.add_argument(
        '--ytitle-fontsize', type=int, required=False, default=None,
        help="Y axis label font size."
    )
    parser.add_argument(
        '--ylim', nargs=2, type=float, required=False,
        help="Numerical limits of the y axis of the main histogram."
    )
    parser.add_argument(
        '--ylim-ratio', nargs=2, type=float, required=False,
        help="Numerical limits of the y axis of the ratio histogram."
    )
    parser.add_argument(
        '--ytitle-ratio-fontsize', type=int, required=False, default=None,
        help="Bottom ratio plot Y axis label font size."
    )
    parser.add_argument(
        '--ratio-label', type=str, nargs='?', const='auto', default=None,
        help="Y axis label of the ratio panel. Omit the flag for 'Ratio'; pass "
        "'--ratio-label' (or '--ratio-label auto') for 'Ratio wrt <first legend "
        "label>'; or pass '--ratio-label VALUE' for a custom label. Since the bare "
        "flag takes an optional value, keep the input ROOT files away from it "
        "(place them first, or after a '--' separator)."
    )
    parser.add_argument(
        "--histograms",
        action="store_true",
        help="Plot histograms instead of points (default: points with error bars).",
    )
    parser.add_argument(
        "--pdf", action="store_true", help="Also save plots in PDF format."
    )
    parser.add_argument(
        "--more-files",
        action="store_true",
        help="Allow more than 10 input files (legend might overlap with the plots, color palette only supports 10 files).",
    )
    parser.add_argument(
        "-n",
        "--nProcessors",
        type=int,
        default=4,
        help="Number of parallel processes to use (default: 4). Set to 0 or negative to use all available threads.",
    )
    parser.add_argument(
        "--web",
        action="store_true",
        help="Create index.php files for web viewing (downloads from CernBox).",
    )
    parser.add_argument(
        "--legend-outside",
        action="store_true",
        help="Force legend to be placed outside the plot area",
    )
    parser.add_argument(
        "--legend-title",
        help="Legend title.",
    )
    parser.add_argument(
        "--legend-location",
        default="upper right",
        choices=list(mlegend.Legend.codes.keys()),
        help="Legend location.",
    )
    parser.add_argument(
        "--legend-title-fontsize",
        required=False,
        default=None,
        type=int,
        help="Legend title fontsize.",
    )
    parser.add_argument(
        "--legend-fontsize",
        required=False,
        default=None,
        type=int,
        help="Legend label fontsize.",
    )
    parser.add_argument(
        "--overlay",
        action="append",
        help="Overlay histograms from different collections into one plot. Each "
        "--overlay is one overlay job written to its own plots/overlay/<name>/ "
        "folder, and is a colon-separated list of 'collections[@file]' specs: "
        "comma-separated collection patterns, each optionally bound to a 1-based "
        "file index (or '*'/omitted for all files). "
        "'--overlay collA:collB' overlays both collections from every file; "
        "'--overlay collA@1:collB@2' overlays collA from file 1 with collB from "
        "file 2. Can be used multiple times for several overlay jobs. "
        "By default also produces the plots that would have been produced "
        "without setting this option. "
        "Use --no-overlay-individual if you want to plot only the combinations "
        "specified in this parameter."
    )
    parser.add_argument(
        "--overlay-legend",
        action="append",
        help="Comma-separated legend labels for an --overlay job, in file-major "
        "order (file 1 collections first, then file 2, ...). Pair one "
        "--overlay-legend with each --overlay, in the same order.",
    )
    parser.add_argument(
        "--overlay-legend-title",
        action="append",
        help="Legend title for an --overlay job's plots. Pair one "
        "--overlay-legend-title with each --overlay, in the same order.",
    )
    parser.add_argument(
        "--no-overlay-individual",
        action="store_true",
        help="Disable creation of individual plots for collections that are overlaid (default: create both overlay and individual plots).",
    )
    parser.add_argument(
        "--complement",
        action="store_true",
        help="Apply the 1-j operation to every bin j of the histogram before plotting it. Useful for converting purity plots into fake rates.",
    )
    parser.add_argument(
        "--rebin",
        type=str,
        default=None,
        help="Rebin histograms before plotting. Pass a single integer to merge "
        "that many existing bins together (e.g. '--rebin 2' merges every 2 "
        "bins into 1), or a comma-separated list of numbers to define new "
        "variable-width bin edges (e.g. '--rebin 0,10,20,50,100'). Must be "
        "a single shell argument (quote it if it contains spaces) since "
        "this flag takes exactly one value - this keeps it from swallowing "
        "the positional ROOT file arguments that follow it.",
    )
    parser.add_argument(
        "--clamp_yerrors",
        action="store_true",
        help="Avoids y errors to lie above 1 or below 0. "
        "To be used when plotting efficiencies. "
        "This is NOT the correct approach, but can be sufficient in many cases. "
        "The DQMGenericClient CMSSW module does not support assymetric errors, "
        "so most DQM files will be incorrect. "
        "If you can, consider instead Clopper-Pearson errors by manually computing "
        "the errors using the numerator and denominator histograms, if available."
    )
    parser.add_argument(
        "--debug",
        action="store_true",
        help="Run in debug mode: converts CPU-parallel calls to sequential steps.",
    )

    args = parser.parse_args()

    if args.operation and args.source:
        parser.error("Cannot use --operation and --source together")

    if not args.operation and not args.source:
        parser.error("Either --source or --operation must be specified")

    # '--ratio-label' takes an optional value (nargs='?'); if it accidentally consumed
    # an input ROOT file, flag it clearly instead of silently dropping that file.
    if args.ratio_label not in (None, "auto") and (
        args.ratio_label.endswith(".root") or os.path.exists(args.ratio_label)
    ):
        parser.error(
            f"--ratio-label consumed '{args.ratio_label}', which looks like an input "
            "file. Use '--ratio-label auto', quote a custom label, or separate the "
            "ROOT files with '--' (e.g. '... --ratio-label -- file1.root file2.root')."
        )

    if args.legend is not None:
        # --legend is a single comma-separated string, e.g. -l "File 1,File 2".
        args.legend = [entry.strip() for entry in args.legend.split(",")]
        if args.source is not None and len(args.legend) != len(args.files):
            print(args.legend)
            print(args.files)
            raise ValueError(f"Number of legend entries ({len(args.legend)}) must match number of files ({len(args.files)})")

    if args.nProcessors < 1:
        args.nProcessors = len(os.sched_getaffinity(0))

    plotter = DQMPlotter(processors=args.nProcessors, debug=args.debug)

    # Ensure the first axis limit lies below the second
    if args.xlim is not None and args.xlim[0] >= args.xlim[1]:
        parser.error(f"--xlim: lower bound ({args.xlim[0]}) must be <= upper bound ({args.xlim[1]})")
    if args.ylim is not None and args.ylim[0] >= args.ylim[1]:
        parser.error(f"--ylim: lower bound ({args.ylim[0]}) must be <= upper bound ({args.ylim[1]})")
    if args.ylim_ratio is not None and args.ylim_ratio[0] >= args.ylim_ratio[1]:
        parser.error(f"--ylim_ratio: lower bound ({args.ylim_ratio[0]}) must be <= upper bound ({args.ylim_ratio[1]})")

    rebin_arg = None
    if args.rebin is not None:
        rebin_tokens = [tok for tok in re.split(r"[,\s]+", args.rebin.strip()) if tok]
        if not rebin_tokens:
            parser.error("--rebin: no values found.")
        try:
            rebin_values = [float(tok) for tok in rebin_tokens]
        except ValueError:
            parser.error(f"--rebin: could not parse '{args.rebin}' as number(s).")

        if len(rebin_values) == 1:
            rebin_value = rebin_values[0]
            if not rebin_value.is_integer():
                parser.error("--rebin: a single value must be an integer number of bins to merge.")
            rebin_arg = int(rebin_value)
        else:
            rebin_arg = rebin_values

    # Flatten list of patterns
    source_patterns = None
    if args.source:
        source_patterns = [p for group in args.source for p in group]

    operation_specs = None
    if args.operation:
        operation = [args.operation]
        operation_specs = [tuple(op) for op in operation]

    # Parse overlay groups
    overlay_groups = (
        plotter._parse_overlay_groups(args.overlay, len(args.files)) if args.overlay else None
    )

    if (args.overlay_legend or args.overlay_legend_title) and not overlay_groups:
        parser.error("--overlay-legend/--overlay-legend-title require --overlay")

    # Attach per-job legend labels and titles, paired with --overlay by order.
    overlay_legends = args.overlay_legend or []
    overlay_titles = args.overlay_legend_title or []
    for job_index, job in enumerate(overlay_groups or []):
        if job_index < len(overlay_legends) and overlay_legends[job_index] is not None:
            labels = [entry.strip() for entry in overlay_legends[job_index].split(",")]
            expected = sum(len(patterns) for patterns in job["file_patterns"])
            if len(labels) != expected:
                parser.error(
                    f"--overlay-legend #{job_index + 1} has {len(labels)} labels but its "
                    f"overlay job produces {expected} series"
                )
            job["legend"] = labels
        if job_index < len(overlay_titles) and overlay_titles[job_index] is not None:
            job["legend_title"] = overlay_titles[job_index]

    plotter.compare_files(
        file_paths=args.files,
        hist_pattern=source_patterns,
        labels=args.legend,
        operation_specs=operation_specs,
        operation_type=args.operation_type,
        output_dir=args.output_dir,
        logy=args.logy,
        normalize=args.normalize,
        cms_text=args.cms_text,
        cms_text_fontsize=args.cms_text_fontsize,
        energy_text=args.energy_text,
        show_energy=not args.no_energy or args.energy_text is not None,
        data=args.data,
        do_ratio=not args.no_ratio,
        grid=not args.no_grid,
        title=args.title,
        xtitle=args.xtitle,
        xtitle_fontsize=args.xtitle_fontsize,
        ytitle=args.ytitle,
        ytitle_fontsize=args.ytitle_fontsize,
        xlim=args.xlim,
        ylim=args.ylim,
        ylim_ratio=args.ylim_ratio,
        ytitle_ratio_fontsize=args.ytitle_ratio_fontsize,
        ratio_label=args.ratio_label,
        plot_histograms=args.histograms,
        save_as_pdf=args.pdf,
        more_files=args.more_files,
        create_web_index=args.web,
        legend_outside=args.legend_outside,
        legend_title=args.legend_title,
        legend_location=args.legend_location,
        legend_title_fontsize=args.legend_title_fontsize,
        legend_fontsize=args.legend_fontsize,
        overlay_groups=overlay_groups,
        overlay_individual=not args.no_overlay_individual,
        complement=args.complement,
        rebin=rebin_arg,
        clamp_yerrors=args.clamp_yerrors,
        multiply_x_values=args.multiply_x_values,
        hline=args.hline,
        debug=args.debug,
    )