| """Stochastic noise model for the LHC simulator.
|
|
|
| All randomness is funneled through a single seeded ``numpy.Generator`` so
|
| episodes are reproducible. The methods are physics-flavoured: Poisson event
|
| counts, Gaussian-smeared masses, log-normal cross-sections, false discovery
|
| helpers, and quality degradation.
|
| """
|
|
|
| from __future__ import annotations
|
|
|
| from typing import List
|
|
|
| import numpy as np
|
|
|
|
|
| class NoiseModel:
|
| """Centralised noise generator for the CERN simulator."""
|
|
|
| def __init__(self, seed: int = 42):
|
| self.rng = np.random.default_rng(seed)
|
|
|
| def reseed(self, seed: int) -> None:
|
| self.rng = np.random.default_rng(seed)
|
|
|
|
|
|
|
| def poisson(self, lam: float) -> int:
|
| return int(self.rng.poisson(max(lam, 0.0)))
|
|
|
| def signal_yield(
|
| self,
|
| cross_section_fb: float,
|
| luminosity_fb: float,
|
| branching: float,
|
| efficiency: float,
|
| trigger_efficiency: float,
|
| ) -> int:
|
| """Expected signal events ~ σ × L × BR × ε_reco × ε_trig + Poisson noise.
|
|
|
| BR = branching ratio of the decay channel.
|
| ε_reco = channel reconstruction efficiency.
|
| ε_trig = trigger acceptance.
|
| """
|
| mu = cross_section_fb * luminosity_fb * branching * efficiency * trigger_efficiency
|
| return self.poisson(mu)
|
|
|
| def background_yield(
|
| self,
|
| baseline_per_fb: float,
|
| luminosity_fb: float,
|
| qcd_strength: float,
|
| trigger_efficiency: float,
|
| ) -> int:
|
| """Expected background events scale linearly with luminosity."""
|
| mu = baseline_per_fb * luminosity_fb * qcd_strength * trigger_efficiency
|
| return self.poisson(mu)
|
|
|
|
|
|
|
| def smear_mass(
|
| self,
|
| true_mass_gev: float,
|
| resolution_gev: float,
|
| scale_offset_gev: float = 0.0,
|
| ) -> float:
|
| return float(self.rng.normal(true_mass_gev + scale_offset_gev, resolution_gev))
|
|
|
| def fit_mass_estimate(
|
| self,
|
| true_mass_gev: float,
|
| n_signal: int,
|
| resolution_gev: float,
|
| scale_offset_gev: float,
|
| ) -> float:
|
| """Fitted mass ≈ true mass + Gaussian error scaling like 1/√N_signal."""
|
| n_eff = max(n_signal, 1)
|
| sigma = resolution_gev / np.sqrt(n_eff)
|
| return float(self.rng.normal(true_mass_gev + scale_offset_gev, sigma))
|
|
|
| def fit_mass_uncertainty(
|
| self,
|
| n_signal: int,
|
| resolution_gev: float,
|
| ) -> float:
|
| """Statistical mass uncertainty from a peak with N_signal events."""
|
| n_eff = max(n_signal, 1)
|
| return float(resolution_gev / np.sqrt(n_eff))
|
|
|
|
|
|
|
| def asimov_significance(
|
| self,
|
| n_signal: int,
|
| n_background: int,
|
| nuisance_inflation: float = 0.0,
|
| ) -> float:
|
| """Asymptotic Asimov-style significance Z = √(2[(s+b) ln(1+s/b) - s]).
|
|
|
| A small nuisance_inflation term in [0,1] shrinks Z to mimic systematic
|
| penalties when calibration / systematics studies are skipped.
|
| """
|
| if n_background <= 0:
|
| return 0.0
|
| s = float(n_signal)
|
| b = float(n_background)
|
| if s <= 0:
|
| return 0.0
|
| term = (s + b) * np.log(1.0 + s / b) - s
|
| z = float(np.sqrt(max(2.0 * term, 0.0)))
|
| return float(z * (1.0 - nuisance_inflation))
|
|
|
|
|
|
|
| def coin_flip(self, p: float) -> bool:
|
| return bool(self.rng.random() < p)
|
|
|
| def jitter(self, mean: float, sigma: float) -> float:
|
| return float(self.rng.normal(mean, sigma))
|
|
|
| def quality_degradation(self, base_quality: float, factors: List[float]) -> float:
|
| q = base_quality
|
| for f in factors:
|
| q *= f
|
| return float(np.clip(q + self.rng.normal(0, 0.02), 0.0, 1.0))
|
|
|
| def sample_qc_metric(
|
| self, mean: float, std: float, clip_lo: float = 0.0, clip_hi: float = 1.0
|
| ) -> float:
|
| return float(np.clip(self.rng.normal(mean, std), clip_lo, clip_hi))
|
|
|
| def histogram(
|
| self,
|
| n_signal: int,
|
| n_background: int,
|
| true_mass_gev: float,
|
| resolution_gev: float,
|
| window_lo_gev: float,
|
| window_hi_gev: float,
|
| n_bins: int = 40,
|
| background_alpha: float = -2.5,
|
| ) -> List[int]:
|
| """Generate a coarse invariant-mass histogram.
|
|
|
| Signal is Gaussian around the (smeared) true mass with width
|
| =resolution; background is a falling power-law shape.
|
| """
|
| if window_hi_gev <= window_lo_gev:
|
| return [0] * n_bins
|
| edges = np.linspace(window_lo_gev, window_hi_gev, n_bins + 1)
|
| centers = 0.5 * (edges[:-1] + edges[1:])
|
|
|
| sig_mu = true_mass_gev
|
| sig_pdf = np.exp(-0.5 * ((centers - sig_mu) / max(resolution_gev, 1e-3)) ** 2)
|
| sig_pdf /= max(sig_pdf.sum(), 1e-9)
|
|
|
| bg_pdf = np.power(np.clip(centers, 1.0, None), background_alpha)
|
| bg_pdf /= max(bg_pdf.sum(), 1e-9)
|
|
|
| sig_counts = self.rng.multinomial(max(n_signal, 0), sig_pdf)
|
| bg_counts = self.rng.multinomial(max(n_background, 0), bg_pdf)
|
| return (sig_counts + bg_counts).astype(int).tolist()
|
|
|
