#!/usr/bin/env python3

"""

A Python script, doing a deep STEP/XCAF -> ROOT TGeo conversion.

For now, all CAD solids are simply meshed. The ROOT geometry is build as a C++ ROOT macro

and facet data is stored in binary form to keep disc space minimal.

NEW (03/2026):

- Optional material/medium emission from a BOM (bill of materials) CSV file.

The CSV is expected to contain lines like:

CAD, Mechanical/Part, <PartNumber>, <Rev>, <Name>, <Mass>, <Material>, ...

- If both a part mass and a CAD volume are available, an effective density is computed

and used in the emitted TGeoMaterial. Otherwise a reasonable default density is used

for a few common materials, or 1.0 g/cm^3 as fallback.

Generates (into --output-folder):

- geom.C (small ROOT macro)

- facets_<VOLNAME>_<LID>.bin for each leaf logical volume (float32 triangles)

Facet file format (little-endian):

uint32 nTriangles

then nTriangles * 9 * float32:

ax ay az bx by bz cx cy cz

VOLNAME is a filename-safe version of the XCAF label name when available (e.g. "nut"),

and LID is the XCAF label entry (e.g. "0:1:1:7" -> "0_1_1_7") to keep filenames unique.

Naming:

- C++ variable names stay based on XCAF label entry (e.g. 0:1:1:7) for uniqueness.

- ROOT object names (TGeoVolume / TGeoTessellated / TGeoVolumeAssembly) use the label's

human name when available (e.g. "nut", "rod-assembly"), falling back to the entry.

Units:

- By default, the script tries to detect the STEP LENGTH unit by scanning the STEP file

header/contents (common patterns like .MILLI. / .CENTI. / .METRE. / INCH / FOOT).

- You can override with --step-unit {auto,mm,cm,m,in,ft}. TGeo expects cm.

Author:

- Sandro Wenzel, CERN (02/2026)

- Material/BOM integration patch (03/2026)

"""

import warnings

warnings.filterwarnings("ignore", message=".*all to deprecated function.*", category=DeprecationWarning)

import argparse

import csv

import json

import math

import re

import struct

from dataclasses import dataclass

from pathlib import Path as _Path

from typing import Dict, List, Optional, Tuple

from OCC.Core.Bnd import Bnd_Box

from OCC.Core.BRepBndLib import brepbndlib

from OCC.Core.BRepMesh import BRepMesh_IncrementalMesh

from OCC.Core.BRep import BRep_Tool

from OCC.Core.TopLoc import TopLoc_Location

from OCC.Core.TopAbs import TopAbs_REVERSED

from OCC.Extend.TopologyUtils import TopologyExplorer

from OCC.Core.STEPCAFControl import STEPCAFControl_Reader

from OCC.Core.TDocStd import TDocStd_Document

from OCC.Core.XCAFDoc import XCAFDoc_DocumentTool

from OCC.Core.IFSelect import IFSelect_RetDone

from OCC.Core.TDF import TDF_Label, TDF_LabelSequence, TDF_Tool

from OCC.Core.TCollection import TCollection_AsciiString

from OCC.Core.gp import gp_Trsf

# volume properties for density calcs (may not be present in older pythonOCC builds)

try:

from OCC.Core.GProp import GProp_GProps

from OCC.Core.BRepGProp import brepgprop_VolumeProperties

_HAS_VOLPROPS = True

except Exception:

_HAS_VOLPROPS = False

# -------------------------------

# STEP/XCAF loading

# -------------------------------

def load_step_with_xcaf(path: str):

doc = TDocStd_Document("pythonocc-doc")

reader = STEPCAFControl_Reader()

reader.SetColorMode(True)

reader.SetNameMode(True)

reader.SetLayerMode(True)

status = reader.ReadFile(path)

if status != IFSelect_RetDone:

raise RuntimeError(f"STEP read failed for: {path}")

reader.Transfer(doc)

shape_tool = XCAFDoc_DocumentTool.ShapeTool(doc.Main())

return doc, shape_tool

def label_id(label: TDF_Label) -> str:

s = TCollection_AsciiString()

TDF_Tool.Entry(label, s)

return s.ToCString()

def label_name(label: TDF_Label) -> str:

# Uses the XCAF/STEP name when present; can be empty.

try:

n = label.GetLabelName()

if n:

return str(n)

except Exception:

pass

return ""

# -------------------------------

# Units

# -------------------------------

