2hdm-a

name

2hdm-a

description

Constraint-first interactive workflow for the 2HDM + pseudoscalar mediator (2HDM+a) DM model (arXiv:2506.19062 §III). Interviews the user for constraint selection (relic density, direct detection, indirect detection), shows the prereq chain with [EXISTS]/[PLANNED] annotations and a time estimate, gates for confirmation, then drives the hand-crafted SARAH model → /madgraph → /maddm for the relic-ready subset. Invoke when the user picks 2HDM+a from /demo or says "run 2hdm-a".

2HDM + Pseudoscalar Mediator

DEMO CAVEAT (fixture-calibrated; not cosmology-aligned): Ωh² ≈ 10.5 is ~88× the Planck value (0.12) — ~4 orders of magnitude off the physics point in arXiv:2506.19062. This demo proves SARAH → SPheno → MadDM pipeline plumbing, not physics accuracy. Expected relic band: [9.9693, 11.0187] (±5% fixture).

Per-model constraint workflow for 2HDM+a (Arcadi & Profumo arXiv:2506.19062 §III). Uses hand-crafted SARAH model (fixtures/sarah_model/); /sarah-build backport pending. All observables from real HEP tools (MadDM Boltzmann solver); no analytic shortcuts.

---

When to invoke

• User selects 2HDM+a from /demo's model picker (Step 3).
• User says "run 2hdm-a", "compute relic density for 2HDM+a", or similar.
• Do NOT invoke if the preflight check in /demo Step 0 failed (missing config.json keys or unresponsive executables).

---

Model metadata

display:
  title: "2HDM + a"
  hook:  "Pseudoscalar mediator, CP-forbidden tree SI, loop-only DD."
dm_candidates:
  - {name: chi, spin: "1/2", notes: "Dirac DM coupled via the CP-odd mediator a."}
plot_axes:
  x: {symbol: "m_a",      range: [10, 1500], units: "GeV", scale: log}
  y: {symbol: "tan_beta", range: [0.5, 50],               scale: log}
multi_component: false
time_overrides:
  dd: {cold: [3.0, 6.0]}

---

Constraints and time estimates

Prereq chains and time ranges are resolved from plugins/hep-ph-toolkit/skills/_shared/constraints.yaml (model 2hdm-a). The table below is generated by running:

python3 plugins/hep-ph-toolkit/skills/_shared/time_budget.py \
    --model 2hdm-a --constraints relic dd id

Constraint	Status	Chain	Cold (hr)	Cached (hr)
Relic density	READY (hand-crafted SARAH model)	fixture copy → /madgraph [EXISTS] → /maddm [EXISTS]	1–2	0.3–0.7
Direct detection	READY-with-caveat — loop-only (CP-forbidden tree SI) via /looptools eval; σ_SI real + EW-anchor-validated (~40%) at the benchmark, single point; box not folded into f_N (≤1.7× upward), SI only; Tier-3 smoke gated off by default, green on a tooled box	fixture copy → /madgraph [EXISTS] → /feynarts [EXISTS] → /formcalc [EXISTS] → /looptools [EXISTS] → /ddcalc [EXISTS]	3–6	0.5–1.0
Indirect detection	READY-with-caveat — parser-only; pull computation pending (v1+)	fixture copy → /madgraph [EXISTS] → /maddm [EXISTS] (generate indirect_detection) → /gamlike [v0 parser]	1–3	0.3–0.7

Note on SPheno: MadDM's internal Boltzmann solver is fed directly from the patched param_card; SPheno spectrum is not required on this path. /spheno-build was removed from the relic-density prereq chain.

All-constraints cold total (overlap-adjusted): 3.5–9.0 hr

Relic density and indirect detection are READY via the fixture path. ID is parser-only: /gamlike v0 emits gamlike/v1 JSON with ⟨σv⟩, channel fractions, and Fermi-LAT likelihood rows, but pull computation (likelihood deltas, exclusion verdicts) is deferred to a future dm-pull skill.

Direct detection is READY-with-caveat. For 2HDM+a the tree-level SI cross-section is CP-forbidden (the mediator a is CP-odd), so the physical signal is loop-only. The /looptools eval runtime now numerically evaluates the FormCalc-reduced one-loop amplitude (charged-Higgs/W box + mediator triangle) into a scattering/v1 JSON for /ddcalc, on the same hand-crafted SARAH fixture used for relic. This produces a real, EW-anchor-validated σ_SI at the benchmark — σ_SI(p)=1.18×10⁻⁴⁸, σ_SI(n)=1.21×10⁻⁴⁸ cm², with the triangle coupling matching an independent analytic EW box+triangle anchor to ~40%. It is a single benchmark point (not a scan, not an experimental exclusion claim); the pseudoscalar-quark box is not folded into f_N (≤1.7× upward if folded later) and σ_SD is null in v1. The Tier-3 FormCalc/LoopTools smoke is gated off by default (HEPPH_RUN_WOLFRAM_TESTS=1) but runs green on a tooled box, and the live /sarah-build renderer UFO backport remains separate, non-blocking debt.

