diff --git a/bin/fccanalysis b/bin/fccanalysis index ed2ff11913e..cebdc234ffb 100755 --- a/bin/fccanalysis +++ b/bin/fccanalysis @@ -16,6 +16,8 @@ from run_analysis import run from run_final_analysis import run_final from do_plots import do_plots from do_combine import do_combine +from fit import run_fit +from extract import run_extract # _____________________________________________________________________________ @@ -125,6 +127,10 @@ def main(): do_plots(parser) elif args.command == 'combine': do_combine(parser) + elif args.command == 'fit': + run_fit(parser) + elif args.command == 'extract': + run_extract(parser) else: run(parser) diff --git a/examples/fcc_ee_zh_mumu_bb.py b/examples/fcc_ee_zh_mumu_bb.py new file mode 100644 index 00000000000..7961c10f459 --- /dev/null +++ b/examples/fcc_ee_zh_mumu_bb.py @@ -0,0 +1,104 @@ +"""Example configuration script for FCC-ee ZH angular/mass fitting. + +Usage Examples: + 1. Generate the text datacard only: + fccanalysis fit examples/fcc_ee_zh_mumu_bb.py -o outputs/datacard.txt + + 2. Generate datacard and run the default AsymptoticLimits engine: + fccanalysis fit examples/fcc_ee_zh_mumu_bb.py -o outputs/datacard.txt -e + + 3. Override the default engine with custom options (e.g. FitDiagnostics): + fccanalysis fit examples/fcc_ee_zh_mumu_bb.py -o outputs/datacard.txt -e -- -M FitDiagnostics +""" + +class Fit: + def __init__(self): + # 1. Global framework configurations + self.autoMCStats = False + + # 2. Path templates for the input shape histograms + self.shapes = { + "*": { + "mumu_bjets_channel": "fcc_ee_zh_shapes.root $CHANNEL/$PROCESS $CHANNEL/$PROCESS_$SYSTEMATIC", + "mumu_inter_channel": "fcc_ee_zh_shapes.root $CHANNEL/$PROCESS $CHANNEL/$PROCESS_$SYSTEMATIC" + } + } + + # 3. Channel definitions based on event selection categories + # Category 1: High purity b-tagged category + # Category 2: Intermediate/loose b-tagged category (to constrain backgrounds) + self.channels = { + "mumu_bjets_channel": { + "observation": -1, # -1 instructs Combine to use asymptotic or toy datasets + "processes": { + "ZH_signal": {"type": "signal", "rate": -1}, + "ZZ_bkg": {"type": "background", "rate": -1}, + "Zjets_bkg": {"type": "background", "rate": -1} + } + }, + "mumu_inter_channel": { + "observation": -1, + "processes": { + "ZH_signal": {"type": "signal", "rate": -1}, + "ZZ_bkg": {"type": "background", "rate": -1}, + "Zjets_bkg": {"type": "background", "rate": -1} + } + } + } + + # 4. Systematics Matrix + # Models how uncertainties affect each process in each distinct channel + self.systematics = { + # Integrated luminosity precision at FCC-ee (highly precise, ~0.1% to 0.5%) + "lumi_FCC": { + "type": "lnN", + "apply_to": { + "ZH_signal": 1.005, + "ZZ_bkg": 1.005, + "Zjets_bkg": 1.005 + } + }, + # Muon tracking and Identification efficiency uncertainty + "muon_eff": { + "type": "lnN", + "apply_to": { + "ZH_signal": 1.01, + "ZZ_bkg": 1.01, + "Zjets_bkg": 1.01 + } + }, + # B-tagging efficiency uncertainty (highly correlated with signal) + "btag_eff": { + "type": "lnN", + "apply_to": { + "ZH_signal": {"mumu_bjets_channel": 1.03, "mumu_inter_channel": 0.98}, + "ZZ_bkg": {"mumu_bjets_channel": 1.03, "mumu_inter_channel": 0.98} + } + }, + # Mistag rate uncertainty for light/charm jets leaking into the b-jet channel + "mistag_light": { + "type": "lnN", + "apply_to": { + "Zjets_bkg": {"mumu_bjets_channel": 1.10, "mumu_inter_channel": 1.04} + } + }, + # Theoretical cross-section uncertainty for electroweak backgrounds + "theory_ZZ_xsec": { + "type": "lnN", + "apply_to": { + "ZZ_bkg": 1.04 + } + }, + # Recoil mass resolution shape systematic + "recoil_res": { + "type": "lnN", + "apply_to": { + "ZH_signal": 1.05, + "ZZ_bkg": 1.05 + } + } + } + + # Optional: Define default command-line arguments to pass to the Combine tool + # The framework will automatically append `--out ` to this. + self.combine_args = "-M FitDiagnostics" diff --git a/man/man1/fccanalysis-extract.1 b/man/man1/fccanalysis-extract.1 new file mode 100644 index 00000000000..f3d34c7984d --- /dev/null +++ b/man/man1/fccanalysis-extract.1 @@ -0,0 +1,32 @@ +.\" Manpage for fccanalysis-extract +.\" Contact FCC-PED-SoftwareAndComputing-Analysis@cern.ch to correct errors or typos. +.TH FCCANALYSIS\-EXTRACT 1 "August 2026" "0.12.0" "fccanalysis-extract man page" + +.SH NAME +\fBfccanalysis\-extract\fR \- extract parameters and uncertainties from Combine ROOT files + +.SH SYNOPSIS +.B fccanalysis extract +[\fB\-h\fR | \fB\-\-help\fR] +.B input_path + +.SH DESCRIPTION +.B fccanalysis\-extract +is a standalone utility that reads ROOT files or output directories generated by Higgs Combine. +If a directory path is provided, it automatically scans for generated output ROOT files (such as +\fIhiggsCombine*.root\fR or \fIfitDiagnostics*.root\fR) within that directory. + +It safely handles grid scans and gracefully falls back to displaying expected and +observed quantiles if the file was generated using the AsymptoticLimits method. + +.SH OPTIONS +.TP +.B input_path +Path to a specific \fIhiggsCombine*.root\fR file, or a directory containing one. +If a directory is provided, the most recently modified Combine ROOT file will be used. +.TP +.BR \-h ", " \-\-help +Prints short help message and exits\&. + +.SH SEE ALSO +fccanalysis(1), fccanalysis-fit(1) diff --git a/man/man1/fccanalysis-fit.1 b/man/man1/fccanalysis-fit.1 new file mode 100644 index 00000000000..1068057ad1a --- /dev/null +++ b/man/man1/fccanalysis-fit.1 @@ -0,0 +1,54 @@ +.\" Manpage for fccanalysis-fit +.\" Contact FCC-PED-SoftwareAndComputing-Analysis@cern.ch to correct errors or typos. +.TH FCCANALYSIS\-FIT 1 "August 2026" "0.12.0" "fccanalysis-fit man page" + +.SH NAME +\fBfccanalysis\-fit\fR \- generate combine datacards and execute statistical fits + +.SH SYNOPSIS +.B fccanalysis fit +[\fB\-h\fR | \fB\-\-help\fR] +[\fB\-o\fR \fIOUTPUT\fR | \fB\-\-output\fR \fIOUTPUT\fR] +[\fB\-b\fR \fIBACKEND\fR | \fB\-\-backend\fR \fIBACKEND\fR] +[\fB\-e\fR | \fB\-\-execute\fR] +.B analysis-script +\-\- +[\fICOMBINE ARGUMENTS\fR] + +.SH DESCRIPTION +.B fccanalysis\-fit +parses an object-oriented Python configuration script containing a \fIFit\fR class +to dynamically generate a text-based datacard for statistical inference. + +If the \fB\-\-execute\fR flag is provided, the script will automatically invoke the +underlying statistical engine (e.g., Higgs Combine) on the generated datacard. +Output ROOT files are automatically shifted to the datacard directory, and +1-sigma crossing uncertainties are automatically extracted and displayed. + +.SH OPTIONS +.TP +.B analysis-script +Path to the Python analysis script containing the \fIFit\fR class. +.TP +.BR \-h ", " \-\-help +Prints short help message and exits\&. +.TP +\fB\-o\fR \fIOUTPUT\fR, \fB\-\-output\fR \fIOUTPUT\fR +Path where the generated datacard.txt will be saved\&. +.TP +\fB\-b\fR \fIBACKEND\fR, \fB\-\-backend\fR \fIBACKEND\fR +Specify the statistical backend to use\&. Currently supported: \fBcombine\fR. +.br +Default value: \fBcombine\fR. +.TP +.BR \-e ", " \-\-execute +Execute the statistical backend after generating the datacard. +.TP +.I COMBINE ARGUMENTS +Command line arguments passed after the \fB\-\-\fR separator will not be parsed +by the framework and will be forwarded directly to the backend Combine engine +(e.g., \fB\-M FitDiagnostics\fR or \fB\-\-algo grid\fR). These override any default +arguments defined in the Python class. + +.SH SEE ALSO +fccanalysis(1), fccanalysis-script(7), fccanalysis-extract(1) diff --git a/man/man1/fccanalysis-run.1 b/man/man1/fccanalysis-run.1 index 7b3e0eb80f8..524d2f2eb92 100644 --- a/man/man1/fccanalysis-run.1 +++ b/man/man1/fccanalysis-run.1 @@ -20,6 +20,7 @@ [\fB\-\-graph\-path\fR \fIGRAPH_PATH\fR] [\fB\-\-use-data-source\fR] [\fB\-\-output\-dir\fR \fIOUTPUT_DIR\fR] +[\fB\-\-output\-format\fR {\fBttree\fR,\fBrntuple\fR}] [\fB\-a\fR \fIANALYSIS_NAME\fR | \fB\-\-analysis\-name\fR \fIANALYSIS_NAME\fR] [\fB\-\-apply\-filepath\-rewrites\fR | \fB\-\-no\-filepath\-rewrites\fR] [\fB\-s\fR \fISAMPLE_NAME\fR | \fB\-\-sample\-name\fR \fISAMPLE_NAME\fR] @@ -78,6 +79,11 @@ Load events through podio::DataSource\&. \fB\-\-output\-dir\fR \fIOUTPUT_DIR\fR Set global output directory for the analysis. .TP +\fB\-\-output\-format\fR {\fBttree\fR, \fBrntuple\fR} +Specify the output ROOT file format. +.br +Default value: \fBttree\fR. +.TP \fB\-a\fR \fIANALYSIS_NAME\fR, \fB\-\-analysis\-name\fR \fIANALYSIS_NAME\fR Set analysis name, mostly used for the output path. .TP @@ -124,7 +130,7 @@ Controls search path for the process dictionaries. The default value is \fI/cvmfs/fcc.cern.ch/FCCDicts/\fR\&. .SH SEE ALSO -fccanalysis(1), fccanalysis-script(7) +fccanalysis(1), fccanalysis-script(7), fccanalysis-fit(1), fccanalysis-extract(1) .SH BUGS Many diff --git a/man/man7/fccanalysis-script.7 b/man/man7/fccanalysis-script.7 index abbbd2d07c3..69012ecfa19 100644 --- a/man/man7/fccanalysis-script.7 +++ b/man/man7/fccanalysis-script.7 @@ -90,6 +90,22 @@ These stages do not require \fIAnalysis\fR class neither \fIbuild_graph\fR function, they have their own set of attributes, please see the examples in the \fIexamples\fR of FCCAnalyses repository or fccanalysis-final-script(7) and fccanalysis-plots-script(7) manual pages. +.TP +\fBFit style analysis\fR +The analysis script needs to contain a \fIFit\fR class of the following +structure to generate Higgs Combine datacards: +.IP +class Fit(): + def __init__(self, cmdline_args): + self.datacard_dir = "outputs/datacards/" + self.combine_args = "\-M MultiDimFit \-\-algo grid" + self.processes = { + "signal": 0, + "background": 1 + } + self.channels = { + "channel_name": { ... } + } .SH ATTRIBUTES In case of running the FCCAnalysis in the staged style user can/should add the @@ -149,6 +165,11 @@ the sub-folders of $FCCDICTSDIR\&. \fBoutput_dir\fR (mandatory) User can specify the directory for the output files\&. .TP +\fBoutput_format\fR (optional) +Specifies the output ROOT file format. Accepted values are \fBttree\fR and \fBrntuple\fR. +.br +Default value: ttree +.TP \fBanalysis_name\fR (optional) Optional name for the analysis .br @@ -232,7 +253,7 @@ sample. This section is under construction. You are invited to help :) .SH SEE ALSO -fccanalysis(1), fccanalysis-run(1) +fccanalysis(1), fccanalysis-run(1), fccanalysis-fit(1), fccanalysis-extract(1) .SH BUGS Many diff --git a/python/combine.py b/python/combine.py new file mode 100644 index 00000000000..78e82241387 --- /dev/null +++ b/python/combine.py @@ -0,0 +1,315 @@ +import os +import sys +import logging +import argparse +import importlib.util +from typing import Dict, Any + +# Use the standard FCCAnalyses logging hierarchy +LOGGER = logging.getLogger('FCCAnalyses.combine') + +class DatacardWriter: + """Writes the structured OOP configuration out to Combine's standard text format.""" + + def __init__(self, datacard_obj): + self.datacard = datacard_obj + self.channels = getattr(datacard_obj, 'channels', {}) + self.systematics = getattr(datacard_obj, 'systematics', {}) + self.shapes = getattr(datacard_obj, 'shapes', {}) + self.autoMCStats = getattr(datacard_obj, 'autoMCStats', False) + + self.metadata = self._get_metadata() + self.col_w = self._calculate_col_width() + + def _get_metadata(self) -> Dict[str, int]: + imax = len(self.channels) + jmax = max( + ( + sum(1 for p in ch.get('processes', {}).values() if p.get('type') == 'background') + for ch in self.channels.values() + ), + default=0, + ) + + kmax = len(self.systematics) + return {'imax': imax, 'jmax': jmax, 'kmax': kmax} + + def _calculate_col_width(self) -> int: + max_len = 10 + for ch_name, ch_data in self.channels.items(): + max_len = max(max_len, len(ch_name)) + for proc_name in ch_data.get('processes', {}).keys(): + max_len = max(max_len, len(proc_name)) + for syst_name in self.systematics.keys(): + max_len = max(max_len, len(syst_name)) + return max_len + 4 + + def generate(self, output_path: str): + with open(output_path, 'w') as f: + self._write_header(f) + self._write_shapes(f) + self._write_observation(f) + self._write_rates(f) + self._write_systematics(f) + self._write_automcstats(f) + LOGGER.info('Successfully generated production datacard: %s', output_path) + + def _write_header(self, f): + f.write(f"imax {self.metadata['imax']} number of channels\n") + f.write(f"jmax {self.metadata['jmax']} number of backgrounds\n") + f.write(f"kmax {self.metadata['kmax']} number of nuisance parameters\n") + f.write("-" * 50 + "\n") + + def _write_shapes(self, f): + if self.shapes: + for process, channels in self.shapes.items(): + for channel, mapping in channels.items(): + f.write(f"shapes {process} {channel} {mapping}\n") + f.write("-" * 50 + "\n") + + def _write_observation(self, f): + bin_names = "".join([f"{name:<{self.col_w}}" for name in self.channels.keys()]) + obs_values = "".join([f"{str(ch.get('observation', 0)):<{self.col_w}}" for ch in self.channels.values()]) + + f.write(f"{'bin':<{self.col_w}}{bin_names}\n") + f.write(f"{'observation':<{self.col_w}}{obs_values}\n") + f.write("-" * 50 + "\n") + + def _write_rates(self, f): + bin_line = [f"{'bin':<{self.col_w}}"] + process_name_line = [f"{'process':<{self.col_w}}"] + process_id_line = [f"{'process':<{self.col_w}}"] + rate_line = [f"{'rate':<{self.col_w}}"] + + for bin_name, channel_data in self.channels.items(): + processes = channel_data.get('processes', {}) + sorted_procs = sorted( + processes.items(), + key=lambda item: 0 if item[1].get('type') == 'signal' else 1 + ) + bkg_counter = 1 + for proc_name, proc_data in sorted_procs: + bin_line.append(f"{bin_name:<{self.col_w}}") + process_name_line.append(f"{proc_name:<{self.col_w}}") + if proc_data.get('type') == 'signal': + process_id_line.append(f"{'0':<{self.col_w}}") + else: + process_id_line.append(f"{str(bkg_counter):<{self.col_w}}") + bkg_counter += 1 + rate_line.append(f"{str(proc_data.get('rate', 0.0)):<{self.col_w}}") + + f.write("".join(bin_line) + "\n") + f.write("".join(process_name_line) + "\n") + f.write("".join(process_id_line) + "\n") + f.write("".join(rate_line) + "\n") + f.write("-" * 50 + "\n") + + def _write_systematics(self, f): + for syst_name, syst_info in self.systematics.items(): + syst_type = syst_info.get("type", "lnN") + apply_to = syst_info.get("apply_to", {}) + + # Start each row with the systematic name and its type + row_cells = [syst_name, syst_type] + + # Loop over channels + for ch_name, ch_info in self.channels.items(): + sorted_procs = sorted( + ch_info.get('processes', {}).items(), + key=lambda item: 0 if item[1].get('type') == 'signal' else 1 + ) + + for proc_name, _ in sorted_procs: + # Check if this systematic applies to the current process + if proc_name in apply_to: + val = apply_to[proc_name] + + # If it's a channel-dependent dictionary, extract the specific channel value + if isinstance(val, dict): + val = val.get(ch_name, "-") + + row_cells.append(str(val)) + else: + row_cells.append("-") + + f.write("\t".join(row_cells) + "\n") + + def _write_automcstats(self, f): + if self.autoMCStats: + f.write("-" * 50 + "\n") + f.write("* autoMCStats 0 1 1\n") + + +def sanitize_and_validate_config(user_datacard) -> None: + """Validates properties and verifies shape histogram existence in ROOT files before backend execution.""" + import ROOT + import os + + channels = getattr(user_datacard, 'channels', {}) + if not channels: + LOGGER.warning("No channels defined in the datacard object.") + + # 1. Standard structural and rate checks + for ch_name, ch_data in channels.items(): + obs = ch_data.get('observation', 0) + if not (isinstance(obs, (int, float)) and (obs >= 0 or obs == -1)): + raise ValueError(f"Sanitization Error: Invalid observation '{obs}' in channel '{ch_name}'. Must be >= 0 or -1.") + + processes = ch_data.get('processes', {}) + for proc_name, proc_data in processes.items(): + p_type = proc_data.get('type') + if p_type not in ['signal', 'background']: + raise ValueError(f"Validation Error: Invalid type '{p_type}' for process '{proc_name}' in channel '{ch_name}'. Must be 'signal' or 'background'.") + + rate = proc_data.get('rate', 0.0) + if not (isinstance(rate, (int, float)) and (rate >= 0 or rate == -1)): + raise ValueError(f"Validation Error: Invalid rate '{rate}' for process '{proc_name}'. Must be >= 0 or -1.") + + valid_syst_types = ['lnN', 'shape', 'gmM', 'gmN'] + systematics = getattr(user_datacard, 'systematics', {}) + shapes_config = getattr(user_datacard, 'shapes', {}) + + # --- PHASE A: Bulk Upfront In-Memory Pre-loading --- + targets_to_load = {} + failed_to_open_files = set() + hist_cache = {} + + # Discover all target histogram paths required across all shape systematics + for s_name, s_data in systematics.items(): + if s_data.get('type') == 'shape': + a_to = s_data.get('apply_to', {}) + for c_name, c_info in channels.items(): + for p_name in c_info.get('processes', {}).keys(): + m_rule = shapes_config.get(p_name, {}).get(c_name) or shapes_config.get('*', {}).get(c_name) + if not m_rule: + continue + rf_path = m_rule.split()[0] + if not os.path.isfile(rf_path): + continue + if p_name in a_to and str(a_to[p_name]) != '-': + b_path = m_rule.split()[1].replace('$CHANNEL', c_name).replace('$PROCESS', p_name) + for var in ['Up', 'Down']: + t_path = f"{b_path}_{s_name}{var}" + if rf_path not in targets_to_load: + targets_to_load[rf_path] = set() + targets_to_load[rf_path].add(t_path) + + # Open files once, read matching objects, decouple from disk, and close immediately + for rf_path, paths in targets_to_load.items(): + hist_cache[rf_path] = {} + r_file = ROOT.TFile.Open(rf_path, "READ") + if not r_file or r_file.IsZombie(): + failed_to_open_files.add(rf_path) + continue + for t_path in paths: + h = r_file.Get(t_path) + if h: + # Disassociate from TFile lifecycle to persist in high-speed RAM + h.SetDirectory(0) + hist_cache[rf_path][t_path] = h + else: + hist_cache[rf_path][t_path] = None + r_file.Close() + + # --- PHASE B: In-Memory Validation Engine --- + # 2. Advanced Shape Systematics Check + for syst_name, syst_data in systematics.items(): + s_type = syst_data.get('type') + if s_type not in valid_syst_types: + raise ValueError(f"Validation Error: Invalid systematic type '{s_type}' for '{syst_name}'. Allowed: {valid_syst_types}") + + # If it's a shape systematic, verify the physical histograms exist + if s_type == 'shape': + apply_to = syst_data.get('apply_to', {}) + + for ch_name, ch_info in channels.items(): + for proc_name in ch_info.get('processes', {}).keys(): + + # 1. Resolve mapping rule per process + mapping_rule = shapes_config.get(proc_name, {}).get(ch_name) or shapes_config.get('*', {}).get(ch_name) + + if not mapping_rule: + continue + + root_file_path = mapping_rule.split()[0] + + if not os.path.isfile(root_file_path): + LOGGER.warning("Shape file '%s' not created yet or using relative build path. Skipping deep object check.", root_file_path) + continue + + # 2. Check if PyROOT failed to open it + if root_file_path in failed_to_open_files: + raise ValueError(f"Validation Error: Failed to open shape ROOT file: {root_file_path}") + + # 3. Check variations + if proc_name in apply_to and str(apply_to[proc_name]) != '-': + + # Resolve the inner path template + base_path = mapping_rule.split()[1].replace('$CHANNEL', ch_name).replace('$PROCESS', proc_name) + + # Combine expects clones with "Up" and "Down" + for variation in ['Up', 'Down']: + target_hist_path = f"{base_path}_{syst_name}{variation}" + + # Retrieve the pre-loaded detached histogram object + hist = hist_cache.get(root_file_path, {}).get(target_hist_path) + + if not hist or not hist.InheritsFrom("TH1"): + raise ValueError(f"Validation Error: Target histogram '{target_hist_path}' missing in '{root_file_path}'") + + for proc_name in ch_info.get('processes', {}).keys(): + if proc_name in apply_to and str(apply_to[proc_name]) != '-': + + # Resolve the inner path template (e.g., $CHANNEL/$PROCESS -> mumu_bjets_channel/ZH_signal) + base_path = mapping_rule.split()[1].replace('$CHANNEL', ch_name).replace('$PROCESS', proc_name) + + # Combine expects clones with "Up" and "Down" appended to the nominal path/name + for variation in ['Up', 'Down']: + target_hist_path = f"{base_path}_{syst_name}{variation}" + + # Retrieve the pre-loaded detached histogram object directly from cache + hist = hist_cache.get(root_file_path, {}).get(target_hist_path) + + if not hist or not hist.InheritsFrom("TH1"): + raise ValueError( + f"Validation Error: Missing required shape histogram!\n" + f" File: {root_file_path}\n" + f" Expected Path: {target_hist_path}\n" + f" Reason: Systematic '{syst_name}' is declared as 'shape' for '{proc_name}' in '{ch_name}'." + ) + +def generate_datacard(anapath: str, output_path: str) -> None: + """Sub-command engine execution block.""" + + LOGGER.info('Loading Combine fit script from: %s', anapath) + + if not os.path.isfile(anapath): + LOGGER.error('Fit script file not found! Aborting...') + sys.exit(3) + + try: + spec = importlib.util.spec_from_file_location('user_fit', anapath) + user_module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(user_module) + except SyntaxError as err: + LOGGER.error('Syntax error encountered in the fit script:\n%s', err) + sys.exit(3) + + if not hasattr(user_module, "Fit"): + LOGGER.error('Fit script must define a class named "Fit"! Aborting...') + sys.exit(3) + + user_datacard = user_module.Fit() + + try: + sanitize_and_validate_config(user_datacard) + except ValueError as err: + LOGGER.error('%s\nAborting...', err) + sys.exit(3) + + LOGGER.info('Successfully sanitized and verified user input configuration.') + + # Run the generator matrix logic + writer = DatacardWriter(user_datacard) + writer.generate(output_path) + return getattr(user_datacard, 'combine_args', '') diff --git a/python/extract.py b/python/extract.py new file mode 100644 index 00000000000..4c7414c07e8 --- /dev/null +++ b/python/extract.py @@ -0,0 +1,131 @@ +import os +import glob +import logging +import argparse + +LOGGER = logging.getLogger('FCCAnalyses.extract') + +def find_crossing(x_vals, y_vals, target=1.0, left=True, best_fit=0.0): + """Finds the linearly interpolated x-value where y crosses the target threshold.""" + pairs = sorted(zip(x_vals, y_vals)) + if left: + filtered = [p for p in pairs if p[0] <= best_fit] + filtered.reverse() + else: + filtered = [p for p in pairs if p[0] >= best_fit] + + for i in range(len(filtered) - 1): + x1, y1 = filtered[i] + x2, y2 = filtered[i+1] + + if (y1 <= target <= y2) or (y2 <= target <= y1): + if y2 == y1: return x1 + return x1 + (target - y1) * (x2 - x1) / (y2 - y1) + return None + +def extract_fit_results(target_path): + """Extracts uncertainties from a specific Combine ROOT file or directory.""" + try: + import ROOT + ROOT.gROOT.SetBatch(True) + ROOT.PyConfig.IgnoreCommandLineOptions = True + except ImportError: + LOGGER.error("PyROOT is not available. Cannot extract uncertainties.") + return + + if os.path.isdir(target_path): + root_files = glob.glob(os.path.join(target_path, "higgsCombine*.root")) + \ + glob.glob(os.path.join(target_path, "fitDiagnostics*.root")) + + if not root_files: + LOGGER.warning("No higgsCombine output ROOT file found in directory: %s", target_path) + return + target_file = max(root_files, key=os.path.getmtime) + else: + target_file = target_path + + if not os.path.exists(target_file): + LOGGER.error("Target file does not exist: %s", target_file) + return + + LOGGER.info("Extracting fit results from: %s", target_file) + f_in = ROOT.TFile.Open(target_file, "READ") + if not f_in or f_in.IsZombie(): + LOGGER.error("Could not open output ROOT file: %s", target_file) + return + + tree = f_in.Get("limit") + if not tree: + LOGGER.error("Could not find 'limit' TTree inside %s", target_file) + f_in.Close() + return + + tree.GetEntry(0) + + if hasattr(tree, "deltaNLL"): + x_vals, y_vals = [], [] + + ignore_branches = {"limit", "limitErr", "syst", "iTheta", "iEta", "quantileExpected", "deltaNLL", "t_cpu", "t_real"} + branch_names = [b.GetName() for b in tree.GetListOfBranches()] + candidates = [b for b in branch_names if b not in ignore_branches] + + if "r" in candidates: + poi_name = "r" + elif candidates: + poi_name = next((c for c in candidates if c != "mh"), candidates[0]) + else: + poi_name = "r" + + for i in range(tree.GetEntries()): + tree.GetEntry(i) + # Relaxed the deltaNLL filter to 1000 so steep curves aren't deleted + if tree.quantileExpected < -1.5 or tree.deltaNLL > 1000: continue + + val = getattr(tree, poi_name, None) + if val is not None: + x_vals.append(val) + y_vals.append(2.0 * tree.deltaNLL) + + if len(x_vals) > 1: + sorted_pairs = sorted(zip(x_vals, y_vals)) + x_vals, y_vals = zip(*sorted_pairs) + min_idx = y_vals.index(min(y_vals)) + best_fit = x_vals[min_idx] + + unc_left = find_crossing(x_vals, y_vals, target=1.0, left=True, best_fit=best_fit) + unc_right = find_crossing(x_vals, y_vals, target=1.0, left=False, best_fit=best_fit) + + err_down = (unc_left - best_fit) if unc_left is not None else float('nan') + err_up = (unc_right - best_fit) if unc_right is not None else float('nan') + + print("\n" + "="*50) + print(f" FIT RESULTS: MultiDimFit ({poi_name})") + print("="*50) + print(f" Best-fit {poi_name} : {best_fit:.6f}") + print(f" -1 Sigma (68% CL): {err_down:+.6f} (Value: {unc_left:.6f})" if unc_left else " -1 Sigma (68% CL): N/A") + print(f" +1 Sigma (68% CL): {err_up:+.6f} (Value: {unc_right:.6f})" if unc_right else " +1 Sigma (68% CL): N/A") + print("="*50 + "\n") + else: + LOGGER.warning("Tree contains only 1 point. Run Combine with a scan (e.g., --algo grid) to extract crossing uncertainties.") + + else: + limits = {} + for i in range(tree.GetEntries()): + tree.GetEntry(i) + limits[round(tree.quantileExpected, 3)] = tree.limit + + print("\n" + "="*50) + print(" FIT RESULTS: AsymptoticLimits (95% CL)") + print("="*50) + if -1.0 in limits: print(f" Observed Limit : r < {limits[-1.0]:.6f}") + if 0.5 in limits: print(f" Expected Median : r < {limits[0.5]:.6f}") + if 0.16 in limits and 0.84 in limits: print(f" Expected 68% Band: {limits[0.16]:.6f} - {limits[0.84]:.6f}") + if 0.025 in limits and 0.975 in limits: print(f" Expected 95% Band: {limits[0.025]:.6f} - {limits[0.975]:.6f}") + print("="*50 + "\n") + + f_in.Close() + +def run_extract(parser: argparse.ArgumentParser) -> None: + """Sub-command entry point for extracting uncertainties.""" + args, _ = parser.parse_known_args() + extract_fit_results(args.input_path) diff --git a/python/fit.py b/python/fit.py new file mode 100644 index 00000000000..99ca07ba2a9 --- /dev/null +++ b/python/fit.py @@ -0,0 +1,92 @@ +import os +import sys +import logging +import argparse +import subprocess +import shutil +import glob +import shlex +from extract import extract_fit_results + +LOGGER = logging.getLogger('FCCAnalyses.fit') + +def run_fit(parser: argparse.ArgumentParser) -> None: + """Sub-command entry point for object-oriented fitting configurations.""" + + args, tool_args = parser.parse_known_args() + + anapath = os.path.abspath(args.script_path) + output_path = args.output + backend = args.backend.lower() + + LOGGER.info('Steering fit configuration towards the "%s" backend...', backend) + + output_dir = os.path.dirname(os.path.abspath(output_path)) + os.makedirs(output_dir, exist_ok=True) + + if backend == 'combine': + from combine import generate_datacard + + class_combine_args = generate_datacard(anapath, output_path) + + if args.execute: + if shutil.which('combine') is None: + LOGGER.error('The "combine" command-line tool cannot be found...') + sys.exit(6) + + LOGGER.info('Launching Combine statistical engine execution on: %s', output_path) + try: + if tool_args and tool_args[0] == '--': + tool_args = tool_args[1:] + + cleaned_args = [arg for arg in tool_args if arg != output_path] + class_args = shlex.split(class_combine_args) if class_combine_args else [] + + # Strip class-level method flags if CLI explicitly passes a method override + has_cli_method = '-M' in cleaned_args or '--method' in cleaned_args + if has_cli_method: + for flag in ['-M', '--method']: + while flag in class_args: + idx = class_args.index(flag) + del class_args[idx:idx+2] + + combined_args = class_args + cleaned_args + + if not any('MINIMIZER_analytic' in arg for arg in combined_args): + combined_args = ['--X-rtd', 'MINIMIZER_analytic'] + combined_args + + # Build the final command vector + if '-M' in combined_args or '--method' in combined_args: + full_command = ['combine'] + combined_args + [output_path] + else: + full_command = ['combine', '-M', 'MultiDimFit'] + combined_args + [output_path] + + LOGGER.info("Executing command: %s", " ".join(full_command)) + subprocess.run(full_command, check=True) + + # --- ROBUST ARTIFACT SHIFTING --- + # Scoop up any ROOT files Combine dropped in the working directory + artifacts = glob.glob("higgsCombine*.root") + glob.glob("fitDiagnostics*.root") + for artifact in artifacts: + src_file = os.path.abspath(artifact) + dest_file = os.path.abspath(os.path.join(output_dir, os.path.basename(artifact))) + + # Only move/replace if source and destination paths differ + if src_file != dest_file: + if os.path.exists(dest_file): + os.remove(dest_file) + shutil.move(src_file, output_dir) + + # Automatically extract and display results from the target output directory + extract_fit_results(output_dir) + + except subprocess.CalledProcessError: + LOGGER.error('Combine statistical fitting execution failed!') + sys.exit(7) + + except KeyboardInterrupt: + LOGGER.info('Fit execution interrupted by user (Ctrl+C). Terminating cleanly...') + sys.exit(0) + else: + LOGGER.error('Backend "%s" is not implemented yet. Supported backends: combine.', backend) + sys.exit(4) \ No newline at end of file diff --git a/python/parsers.py b/python/parsers.py index ad0c4041bf7..f97198e978f 100644 --- a/python/parsers.py +++ b/python/parsers.py @@ -301,6 +301,18 @@ def setup_run_parser_combine(parser): ''' parser.add_argument('script_path', help="path to the combine script") +# _____________________________________________________________________________ +def setup_run_parser_fit(parser): + ''' + Define command line arguments for the fit sub-command. + ''' + parser.add_argument('script_path', help="path to the object-oriented fit script") + parser.add_argument('-o', '--output', type=str, default="generated_datacard.txt", + help="path to save the output text datacard") + parser.add_argument('-b', '--backend', type=str, choices=['combine', 'pyhf'], default='combine', + help="fitting backend infrastructure to target") + parser.add_argument('-e', '--execute', action='store_true', + help="automatically execute the fit calculation backend after generating the model") # _____________________________________________________________________________ def setup_subparsers(topparser): @@ -333,6 +345,19 @@ def setup_subparsers(topparser): parser_run_combine = topparser.add_parser( 'combine', help="prepare combine cards to run basic template fits") + parser_run_fit = topparser.add_parser( + 'fit', + help="generate combine datacards using object-oriented Python config." + "Use '-- [args]' to forward custom flags directly to the backend engine.") + parser_extract = topparser.add_parser( + 'extract', + help='Extract fit parameters and uncertainties directly from Combine ROOT files' + ) + parser_extract.add_argument( + 'input_path', + type=str, + help='Path to a higgsCombine ROOT file, or a directory containing one.' + ) # Register sub-parsers setup_build_parser(parser_build) @@ -343,3 +368,4 @@ def setup_subparsers(topparser): setup_run_parser_final(parser_run_final) setup_run_parser_plots(parser_run_plots) setup_run_parser_combine(parser_run_combine) + setup_run_parser_fit(parser_run_fit) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d72ee367367..097976b2ea0 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -21,3 +21,16 @@ if(WITH_PODIO_DATASOURCE) add_integration_test("examples/data_source/stages_source.py") add_integration_test("examples/data_source/analysis_stage1.py") endif() + +# Python Integration Tests for Combine Backend (Fit & Extract) +find_program(PYTHON_EXECUTABLE NAMES python3 python) + +if(PYTHON_EXECUTABLE) + add_test( + NAME test_combine_backend + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/python/test_combine_backend.py + ) + + # Ensure test data assets are fully downloaded via CI before running this test + set_tests_properties(test_combine_backend PROPERTIES WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}) +endif() diff --git a/tests/python/test_combine_backend.py b/tests/python/test_combine_backend.py new file mode 100644 index 00000000000..5aa776775c3 --- /dev/null +++ b/tests/python/test_combine_backend.py @@ -0,0 +1,82 @@ +import os +import subprocess +import shutil +import sys + +# Dynamically calculate the repository root directory +SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) +REPO_ROOT = os.path.abspath(os.path.join(SCRIPT_DIR, "..", "..")) + +def test_combine_backend_execution(): + """Integration test to verify datacard generation, fit execution, and result extraction.""" + + # Anchor all relative tracks directly to the dynamic REPO_ROOT + test_config = os.path.join(REPO_ROOT, "examples", "fcc_ee_zh_mumu_bb.py") + tmp_root = os.environ.get("TMPDIR") or os.environ.get("TEMP") or "/tmp" + test_output_dir = os.path.join(tmp_root, f"fccanalyses_test_integration_mumu_{os.getpid()}") + test_datacard = os.path.join(test_output_dir, "datacard.txt") + + if not os.path.exists(test_config): + print(f"----> ERROR: Target config missing at: {test_config}") + sys.exit(1) + + if os.path.exists(test_output_dir): + shutil.rmtree(test_output_dir) + + # Target the local script asset directly to avoid pulling the global CVMFS version + local_fccanalysis = os.path.join(REPO_ROOT, "bin", "fccanalysis") + + # Explicitly preserve and inject local paths into PYTHONPATH for isolated runners + test_env = os.environ.copy() + local_python_dir = os.path.join(REPO_ROOT, "python") + local_bin_dir = os.path.join(REPO_ROOT, "bin") + current_python_path = test_env.get("PYTHONPATH", "") + + if current_python_path: + test_env["PYTHONPATH"] = f"{local_python_dir}:{local_bin_dir}:{current_python_path}" + else: + test_env["PYTHONPATH"] = f"{local_python_dir}:{local_bin_dir}" + + # Test `fccanalysis fit` + command_fit = [ + sys.executable, local_fccanalysis, "fit", test_config, + "-o", test_datacard, + "-e", + "--", "-M", "FitDiagnostics" + ] + + print(f"----> Running fit verification command: {' '.join(command_fit)}") + result_fit = subprocess.run(command_fit, capture_output=True, text=True, cwd=REPO_ROOT, env=test_env) + + if result_fit.returncode == 6: + print("----> WARNING: 'combine' tool not found in this environment. Skipping fit & extraction execution test.") + sys.exit(0) + elif result_fit.returncode != 0: + print(f"----> ERROR: Framework fit execution failed!\nSTDOUT:\n{result_fit.stdout}\nSTDERR:\n{result_fit.stderr}") + sys.exit(1) + + if not os.path.exists(test_datacard): + print("----> ERROR: Datacard file asset was not generated successfully!") + sys.exit(1) + + print("----> INFO: Datacard generation & fit execution PASSED.") + + # Test `fccanalysis extract` + command_extract = [ + sys.executable, local_fccanalysis, "extract", test_output_dir + ] + + print(f"----> Running extraction verification command: {' '.join(command_extract)}") + result_extract = subprocess.run(command_extract, capture_output=True, text=True, cwd=REPO_ROOT, env=test_env) + + if result_extract.returncode != 0: + print(f"----> ERROR: Framework extraction execution failed!\nSTDOUT:\n{result_extract.stdout}\nSTDERR:\n{result_extract.stderr}") + sys.exit(1) + + print("----> INFO: Result extraction PASSED.") + +if __name__ == "__main__": + print("----> INFO: Starting Combine backend integration test suite (fit + extract)...") + test_combine_backend_execution() + print("----> INFO: Integration test suite PASSED successfully!") + sys.exit(0) \ No newline at end of file