def step_unit_scale_to_cm(step_unit: str) -> float:

step_unit = (step_unit or "auto").lower()

if step_unit == "mm":

return 0.1

if step_unit == "cm":

return 1.0

if step_unit == "m":

return 100.0

if step_unit == "in":

return 2.54

if step_unit == "ft":

return 30.48

raise ValueError(f"Unknown --step-unit {step_unit} (use auto, mm, cm, m, in, ft)")

def detect_step_length_unit(step_path: str) -> str:

"""

Heuristic unit detection by scanning STEP file text for common unit tokens.

This avoids relying on OCCT APIs that can vary across pythonOCC builds.

Returns one of: mm, cm, m, in, ft. Defaults to mm if uncertain.

"""

p = _Path(step_path)

# STEP can be huge: read only the first few MB; units are near the header.

max_bytes = 4 * 1024 * 1024

data = p.open("rb").read(max_bytes).decode("latin-1", errors="ignore").upper()

if ".MILLI." in data:

return "mm"

if ".CENTI." in data:

return "cm"

if ".METRE." in data or ".METER." in data:

return "m"

if "INCH" in data:

return "in"

if "FOOT" in data or "FEET" in data:

return "ft"

# Conservative default for mechanical CAD STEP is mm

return "mm"

# -------------------------------

# Triangulation helpers

# -------------------------------

def _scale_triangles(triangles, s: float):

if s == 1.0:

return triangles

out = []

for (a, b, c) in triangles:

out.append((

(a[0] * s, a[1] * s, a[2] * s),

(b[0] * s, b[1] * s, b[2] * s),

(c[0] * s, c[1] * s, c[2] * s),

))

return out

def triangulate_asbbox(shape, scale_to_cm: float = 1.0):

box = Bnd_Box()

brepbndlib.Add(shape, box)

xmin, ymin, zmin, xmax, ymax, zmax = box.Get()

p000 = (xmin, ymin, zmin)

p001 = (xmin, ymin, zmax)

p010 = (xmin, ymax, zmin)

p011 = (xmin, ymax, zmax)

p100 = (xmax, ymin, zmin)

p101 = (xmax, ymin, zmax)

p110 = (xmax, ymax, zmin)

p111 = (xmax, ymax, zmax)

triangles = [

(p000, p100, p110), (p000, p110, p010),

(p001, p111, p101), (p001, p011, p111),

(p000, p101, p100), (p000, p001, p101),

(p010, p110, p111), (p010, p111, p011),

(p000, p010, p011), (p000, p011, p001),

(p100, p101, p111), (p100, p111, p110),

]

return _scale_triangles(triangles, scale_to_cm)

def triangulate_CAD_solid(my_solid, meshparam, scale_to_cm: float = 1.0):

lin_defl = float(meshparam.get("lin_defl", 0.1))

ang_defl = float(meshparam.get("ang_defl", 0.1))

parallel = True

try:

BRepMesh_IncrementalMesh(my_solid, lin_defl, False, ang_defl, bool(parallel))

except TypeError:

BRepMesh_IncrementalMesh(my_solid, lin_defl, False, ang_defl)

triangles = []

for face in TopologyExplorer(my_solid).faces():

loc = TopLoc_Location()

triangulation = BRep_Tool.Triangulation(face, loc)

if triangulation is None:

continue

trsf = loc.Transformation()

reverse = (face.Orientation() == TopAbs_REVERSED)

for i in range(1, triangulation.NbTriangles() + 1):

tri = triangulation.Triangle(i)

n1, n2, n3 = tri.Get()

p1 = triangulation.Node(n1).Transformed(trsf)

p2 = triangulation.Node(n2).Transformed(trsf)

p3 = triangulation.Node(n3).Transformed(trsf)

if reverse:

p2, p3 = p3, p2

triangles.append((

(p1.X(), p1.Y(), p1.Z()),

(p2.X(), p2.Y(), p2.Z()),

(p3.X(), p3.Y(), p3.Z()),

))

return _scale_triangles(triangles, scale_to_cm)

# -------------------------------

# Volume helpers (for density)

# -------------------------------

def volume_cm3_of_shape(shape, scale_to_cm: float) -> float:

"""Compute CAD solid volume in cm^3 (using STEP->cm scale)."""

if _HAS_VOLPROPS:

try:

props = GProp_GProps()

brepgprop_VolumeProperties(shape, props)