---

Flow

---

Step 1 — DM-candidate declaration

Print verbatim:

For 2HDM + Pseudoscalar Mediator, the DM candidate is:

• chi — Dirac fermion, coupled to the visible sector via the CP-odd pseudoscalar mediator a.

This is a single-candidate model; relic, DD, and ID rates are computed directly for chi.

Reference Lagrangian — this is the physics content the pipeline reproduces:

-L ⊃ M_chi  chiL·chiR + h.c.                      # Dirac DM mass
   + i g_chi  a · chibar·chi                       # DM–mediator Yukawa (CP-odd, pure pseudoscalar)
   + Type-II 2HDM Yukawas (u ↔ H2; d, e ↔ H1)     # softly-broken Z2 structure
   + m12² (H1†·H2) + h.c.                          # soft Z2-breaking mass term
   + λ_P  (conj[H1]·H2) a + h.c.                   # CP-odd portal (H1–H2 ↔ a)

Fields: H1/H2 (SU(2)_L Higgs doublets, Y = +½), a (CP-odd real gauge singlet, renamed a0s in SARAH to avoid collision with the photon symbol A), chiL/chiR (Dirac DM Weyl components via the PortalDM idiom — two separate LH Weyl fields, the RH one wrapped in conj[...]). Parameters: mchi (DM Dirac mass), Ma (pseudoscalar mass), gchi (DM–mediator coupling), tan β = vu/vd, lamP (portal coupling).

Because a is CP-odd, tree-level SI DM–nucleon scattering is CP-forbidden (σ_SI_tree ≈ 0); the leading DD signal is loop-induced. Relic density proceeds through the s-channel a-resonance at m_a ≈ 2 m_chi, with the width controlled by the tan β-dependent a-SM fermion coupling.

