File size: 5,974 Bytes
2b0bffa | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 | """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)
# ── counting / Poisson statistics ─────────────────────────────────
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)
# ── mass smearing ──────────────────────────────────────────────────
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))
# ── significance ───────────────────────────────────────────────────
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))
# ── helpers ─────────────────────────────────────────────────────────
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()
|