# volume returned in STEP length units^3

v = float(props.Mass())

return v * (scale_to_cm ** 3)

except Exception:

pass

# Fallback: bounding-box volume (rough but always defined)

box = Bnd_Box()

brepbndlib.Add(shape, box)

xmin, ymin, zmin, xmax, ymax, zmax = box.Get()

dx, dy, dz = (xmax - xmin) * scale_to_cm, (ymax - ymin) * scale_to_cm, (zmax - zmin) * scale_to_cm

return max(dx, 0.0) * max(dy, 0.0) * max(dz, 0.0)

# -------------------------------

# Naming helpers

# -------------------------------

def sanitize_cpp_name(s: str) -> str:

safe = re.sub(r"[^0-9a-zA-Z]", "_", s)

if not safe:

safe = "x"

if not (safe[0].isalpha() or safe[0] == "_"):

safe = "_" + safe

return safe

def sanitize_filename(s: str) -> str:

safe = re.sub(r"[^0-9a-zA-Z]", "_", s)

return safe or "x"

# -------------------------------

# Binary facet IO

# -------------------------------

def write_facets_bin(path: _Path, triangles):

path.parent.mkdir(parents=True, exist_ok=True)

with open(path, "wb") as f:

f.write(struct.pack("<I", len(triangles)))

for (a, b, c) in triangles:

f.write(struct.pack(

"<9f",

float(a[0]), float(a[1]), float(a[2]),

float(b[0]), float(b[1]), float(b[2]),

float(c[0]), float(c[1]), float(c[2]),

))

# -------------------------------

# BOM / material mapping

# -------------------------------

@dataclass(frozen=True)

class BomEntry:

part_number: str

revision: str

name: str

mass_value: float # as in CSV

material: str

@property

def part_number_key(self) -> str:

return (self.part_number or "").strip()

@property

def name_key(self) -> str:

return (self.name or "").strip()

def _to_float(s: str) -> Optional[float]:

try:

if s is None:

return None

s = str(s).strip()

if not s:

return None

return float(s)

except Exception:

return None

def read_bom_csv(csv_path: str) -> List[BomEntry]:

"""

Reads a BOM CSV in the format provided by design team.

We look for rows whose first column is 'CAD' and second is 'Mechanical/Part'.

Columns (0-based):

0 CAD

1 type

2 part number

3 revision

4 name/description

5 mass

6 material

"""

entries: List[BomEntry] = []

with open(csv_path, newline="", encoding="utf-8", errors="ignore") as f:

reader = csv.reader(f)

for row in reader:

if not row:

continue

if len(row) < 7:

continue

if row[0].strip() != "CAD":

continue

if row[1].strip() != "Mechanical/Part":

continue

part_no = (row[2] or "").strip()

rev = (row[3] or "").strip()

name = (row[4] or "").strip()

mass = _to_float(row[5])

mat = (row[6] or "").strip()

if not (part_no or name):

continue

if mass is None:

mass = float("nan")

if not mat:

mat = "Default"

entries.append(BomEntry(part_no, rev, name, float(mass), mat))

return entries

def normalize_material_name(mat: str) -> str:

"""

Normalizes a BOM material string for matching / caching.

Note: We keep the *original* string for ROOT object names; this is only used

internally for robust matching and dictionary keys.

"""

mat = (mat or "Default").strip()

mat = re.sub(r"\s+", " ", mat)

return mat

def _norm_tokens(s: str) -> List[str]:

s = (s or "").lower()

# common grade/format noise

s = re.sub(r"\(.*?\)", " ", s)

s = s.replace("en aw", " ")

s = s.replace("en-aw", " ")

s = s.replace("en", " ")

s = s.replace("aw", " ")

s = s.replace("_", " ").replace("-", " ")

s = re.sub(r"[^a-z0-9]+", " ", s)

s = re.sub(r"\s+", " ", s).strip()

if not s:

return []

toks = s.split(" ")

# small synonym normalization

syn = {

"alu": "al",

"aluminium": "aluminum",

"silicium": "silicon",

"inox": "stainless",

"ss": "stainless",

"cu": "copper",

"fe": "iron",

"ptfe": "teflon",

"ti": "titanium",

"be": "beryllium",

}

# Expand common element symbols to names and vice-versa so that e.g. "G4_Si" can match "silicon".