Model source note: The SARAH model driving this demo was hand-crafted using the PortalDM idiom from SARAH/Models/SM+VL/PortalDM/ (Staub's canonical Dirac-singlet reference). It lives at plugins/hep-ph-toolkit/skills/2hdm-a/fixtures/sarah_model/. The live /sarah-build renderer backport is tracked separately; relic_approx: false (this is a genuine MadDM Boltzmann solve), model_source: "hand_crafted_sarah_model".

---

Step 2 — Constraint multi-select

Ask the user which constraints to compute:

{
  "question": "Which constraints do you want computed for this model?",
  "options": [
    {"id": "relic",    "label": "Relic density",            "description": "Ω h² via MadDM (a-resonance region) — READY"},
    {"id": "dd",       "label": "Direct detection",         "description": "Loop-only σ_SI via MadGraph + FeynArts/FormCalc + /looptools eval + DDCalc — READY-with-caveat (σ_SI real + EW-anchor-validated at the benchmark, single point; box not folded into f_N (≤1.7× upward), SI only; Tier-3 smoke deferred/gated off by default)"},
    {"id": "id",       "label": "Indirect detection (parser-only)", "description": "MadDM `generate indirect_detection` → /gamlike v0 → gamlike/v1 JSON with ⟨σv⟩, channels, and Fermi-LAT likelihood rows. Pull computation (exclusion verdicts) deferred to dm-pull v1+."},
    {"id": "collider", "label": "Collider (coming soon)",   "description": "Placeholder — execution is a no-op"}
  ],
  "allowMultiple": true,
  "required": true
}

Validation: at least one non-collider option must be selected. If the user selects collider only, re-ask with the message: "Collider is a placeholder in this iteration; nothing would run. Please also select relic, DD, or ID."

---

Step 3 — Time estimate + prereq resolve + gate

Run (or reason from) plugins/hep-ph-toolkit/skills/_shared/time_budget.py for the selected constraints and print the chain table. For the selection [relic, dd, id] the output is:

Planned chain for 2HDM + a:

  Relic density       [READY — hand-crafted SARAH model]
    fixture copy → /madgraph [EXISTS] → /maddm [EXISTS]
    cold: 1–2 hr   cached: 20–40 min
    Note: uses committed SARAH model fixture (not /sarah-build renderer).
    model_source = "hand_crafted_sarah_model"; relic_approx = false.
    MadDM's internal Boltzmann solver is fed directly from the patched param_card;
    SPheno spectrum is not required on this path.

  Direct detection    [READY-with-caveat — loop-only via /looptools eval; σ_SI real + EW-anchor-validated, single point; box not folded into f_N (≤1.7× upward), SI only; Tier-3 smoke gated off by default]
    fixture copy → /madgraph [EXISTS]
      → /feynarts [EXISTS] → /formcalc [EXISTS] → /looptools [EXISTS] → /ddcalc [EXISTS]
    cold: 3–6 hr   cached: 30–60 min
    Note: tree-level SI is CP-forbidden (σ_SI_tree ≈ 0); loop-only DD is the primary signal.
    /looptools eval evaluates the charged-Higgs/W box → scattering/v1 for /ddcalc.

  Indirect detection  [READY-with-caveat — parser-only; pull computation pending (v1+)]
    fixture copy → /madgraph [EXISTS] → /maddm [EXISTS] (generate indirect_detection)
      → /gamlike [v0 parser]
    cold: 1–3 hr   cached: 20–40 min
    Note: emits structured ⟨σv⟩ + channels + Fermi-LAT likelihood rows.
    No exclusion verdict — pull computation is the v1+ deliverable.

Overlap-adjusted totals (shared prereqs counted once):
  selected + ready : cold ~1–2 hr,   cached ~20–40 min
  selected total   : cold ~3.3–8 hr, cached ~1–2 hr  (if all prereqs existed)

Adapt the printed block to the user's actual selection (not always all three).

If any selected constraint is COMING SOON:

{
  "question": "Some selected constraints have unimplemented prereqs. How to proceed?",
  "options": [
    {"id": "run_ready", "label": "Run available (drop blocked)", "description": "Run the ready constraints; skip blocked ones with a note in the final summary."},
    {"id": "back",      "label": "Back",                          "description": "Re-select constraints."},
    {"id": "cancel",    "label": "Cancel",                        "description": "Exit cleanly."}
  ],
  "allowMultiple": false,
  "required": true
}

If all selected are READY:

{
  "question": "Run it? Total cold-run estimate: {cold_total} hr.",
  "options": [
    {"id": "go",     "label": "Run it", "description": "Execute the chain in order."},
    {"id": "back",   "label": "Back",   "description": "Re-select constraints."},
    {"id": "cancel", "label": "Cancel", "description": "Exit cleanly."}
  ],
  "allowMultiple": false,
  "required": true
}

On cancel: print "Cancelled." and exit. Do NOT write summary.json.

On back: return to Step 2.

On run_ready or go: proceed to Step 4 with only the READY subset.

---

Step 4 — Execute

Execute the READY subset of the selected constraints in order. Relic density runs end-to-end (Boltzmann solve via MadDM). Indirect detection runs as a parser-only branch (MadDM generate indirect_detection → /gamlike v0; no pull computation). DD runs the loop chain (/feynarts → /formcalc → /looptools eval → /ddcalc) to a real σ_SI when selected — see Step 4f.

---

Step 4 — Relic density branch

Overview: drive the hand-crafted SARAH model → SARAH → UFO → MadGraph5+MadDM to compute Omega h² for chi (Dirac DM) at the off-resonance benchmark (Mchi=100, Ma=400, gchi=1.0, tan β=10). At this point (Ma=400 >> 2·Mchi=200 GeV, deeply off-resonance) the dominant channels are chi chi~ → W±H∓ (~99% combined); see benchmarks/off-resonance-2hdma-004/expectations.json for per-channel fractions.

---

4a. Deploy the hand-crafted SARAH model

The working SARAH model is committed at:

plugins/hep-ph-toolkit/skills/2hdm-a/fixtures/sarah_model/

(4 files: TwoHdmAfix.m, parameters.m, particles.m, SPheno.m)

Resolve $SARAH_ROOT from config.json key sarah_path. Copy the fixture directory into $SARAH_ROOT/Models/TwoHdmAfix/ if that directory does not already exist or if any of the 4 .m files differ from the fixture.

Concretely:

SARAH_ROOT=$(python3 -c "import json; print(json.load(open('config.json'))['sarah_path'])")
DEST="$SARAH_ROOT/Models/TwoHdmAfix"
SRC="plugins/hep-ph-toolkit/skills/2hdm-a/fixtures/sarah_model"
if [ ! -d "$DEST" ] || ! diff -q "$SRC/TwoHdmAfix.m" "$DEST/TwoHdmAfix.m" > /dev/null 2>&1; then
    mkdir -p "$DEST"
    cp "$SRC"/*.m "$DEST/"
    echo "Deployed TwoHdmAfix to $DEST"
else
    echo "TwoHdmAfix already up to date at $DEST"
fi

---

4b. Run SARAH to produce the UFO

Run SARAH directly via wolframscript with the SARAH mutex lock to avoid concurrent SARAH sessions:

Cross-platform note. SARAH must run under a global mutex to prevent concurrent kernels from corrupting the shared $SARAH_ROOT/Output/<Model>/ tree. We use the repo-vendored bin/flock_run.sh (Python fcntl.flock, no system install) instead of util-linux flock(1), which Apple's BSD base does not ship. See bin/README.md for the helper's contract.

SARAH_PATH=$(jq -r '.sarah_path' "${XDG_CONFIG_HOME:-$HOME/.config}/hephaestus/config.json")
[ -n "$SARAH_PATH" ] && [ "$SARAH_PATH" != "null" ] || { echo "ERROR: SARAH_PATH not set in config.json" >&2; exit 1; }
HEPPH_STATE_ROOT="${HEPPH_STATE_ROOT:-$HOME/.local/share/hephaestus}"
SARAH_LOCK="$HEPPH_STATE_ROOT/.locks/sarah_global.lock"
REPO_ROOT="$(git rev-parse --show-toplevel 2>/dev/null || pwd)"
"$REPO_ROOT/bin/flock_run.sh" "$SARAH_LOCK" 120 -- \
  wolframscript -code "AppendTo[\$Path, \"$SARAH_PATH\"]; <<SARAH\`; Start[\"TwoHdmAfix\"]; MakeUFO[]; Quit[]"

$SARAH_PATH is read from $HOME/.config/hephaestus/config.json (XDG path, key: sarah_path) and injected into $Path before loading SARAH. SARAH runs MakeUFO[] and emits the UFO to $SARAH_ROOT/Output/TwoHdmAfix/EWSB/UFO/.

The != "null" guard on $SARAH_PATH is intentional: jq -r returns the literal string "null" when the key is absent or null — the plain [ -n "$SARAH_PATH" ] test would pass on this string.

If flock_run.sh exits 124 (timeout after 120 s), log a blocker issue with fix_owner_hint: tool_install and abort — another SARAH session may legitimately be running. Check pgrep wolframscript before retrying.

The equivalent Mathematica session (for reference):

<<SARAH`
Start["TwoHdmAfix"];
MakeUFO[];
Quit[];

Expected output directory: $SARAH_ROOT/Output/TwoHdmAfix/EWSB/UFO/

Sanity check (required before proceeding):

grep -c "chi" "$SARAH_ROOT/Output/TwoHdmAfix/EWSB/UFO/vertices.py"

Must return ≥ 3 (three chi–Ah vertices: chi chibar Ah1, chi chibar Ah2, chi chibar Ah3). If it returns 0, the model was not deployed correctly — re-run Step 4a.

The UFO uses OutputName -> "chi" for the Dirac DM particle (PDG 9989932) and OutputName -> "Ah" for the CP-odd scalars (PDG 36, 9931569, 9949515). MadDM will use chi as the DM candidate name.

---

4c. MadGraph5 + MadDM: relic density

Drive /maddm for the single off-resonance benchmark point (Mchi=100, Ma=400, gchi=1.0, tan β=10).

Two-phase invocation (required for param_card overlay). Cards/param_card.dat is written by output and consumed by launch. In between, patch the card — this is the critical step that enabled the first finite Ωh² result (iter-8 fix: PHASES[1] = 1.0 restores all DM couplings from zero).

import sys
sys.path.insert(0, "plugins/hep-ph-toolkit/skills/maddm/scripts")  # sys.path landmine fix: maddm_run lives here
from scripts.maddm_run import generate_maddm_script
from pathlib import Path
import subprocess, sys

ufo_path = f"{sarah_root}/Output/TwoHdmAfix/EWSB/UFO"
out_dir  = Path("./demo_output/2hdm-a/maddm_run/")
mg5_bin  = config["madgraph_path"] + "/bin/mg5_aMC"

# Phase 1 script: import, define DM, generate relic_density, output
setup_script = f"""\
import model {ufo_path}
define darkmatter chi
generate relic_density
output {out_dir}
exit
"""

# Phase 2 script: launch
launch_script = f"""\
launch -f
"""

(out_dir / "..").mkdir(parents=True, exist_ok=True)
(out_dir.parent / "setup.mg5").write_text(setup_script)
(out_dir.parent / "launch.mg5").write_text(launch_script)

# Phase 1: output (writes Cards/param_card.dat)
subprocess.run([mg5_bin, "--mode=maddm",
                str(out_dir.parent / "setup.mg5")],
               cwd=out_dir.parent, check=True)

# Phase 2: patch the param_card BEFORE launch
# This is iter-8's key fix: sets PHASES[1]=1.0 so pchiR≠0 and DM couplings are live.
subprocess.run([
    sys.executable,
    "plugins/hep-ph-toolkit/skills/2hdm-a/scripts/patch_paramcard.py",
    str(out_dir),
    "--Mchi", "100.0",
    "--Ma-Ah2", "400.0",   # off-resonance (Ma ≈ 4 × Mchi)
    "--gchi", "1.0",
    "--tan-beta", "10.0",
], check=True)

# Phase 3: launch (reads patched card)
subprocess.run([mg5_bin, "--mode=maddm",
                str(out_dir.parent / "launch.mg5")],
               cwd=str(out_dir), check=True)

mg5_aMC --mode=maddm is required — bare mg5_aMC <script> loads the base interpreter without the MadDM plugin and generate relic_density raises InvalidCmd.

chi is a Dirac fermion (UFO particles.py declares it with antiname = 'chibar' and self_conj = False), so MadDM uses EffOperators/COMPLEX for the EFT re-import.

Parse output. After MadDM finishes, the relic density result is in MadDM_results.txt. Use /gamlike [v0 — parser only; pull-computation v1+] (plugins/hep-ph-toolkit/skills/gamlike/scripts/parse_maddm_results.py) to parse it into a stable gamlike/v1 JSON document (see /gamlike SKILL.md for the full schema and invocation pattern):

import json, subprocess, sys
from pathlib import Path

# /gamlike v0 — call the canonical parser.
# Path committed as stable in plugins/hep-ph-toolkit/skills/gamlike/SKILL.md §"Stable invocation path"
# (D-DPATH). Demo runs from repo root, so we use that as the working dir.
maddm_results_path = out_dir / "output" / "run_01" / "MadDM_results.txt"
gamlike_json_path  = out_dir / "gamlike.json"

subprocess.run([
    sys.executable,
    "plugins/hep-ph-toolkit/skills/gamlike/scripts/parse_maddm_results.py",
    str(maddm_results_path),
    "--out", str(gamlike_json_path),
], check=True)

gamlike = json.loads(gamlike_json_path.read_text())

results = {}
results["Omegah2"] = gamlike["relic"]["Omegah2"]

# Flatten nested channels (D4): /gamlike emits {<initial>: {<final>: <pct>}}.
# 2hdm-a's relic.json schema commits to a flat {<final>: <pct>} dict for
# back-compat. We flatten here because 2HDM+a has exactly one initial state (chichibar).
flat_channels = {}
for init_state, finals in gamlike["relic"]["channels"].items():
    for k, v in finals.items():
        if v is not None:
            flat_channels[k] = v

# Convert percent values to unit fractions (D3):
# /gamlike emits raw pct; consumer normalizes.
total = sum(flat_channels.values()) or 1.0
fractions = {k: v / total for k, v in flat_channels.items()}

# Fail-fast gate: fractions must sum to [0.99, 1.01].
gate_check = {
    "channels_sum_in_unity_range": 0.99 <= sum(fractions.values()) <= 1.01,
}
if not gate_check["channels_sum_in_unity_range"]:
    raise ValueError(f"channel_fractions out of [0.99,1.01]: sum={sum(fractions.values())}")

results["channel_percentages"] = flat_channels
results["channel_fractions"]   = fractions
results["gate_check"]          = gate_check
results["sigmav_channels"]     = flat_channels  # legacy alias (raw %)

# WS2: emit summary["dm_indirect_detection_status"] = "parser-only"

See benchmarks/off-resonance-2hdma-004/expectations.json — channel fractions are parameter-point properties, not model properties.

For on-resonance points (Ma ≈ 2 × Mchi), expect Romberg integration to hang (>20 min). Debug off-resonance first; scan toward resonance only after the pipeline is confirmed working.

Record the result to ./demo_output/2hdm-a/relic.json:

{
  "m_chi":         100.0,
  "m_a":           400.0,
  "tan_beta":       10.0,
  "gchi":            1.0,
  "omega_h2":      <value from MadDM>,
  "on_resonance":  false,
  "status":        "ok" | "failed",
  "relic_approx":  false,
  "model_source":  "hand_crafted_sarah_model",
  "model_fixture": "plugins/hep-ph-toolkit/skills/2hdm-a/fixtures/sarah_model/",
  "maddm_results": "./demo_output/2hdm-a/maddm_run/output/run_01/MadDM_results.txt",
  "channel_percentages": {"bbx": <from MadDM>, "...": "..."},
  "channel_fractions":   {"bbx": <from MadDM>, "...": "..."},
  "gate_check": {
    "channels_sum_in_unity_range": true
  }
}

relic_approx: false — this is a real MadDM Boltzmann result, not an analytic approximation. model_source: "hand_crafted_sarah_model" — consumers should know the UFO came from the fixture, not the live renderer.

Observed relic density. The observed value from Planck is Omega h² ≈ 0.120 ± 0.001. The off-resonance benchmark is intentionally overabundant; the Planck band sits near the a-resonance at m_a ≈ 2 m_chi.

---

4d. Plotting

Consult /hep-plotting to produce ./demo_output/2hdm-a/summary.{pdf,png}.

Since the demo runs one benchmark point, the figure is an annihilation-channel breakdown bar chart (same pattern as singlet-doublet/SKILL.md Step 4d):

from styles.hep_ph_style import set_hep_context, check_overlaps
import matplotlib.pyplot as plt

palette = set_hep_context("analytic")

channels = sorted(results["sigmav_channels"].items(),
                  key=lambda kv: kv[1], reverse=True)[:10]
labels   = [c for c, _ in channels]
percents = [p for _, p in channels]

fig, ax = plt.subplots(figsize=(85/25.4, 63.75/25.4))
ax.barh(labels, percents, color=palette["data"])
ax.invert_yaxis()
ax.set_xlabel(r"$\langle\sigma v\rangle$ contribution [%]")
ax.text(0.02, 0.98,
        f"2HDM+a · χ m={results.get('m_chi', 100):.0f} GeV, "
        f"$m_a$={results.get('m_a', 400):.0f} GeV\n"
        f"Ω h² = {results['Omegah2']:.3f} (Planck 0.120)\n"
        "arXiv:2506.19062 §III",
        transform=ax.transAxes, va="top", ha="left",
        fontsize=7, color=palette["deemph"])

issues = check_overlaps(fig)
assert not issues, f"Overlaps detected: {issues}"

fig.tight_layout()
fig.savefig("./demo_output/2hdm-a/summary.pdf", bbox_inches="tight")
fig.savefig("./demo_output/2hdm-a/summary.png", dpi=200, bbox_inches="tight")

Key requirements per HEP convention:

• Analytic style via set_hep_context("analytic") from styles/hep_ph_style.py.
• Bars in black via palette["data"] (per project memory feedback_data_point_color).
• Figure size 85mm × 63.75mm (single-column 4:3); no legend box; left + bottom axes only.
• Save to ./demo_output/2hdm-a/summary.{pdf,png}.

---

4e. Indirect detection branch (READY-with-caveat)

Overview: drive MadDM a second time at the off-resonance benchmark with generate indirect_detection, parse via /gamlike v0 into a gamlike/v1 JSON document, and write id.json. The v0 parser surfaces ⟨σv⟩, per-channel fractions, peak/no-peak metadata, Fermi-LAT likelihood rows, and (when MadDM emits them) PPPC4DMID-style spectra. Pull computation — converting MadDM's Fermi_Likelihood(...) rows into Δχ² exclusion verdicts — is deferred to a future dm-pull skill.

import json, shutil, subprocess, sys
from pathlib import Path

id_out_dir = Path("./demo_output/2hdm-a/maddm_run_id/")
shutil.rmtree(id_out_dir, ignore_errors=True)  # MG5 `output` refuses to overwrite
id_out_dir.parent.mkdir(parents=True, exist_ok=True)

# Reuse the fixture UFO from 4b; only the MadDM target changes.
setup_id = f"""\
import model {ufo_path}
define darkmatter chi
generate indirect_detection
output {id_out_dir}
exit
"""
launch_id = "launch -f\n"
(id_out_dir.parent / "setup_id.mg5").write_text(setup_id)
(id_out_dir.parent / "launch_id.mg5").write_text(launch_id)

# Phase 1: write Cards/param_card.dat
subprocess.run([mg5_bin, "--mode=maddm", str(id_out_dir.parent / "setup_id.mg5")],
               cwd=id_out_dir.parent, check=True)

# Phase 2: patch param_card (same off-resonance benchmark + PHASES[1]=1.0 fix)
subprocess.run([
    sys.executable,
    "plugins/hep-ph-toolkit/skills/2hdm-a/scripts/patch_paramcard.py",
    str(id_out_dir),
    "--Mchi", "100.0",
    "--Ma-Ah2", "400.0",
    "--gchi", "1.0",
    "--tan-beta", "10.0",
], check=True)

# Phase 3: launch
subprocess.run([mg5_bin, "--mode=maddm", str(id_out_dir.parent / "launch_id.mg5")],
               cwd=str(id_out_dir), check=True)

# Parse MadDM ID output via /gamlike v0
maddm_id_results = id_out_dir / "output" / "run_01" / "MadDM_results.txt"
gamlike_id_json  = id_out_dir / "gamlike.json"
subprocess.run([
    sys.executable,
    "plugins/hep-ph-toolkit/skills/gamlike/scripts/parse_maddm_results.py",
    str(maddm_id_results),
    "--out", str(gamlike_id_json),
], check=True)
gamlike_id = json.loads(gamlike_id_json.read_text())

# Pull the indirect-detection blocks out of the gamlike/v1 document.
indirect       = gamlike_id.get("indirect", {})
spectral       = gamlike_id.get("spectral", {})
fluxes_source  = gamlike_id.get("fluxes_source", {})

id_doc = {
    "schema_version":               "id/v0",
    "m_chi":                        100.0,
    "m_a":                          400.0,
    "tan_beta":                     10.0,
    "gchi":                         1.0,
    "indirect":                     indirect,
    "spectral":                     spectral,
    "fluxes_source":                fluxes_source,
    "dm_indirect_detection_status": "parser-only",
    "pull_computation_pending":     True,
    "status":                       "ok" if indirect.get("present") else "failed",
    "maddm_results":                str(maddm_id_results.resolve()),
    "gamlike_json":                 str(gamlike_id_json.resolve()),
}
Path("./demo_output/2hdm-a/id.json").write_text(json.dumps(id_doc, indent=2))

Caveat banner. Every consumer of id.json MUST surface dm_indirect_detection_status: "parser-only". Do not present indirect.global.Total_xsec as an exclusion-verdict input; the Fermi-LAT and IceCube pull computations are not implemented. Conditional-emission gates G13/G16 etc. (see /gamlike SKILL §"Conditional emission gates") may suppress fields with FIELD_GATED warnings — propagate them rather than coercing to zero.

4f. Direct detection branch (READY-with-caveat — loop-only)

For 2HDM+a the tree-level SI cross-section is CP-forbidden (the mediator a is CP-odd, so σ_SI_tree ≈ 0), so loop is the only physical channel — there is no tree-DD shortcut analogous to singlet-doublet's MadDM generate direct_detection path. The loop-DD chain runs on the same hand-crafted SARAH fixture:

/madgraph → /feynarts generate → /formcalc reduce → /looptools eval → /ddcalc run --sigma-json

/feynarts generates the one-loop amplitude (charged-Higgs/W box A⁰H⁺W⁻ + mediator triangle), /formcalc reduce reduces it to FormCalc-native PV heads (amp_reduced.m + amp_reduced.meta.json), and /looptools eval numerically evaluates the PV integrals and emits a scattering/v1 JSON (sigma_si_proton_cm2, sigma_si_neutron_cm2; σ_SD null in v1):

python plugins/hep-ph-toolkit/skills/looptools/scripts/run_looptools.py eval \
    --amp-reduced <formcalc_out>/amp_reduced.m \
    --point ./demo_output/2hdm-a/maddm_run/Cards/param_card.dat \
    --output-dir ./demo_output/2hdm-a/looptools_out \
    --form-factors default_2018

Feed looptools_out/scattering.json into /ddcalc run --sigma-json for the DD comparison. The σ_SI is a real, EW-anchor-validated number at this benchmark (σ_SI(p)=1.18×10⁻⁴⁸, σ_SI(n)=1.21×10⁻⁴⁸ cm²; triangle coupling ~40% vs the independent analytic EW box+triangle anchor). It is a single benchmark point — not a scan, and not an experimental exclusion claim; the pseudoscalar-quark box is not folded into f_N (≤1.7× upward if folded later) and σ_SD is null in v1. The Tier-3 FormCalc/LoopTools smoke (HEPPH_RUN_WOLFRAM_TESTS=1) is gated off by default but runs green on a tooled box. Stamp model_source: "hand_crafted_sarah_model" and tag the DD row [loop, validated, single-point]. The live /sarah-build renderer UFO backport for the loop chain remains separate, non-blocking debt.

---

4g. Write summary.json

At the end of Step 4 (after all READY constraints have run), write ./demo_output/2hdm-a/summary.json conforming to plugins/hep-ph-toolkit/skills/2hdm-a/summary.schema.json (which $refs _shared/summary.core.schema.json).

Create the output directory first:

mkdir -p ./demo_output/2hdm-a/

Then write the file. Example for the default READY subset (relic + parser-only ID):

{
  "schema_version": "1",
  "model":         "2hdm-a",
  "run_at":        "<ISO-8601 timestamp>",
  "ran":           ["relic", "id"],
  "relic_approx":  false,
  "model_source":  "hand_crafted_sarah_model",
  "model_fixture": "plugins/hep-ph-toolkit/skills/2hdm-a/fixtures/sarah_model/",
  "dm_indirect_detection_status": "parser-only",
  "pull_computation_pending":     true,
  "skipped_constraints": [
    {"id": "dd", "reason": "not selected this run; loop-only DD is READY-with-caveat via /looptools eval (charged-Higgs/W box → scattering/v1 → /ddcalc), σ_SI real + EW-anchor-validated at the benchmark (single point), Tier-3 smoke gated off by default"}
  ],
  "artifacts_dir": "./demo_output/2hdm-a/",
  "headline": "Relic + parser-only ID at off-resonance benchmark (Mchi=100, Ma=400, gchi=1, tan β=10): Ω h² ≈ 10.494; ID artifacts via /gamlike v0 (no exclusion verdict yet). DD available via /looptools eval (loop-only; real EW-anchor-validated σ_SI at the benchmark, single point; box not folded into f_N (≤1.7× upward), SI only)."
}

Adapt ran, skipped_constraints, and headline to the actual constraints that ran.

On Cancel (user chose Cancel at Step 3 gate): do NOT write summary.json. /demo's closing block will degrade gracefully to "2hdm-a interview was cancelled.".

---

Error paths

Skill step	Error code	Trigger	Action
Step 4a deploy	FIXTURE_COPY_FAILED	sarah_path not set in config or target unwritable	Run /install to configure SARAH, then retry
Step 4b SARAH run	WOLFRAM_KERNEL_ABSENT	wolfram_engine_path not set in config	Run /install to configure Wolfram Engine, then retry
Step 4b SARAH run	SARAH_ABSENT	sarah_path not set in config	Run /install to configure SARAH, then retry
Step 4b SARAH run	VERTICES_ZERO	grep -c "chi" vertices.py returns 0	Re-deploy fixture (Step 4a); check that TwoHdmAfix.m was copied correctly
Step 4c MadDM	YUKAWA_COLLISION	mdl_ryu211 is not defined	Check that parameters.m uses YUKAWAU/YUKAWAD/YUKAWAE (not all at YUKAWA [1,1]); fixture version already has this fixed
Step 4c MadDM	Omegah2 = -1.0 sentinel	DM couplings are zero	Check patch_paramcard.py ran and set PHASES[1]=1.0; grep param_card.dat for PHASES block
Step 4c MadDM	Romberg hang >20 min	On-resonance point (Ma ≈ 2*Mchi)	Use off-resonance benchmark first (Ma=400 for Mchi=100); resonance scan requires MadDM statistics tuning
Step 4c MadDM	UFO import fails	UFO path wrong or chi PDG id unknown	Confirm SARAH ran successfully in Step 4b; check $SARAH_ROOT/Output/TwoHdmAfix/EWSB/UFO/ exists
Plotting	styles.hep_ph_style import error	styles/ not on PYTHONPATH	Ensure the repo root is on PYTHONPATH

---

File map

File	Produced by	Contents
plugins/hep-ph-toolkit/skills/2hdm-a/fixtures/sarah_model/*.m	committed fixture	Hand-crafted SARAH model (TwoHdmAfix): 4 Mathematica files
plugins/hep-ph-toolkit/skills/2hdm-a/scripts/patch_paramcard.py	committed script	param_card patcher: sets PHASES[1]=1.0, identity mixing matrices, BSM widths ≥ 1 GeV
$SARAH_ROOT/Models/TwoHdmAfix/	Step 4a	Deployed SARAH model (copied from fixture)
$SARAH_ROOT/Output/TwoHdmAfix/EWSB/UFO/	Step 4b (SARAH)	MadGraph UFO model files
./demo_output/2hdm-a/maddm_run/Cards/param_card.dat	Step 4c (MadDM output + patcher)	Patched param_card with correct PHASES, masses, widths
./demo_output/2hdm-a/maddm_run/output/run_01/MadDM_results.txt	Step 4c (MadDM launch)	Relic density + channel fractions
./demo_output/2hdm-a/relic.json	Step 4c	Omega h² + channel breakdown + metadata
./demo_output/2hdm-a/summary.png	Step 4d	Annihilation-channel bar chart
./demo_output/2hdm-a/maddm_run_id/	Step 4e	MadDM indirect-detection session output
./demo_output/2hdm-a/maddm_run_id/gamlike.json	Step 4e	Parsed gamlike/v1 JSON (⟨σv⟩, channels, Fermi-LAT likelihood rows, spectra)
./demo_output/2hdm-a/id.json	Step 4e	Per-model ID document with dm_indirect_detection_status: "parser-only" flag
./demo_output/2hdm-a/summary.json	Step 4g	Per-model summary consumed by /demo closing block