elem_alias = {

"h": "hydrogen", "he": "helium", "c": "carbon", "n": "nitrogen", "o": "oxygen",

"al": "aluminum", "si": "silicon", "fe": "iron", "cu": "copper", "be": "beryllium",

"mg": "magnesium", "mn": "manganese", "cr": "chromium", "ni": "nickel", "zn": "zinc",

"ti": "titanium", "w": "tungsten", "pb": "lead", "sn": "tin",

}

name_to_sym = {v: k for k, v in elem_alias.items()}

out: List[str] = []

for t in toks:

t2 = syn.get(t, t)

out.append(t2)

if t2 in elem_alias:

out.append(elem_alias[t2])

if t2 in name_to_sym:

out.append(name_to_sym[t2])

# de-dup while preserving order

seen = set()

out2: List[str] = []

for t in out:

if t and t not in seen:

seen.add(t)

out2.append(t)

return out2

def _density_score(rho_part: Optional[float], rho_ref: Optional[float]) -> float:

if rho_part is None or rho_ref is None or not (rho_part > 0.0) or not (rho_ref > 0.0):

return 0.0

# symmetric score in log-space; 1.0 is perfect match

d = abs(math.log(rho_ref / rho_part))

return 1.0 / (1.0 + d)

def _token_score(tokens_a: List[str], tokens_b: List[str]) -> float:

if not tokens_a or not tokens_b:

return 0.0

sa = set(tokens_a)

sb = set(tokens_b)

inter = len(sa & sb)

union = len(sa | sb)

if union == 0:

return 0.0

return inter / union

def load_g4_nist_db(json_path: str) -> Dict[str, dict]:

"""

Loads a JSON dump created by the 'nist_export_all' tool.

Returns a dict: nist_name -> material record.

"""

with open(json_path, "r", encoding="utf-8") as f:

data = json.load(f)

mats = data.get("materials", {})

if not isinstance(mats, dict) or not mats:

raise RuntimeError(f"G4 NIST DB JSON seems empty or malformed: {json_path}")

return mats

# Minimal periodic table for parsing custom alloys not present in NIST.

# Values: Z (atomic number), A (g/mol)

_ELEMENT_TABLE = {

"H": (1, 1.00794),

"C": (6, 12.0107),

"N": (7, 14.0067),

"O": (8, 15.9994),

"Al": (13, 26.9815385),

"Si": (14, 28.0855),

"Fe": (26, 55.845),

"Cu": (29, 63.546),

"Be": (4, 9.0121831),

"Mg": (12, 24.305),

"Mn": (25, 54.938044),

"Cr": (24, 51.9961),

"Ni": (28, 58.6934),

"Zn": (30, 65.38),

"Ti": (22, 47.867),

"W": (74, 183.84),

"Pb": (82, 207.2),

"Sn": (50, 118.71),

}

@dataclass

class ResolvedMaterial:

bom_name: str

nist_name: Optional[str] # e.g. "G4_Al"

score: float

rho_used_g_cm3: Optional[float] # density used in ROOT definition

radlen_cm: Optional[float]

intlen_cm: Optional[float]

elements: Optional[List[dict]] # list of {symbol,Z,A_g_mol,mass_fraction}

note: str # for comments in geom.C (warnings/FIXME)

@dataclass

class MatMatchConfig:

# Minimum combined score to accept a match.

min_score: float = 0.35

# If (best - second_best) < ambiguity_delta, treat as ambiguous/unresolved.

ambiguity_delta: float = 0.05

# Weights for the combined score = w_token * token_score + w_density * density_score

w_token: float = 0.75

w_density: float = 0.25

# Optional hard filter on density proximity (in log-space). If <=0, disabled.

# Example: max_log_density_diff=0.8 means accept within exp(0.8)~2.2x in either direction.

max_log_density_diff: float = 0.0

# Penalize compound matches (oxide/dioxide/carbide/...) when BOM doesn't mention those tokens.

compound_penalty: float = 0.25

def resolve_bom_material(

bom_material: str,

rho_part_g_cm3: Optional[float],

g4db: Optional[Dict[str, dict]],

cfg: MatMatchConfig,

) -> ResolvedMaterial:

"""

Resolves an arbitrary BOM material string to a Geant4 NIST material name using:

- exact key match (BOM already uses e.g. "G4_Al")

- token overlap scoring on names

- density proximity scoring (if rho_part_g_cm3 available)

If unresolved/ambiguous, tries to parse element symbols from the BOM string (e.g. "Cu Be")

and emits a placeholder mixture (equal mass fractions) annotated with FIXME.

"""

raw_bom_material = (bom_material or "").strip()

bom_material = normalize_material_name(bom_material)

if not g4db:

return ResolvedMaterial(

bom_name=bom_material,

nist_name=None,

score=0.0,

rho_used_g_cm3=rho_part_g_cm3,

radlen_cm=None,

intlen_cm=None,

elements=None,

note="FIXME: No Geant4 NIST DB provided; using dummy material.",

)

# Trivial: BOM already provides an exact Geant4 material key

if bom_material in g4db:

rec = g4db[bom_material]

rho_ref = rec.get("density_g_cm3")

# Use NIST density for emission; CAD-derived density is used only for matching.

rho_used = rho_ref

rad = rec.get("radlen_cm")

itl = rec.get("intlen_cm")

return ResolvedMaterial(

bom_name=bom_material,

nist_name=bom_material,

score=1.0,

rho_used_g_cm3=rho_used,

radlen_cm=rad,

intlen_cm=itl,

elements=rec.get("elements", []),

note="Resolved by exact Geant4 NIST name from BOM.",

)

bom_toks = _norm_tokens(bom_material)

if not bom_toks:

return ResolvedMaterial(

bom_name=bom_material,

nist_name=None,

score=0.0,

rho_used_g_cm3=rho_part_g_cm3,

radlen_cm=None,

intlen_cm=None,

elements=None,

note="FIXME: Empty/unknown BOM material string; using dummy material.",

)

def _build_custom_from_elements(note_prefix: str) -> Optional[ResolvedMaterial]:

s = raw_bom_material

if not s:

return None

symbols = set(re.findall(r"\b([A-Z][a-z]?)\b", s))

name_to_symbol = {

"aluminum": "Al", "aluminium": "Al", "silicon": "Si", "iron": "Fe", "copper": "Cu",

"beryllium": "Be", "magnesium": "Mg", "manganese": "Mn", "chromium": "Cr", "nickel": "Ni",

"zinc": "Zn", "titanium": "Ti", "tungsten": "W", "lead": "Pb", "tin": "Sn",

}

for t in bom_toks:

if t in name_to_symbol:

symbols.add(name_to_symbol[t])

symbols = [sym for sym in sorted(symbols) if sym in _ELEMENT_TABLE]

if not symbols:

return None

frac = 1.0 / float(len(symbols))

elems: List[dict] = []

for sym in symbols:

Z, A = _ELEMENT_TABLE[sym]

elems.append({"symbol": sym, "Z": Z, "A_g_mol": A, "mass_fraction": frac})

return ResolvedMaterial(

bom_name=bom_material,

nist_name=None,

score=0.0,

rho_used_g_cm3=rho_part_g_cm3,

radlen_cm=None,

intlen_cm=None,

elements=elems,

note=f"FIXME: {note_prefix} No suitable Geant4 NIST material. Emitting placeholder mixture from parsed elements {symbols} with equal mass fractions; please adjust fractions/material.",

)

best = (None, -1.0, 0.0, 0.0) # (nist_name, score, dens_score, token_score)

second = (None, -1.0, 0.0, 0.0)

bom_has_compound = any(t in bom_toks for t in (

"oxide", "dioxide", "carbide", "nitride", "fluoride", "chloride",

"sulfate", "phosphate", "glass", "dioxyde"

))

for nist_name, rec in g4db.items():

nist_toks = _norm_tokens(nist_name)

ts = _token_score(bom_toks, nist_toks)

if ts <= 0.0:

continue

ds = _density_score(rho_part_g_cm3, rec.get("density_g_cm3"))

# Optional hard density filter

if cfg.max_log_density_diff and cfg.max_log_density_diff > 0.0 and rho_part_g_cm3 and rec.get("density_g_cm3"):

try:

if abs(math.log(float(rec.get("density_g_cm3")) / float(rho_part_g_cm3))) > cfg.max_log_density_diff:

continue

except Exception:

pass

nist_has_compound = any(t in nist_toks for t in (

"oxide", "dioxide", "carbide", "nitride", "fluoride", "chloride",

"sulfate", "phosphate", "glass", "dioxyde"

))

compound_pen = cfg.compound_penalty if (nist_has_compound and not bom_has_compound) else 0.0

score = cfg.w_token * ts + cfg.w_density * ds - compound_pen

if score > best[1]:

second = best

best = (nist_name, score, ds, ts)

elif score > second[1]:

second = (nist_name, score, ds, ts)

nist_best, score_best, ds_best, ts_best = best

nist_second, score_second, _, _ = second

if nist_best is None or score_best < cfg.min_score:

custom = _build_custom_from_elements("Could not resolve with enough confidence.")

if custom is not None:

return custom

return ResolvedMaterial(

bom_name=bom_material,

nist_name=None,

score=float(score_best if score_best > 0 else 0.0),

rho_used_g_cm3=rho_part_g_cm3,

radlen_cm=None,

intlen_cm=None,

elements=None,

note="FIXME: Could not resolve BOM material to a Geant4 NIST material with enough confidence; using dummy material.",

)

if score_second > 0 and (score_best - score_second) < cfg.ambiguity_delta:

custom = _build_custom_from_elements(

f"Ambiguous material match (best '{nist_best}' score={score_best:.3f}, second '{nist_second}' score={score_second:.3f})."

)

if custom is not None:

return custom

return ResolvedMaterial(

bom_name=bom_material,

nist_name=None,

score=float(score_best),

rho_used_g_cm3=rho_part_g_cm3,

radlen_cm=None,

intlen_cm=None,

elements=None,

note=f"FIXME: Ambiguous material match (best '{nist_best}' score={score_best:.3f}, second '{nist_second}' score={score_second:.3f}); using dummy material.",

)

rec = g4db[nist_best]

rho_ref = rec.get("density_g_cm3")

# Use NIST density for emission; CAD-derived density is used only for matching.

rho_used = rho_ref

rad = rec.get("radlen_cm")

itl = rec.get("intlen_cm")

return ResolvedMaterial(

bom_name=bom_material,

nist_name=nist_best,

score=float(score_best),

rho_used_g_cm3=rho_used,

radlen_cm=rad,

intlen_cm=itl,

elements=rec.get("elements", []),

note=f"Resolved to '{nist_best}' (token={ts_best:.3f}, density={ds_best:.3f}, score={score_best:.3f}).",

)

def build_volume_to_material_map(

bom_entries: List[BomEntry],

def_names: Dict[str, str],

) -> Dict[str, BomEntry]:

"""

Builds a mapping def_lid -> BomEntry by matching the XCAF display name to:

- exact part_number match

- exact description/name match

- substring match on part_number within the XCAF name

This is heuristic; if nothing matches we keep no assignment for that volume.

"""

# lookup tables

by_part: Dict[str, BomEntry] = {}

by_name: Dict[str, BomEntry] = {}

for e in bom_entries:

if e.part_number_key:

by_part[e.part_number_key] = e

if e.name_key and e.name_key not in by_name:

by_name[e.name_key] = e

out: Dict[str, BomEntry] = {}

for lid, disp in def_names.items():

key = (disp or "").strip()

if not key:

continue

# 1) exact part number

if key in by_part:

out[lid] = by_part[key]

continue

# 2) exact name/description

if key in by_name:

out[lid] = by_name[key]

continue

# 3) substring match on any part number

for pn, e in by_part.items():

if pn and pn in key:

out[lid] = e

break

return out

# -------------------------------

# C++ emission helpers

# -------------------------------

def trsf_to_tgeo(trsf: gp_Trsf, name: str, scale_to_cm: float) -> str:

m = trsf.GetRotation().GetMatrix()

t = trsf.TranslationPart()

return f"""

Double_t {name}_m[9] = {{

{m.Value(1,1)}, {m.Value(1,2)}, {m.Value(1,3)},

{m.Value(2,1)}, {m.Value(2,2)}, {m.Value(2,3)},

{m.Value(3,1)}, {m.Value(3,2)}, {m.Value(3,3)}

}};

TGeoRotation *{name}_rot = new TGeoRotation();

{name}_rot->SetMatrix({name}_m);

TGeoCombiTrans *{name} = new TGeoCombiTrans({t.X()*scale_to_cm}, {t.Y()*scale_to_cm}, {t.Z()*scale_to_cm}, {name}_rot);

"""

def emit_cpp_prelude() -> str:

return """#include <TGeoManager.h>

#include <TFile.h>

#include <fstream>

#include <functional>

#include <stdexcept>

#include <string>

static void LoadFacets(const std::string& file, TGeoTessellated* solid, bool check=false)

{

std::ifstream in(file, std::ios::binary);

if (!in) throw std::runtime_error("Cannot open facet file: " + file);

uint32_t nTri = 0;

in.read(reinterpret_cast<char*>(&nTri), sizeof(nTri));

if (!in) throw std::runtime_error("Bad facet header in: " + file);

for (uint32_t i=0;i<nTri;i++) {

float v[9];

in.read(reinterpret_cast<char*>(v), sizeof(v));

if (!in) throw std::runtime_error("Unexpected EOF in: " + file);

solid->AddFacet(TGeoTessellated::Vertex_t(v[0],v[1],v[2]),

TGeoTessellated::Vertex_t(v[3],v[4],v[5]),

TGeoTessellated::Vertex_t(v[6],v[7],v[8]));

}

solid->CloseShape(check, true);

}

"""

def emit_materials_cpp(

used_materials: Dict[str, ResolvedMaterial],

# key: BOM material string as used in CSV after normalization

) -> Tuple[str, Dict[str, str]]:

"""

Emits C++ code defining TGeoMaterial/TGeoMixture + TGeoMedium for all used materials.

- If a material resolved to a Geant4 NIST entry, emit a physically correct mixture

(element mass fractions) and set RadLen/IntLen (from Geant4) when available.

- If unresolved/ambiguous, emit a dummy material and annotate with FIXME comments.

"""

cpp: List[str] = []

cpp.append(" // Default material/medium (placeholder; can be replaced later)")

cpp.append(" TGeoMaterial *mat_Default = new TGeoMaterial(\"Default\", 0., 0., 0.);")

cpp.append(" TGeoMedium *med_Default = new TGeoMedium(\"Default\", 1, mat_Default);")

cpp.append("")

emitted_el: Dict[str, str] = {}

def _emit_element(el: dict) -> str:

sym = el.get("symbol", "X")

Z = int(el.get("Z", 0))

A = float(el.get("A_g_mol", 0.0))

if sym in emitted_el:

return emitted_el[sym]

safe = sanitize_cpp_name(sym)

var = f"el_{safe}"

cpp.append(f" TGeoElement *{var} = new TGeoElement(\"{sym}\", \"{sym}\", {Z}, {A:.10g});")

emitted_el[sym] = var

return var

medium_var: Dict[str, str] = {"Default": "med_Default"}

next_id = 2

for bom_mat in sorted(used_materials.keys(), key=lambda s: s.lower()):

rm = used_materials[bom_mat]

safe = sanitize_cpp_name(bom_mat)

base = safe

k = 2

while f"med_{safe}" in medium_var.values():

safe = f"{base}_{k}"

k += 1

rho = rm.rho_used_g_cm3 if (rm.rho_used_g_cm3 and rm.rho_used_g_cm3 > 0.0) else 0.0

cpp.append(f" // BOM material: {rm.bom_name}")

cpp.append(f" // {rm.note}")

if rm.elements:

elems = rm.elements

if len(elems) == 1 and abs(float(elems[0].get('mass_fraction', 1.0)) - 1.0) < 1e-6:

el = elems[0]

A = float(el.get("A_g_mol", 0.0))

Z = float(el.get("Z", 0))

cpp.append(f" TGeoMaterial *mat_{safe} = new TGeoMaterial(\"{bom_mat}\", {A:.10g}, {Z:.10g}, {rho:.10g});")

else:

cpp.append(f" TGeoMixture *mat_{safe} = new TGeoMixture(\"{bom_mat}\", {len(elems)}, {rho:.10g});")

for el in elems:

elvar = _emit_element(el)

w = float(el.get("mass_fraction", 0.0))

cpp.append(f" mat_{safe}->AddElement({elvar}, {w:.10g});")

if rm.radlen_cm is not None and rm.intlen_cm is not None:

cpp.append(f" mat_{safe}->SetRadLen({float(rm.radlen_cm):.10g}, {float(rm.intlen_cm):.10g});")

elif rm.radlen_cm is not None:

cpp.append(f" mat_{safe}->SetRadLen({float(rm.radlen_cm):.10g});")

else:

cpp.append(" // FIXME: Unresolved material. Replace with a proper TGeoMaterial/TGeoMixture.")

cpp.append(f" TGeoMaterial *mat_{safe} = new TGeoMaterial(\"{bom_mat}\", 0., 0., {rho:.10g});")

cpp.append(f" TGeoMedium *med_{safe} = new TGeoMedium(\"{bom_mat}\", {next_id}, mat_{safe});")

cpp.append("")

medium_var[bom_mat] = f"med_{safe}"

next_id += 1

return "\n".join(cpp), medium_var

def emit_tessellated_cpp(lid: str, vol_display_name: str, facet_abspath: str, ntriangles: int, medium_var: str) -> str:

safe = sanitize_cpp_name(lid)

shape_name = vol_display_name if vol_display_name else lid

if ntriangles <= 0:

out = []

out.append(f' TGeoBBox *solid_{safe} = new TGeoBBox("{shape_name}", 0.001, 0.001, 0.001);')

out.append(f' TGeoVolume *vol_{safe} = new TGeoVolume("{shape_name}", solid_{safe}, {medium_var});')

return "\n".join(out)

out = []

out.append(f' TGeoTessellated *solid_{safe} = new TGeoTessellated("{shape_name}", {ntriangles});')

out.append(f' LoadFacets("{facet_abspath}", solid_{safe}, check);')

out.append(f' TGeoVolume *vol_{safe} = new TGeoVolume("{shape_name}", solid_{safe}, {medium_var});')

return "\n".join(out)

def emit_assembly_cpp(lid: str, asm_display_name: str) -> str:

safe = sanitize_cpp_name(lid)

name = asm_display_name if asm_display_name else lid

return f' TGeoVolumeAssembly *asm_{safe} = new TGeoVolumeAssembly("{name}");'

# -------------------------------

# Definition graph extraction

# -------------------------------

logical_volumes: Dict[str, list] = {} # def_lid -> triangles

def_names: Dict[str, str] = {} # def_lid -> human display name (may be "")

def_volumes_cm3: Dict[str, float] = {} # def_lid -> volume in cm^3 (leaf only)

assemblies = set() # def_lid

placements = [] # (parent_def_lid, child_def_lid, gp_Trsf local)

top_defs = set() # top definition lids

visited_defs = set() # expanded defs

def cpp_var_for_def(lid: str) -> str:

safe = sanitize_cpp_name(lid)

return f"asm_{safe}" if lid in assemblies else f"vol_{safe}"

def expand_definition(def_label: TDF_Label, shape_tool, meshparam=None, scale_to_cm: float = 1.0):

def_lid = label_id(def_label)

if def_lid in visited_defs:

return

visited_defs.add(def_lid)

nm = label_name(def_label)

if nm and def_lid not in def_names:

def_names[def_lid] = nm

elif def_lid not in def_names:

def_names[def_lid] = ""

children = TDF_LabelSequence()

shape_tool.GetComponents(def_label, children)

has_children = children.Length() > 0

if has_children or shape_tool.IsAssembly(def_label):

assemblies.add(def_lid)

for i in range(children.Length()):

child = children.Value(i + 1)

if shape_tool.IsReference(child):

referred = TDF_Label()

shape_tool.GetReferredShape(child, referred)

child_def_lid = label_id(referred)

loc = shape_tool.GetLocation(child)

trsf = loc.Transformation()

placements.append((def_lid, child_def_lid, trsf))

expand_definition(referred, shape_tool, meshparam=meshparam, scale_to_cm=scale_to_cm)

else:

child_def_lid = label_id(child)

placements.append((def_lid, child_def_lid, gp_Trsf()))

expand_definition(child, shape_tool, meshparam=meshparam, scale_to_cm=scale_to_cm)

return

if shape_tool.IsSimpleShape(def_label):

if def_lid not in logical_volumes:

shape = shape_tool.GetShape(def_label)

# store volume (for density estimation)

try:

def_volumes_cm3[def_lid] = volume_cm3_of_shape(shape, scale_to_cm=scale_to_cm)

except Exception:

def_volumes_cm3[def_lid] = 0.0

do_meshing = (meshparam is not None) and meshparam.get("do_meshing", None) is True

logical_volumes[def_lid] = triangulate_CAD_solid(shape, meshparam=meshparam, scale_to_cm=scale_to_cm) if do_meshing else triangulate_asbbox(shape, scale_to_cm=scale_to_cm)

return

assemblies.add(def_lid)

def extract_graph(step_path: str, meshparam=None, scale_to_cm: float = 1.0):

global logical_volumes, def_names, def_volumes_cm3, assemblies, placements, top_defs, visited_defs

logical_volumes = {}

def_names = {}

def_volumes_cm3 = {}