Validation of Gain with Monte Carlo

Scott Prahl

July 2025

This notebook compares the calculated gain values against those from a Monte Carlo simulation of light bouncing in spheres.

Since these are stochastic tests, it is likely that a few of the results will be outside the 90% confidence interval.

import numpy as np
import scipy.stats as stats

import iadpython

%config InlineBackend.figure_format='retina'


def p_value(x, x_std_err, y):
    """p-value that x matches known y."""
    # Calculate the z-score
    z_score = (y - x) / x_std_err

    # Calculate the probability from the z-score
    p = 2 * (1 - stats.norm.cdf(abs(z_score)))

    return p


def calculate_confidence_interval(data, confidence=0.90):
    """Compute `calculate_confidence_interval`.

    Args:
        data: Input parameter.
        confidence: Input parameter.

    Returns:
        Computed result from `calculate_confidence_interval`.
    """
    n = len(data)
    mean = np.mean(data)
    std_dev = np.std(data, ddof=1)  # ddof=1 returns the sample standard deviation
    std_err = std_dev / np.sqrt(n)

    # Find the critical value for the 90% confidence interval
    t_critical = stats.t.ppf((1 + confidence) / 2, df=n - 1)

    margin_of_error = t_critical * std_err
    confidence_interval = (mean - margin_of_error, mean + margin_of_error)

    return mean, confidence_interval


def test_gain(s, N=100000):
    """Compute `test_gain`.

    Args:
        s: Input parameter.
        N: Input parameter.
    """
    analytic_gain = s.gain()

    gains, _ = s.do_N_photons_gain_array(N, num_trials=20)
    _, ci = calculate_confidence_interval(gains)

    status = "bad"
    if ci[0] <= analytic_gain <= ci[1]:
        status = "good"

    print("upper 90%% ci       = %.3f" % ci[1])
    print("analytic gain      = %.3f (%s)" % (analytic_gain, status))
    print("lower 90%% ci       = %.3f" % ci[0])


def default_sphere():
    """Compute `default_sphere`.

    Returns:
        Computed result from `default_sphere`.
    """
    R = 30
    d_sphere = 2 * R
    d_sample = 20
    d_third = 15
    d_detector = 10
    s = iadpython.Sphere(
        d_sphere,
        d_sample,
        d_third=d_third,
        d_detector=d_detector,
        r_detector=0.5,
        r_wall=0.75,
        r_std=0.8,
    )
    return s

Check no baffle calculations

case 1, third = empty, sample = empty, no baffle

print("case 1, third = empty, sample = empty, no baffle")
s = default_sphere()
s.sample.uru = 0.0
s.third.uru = 0.00
print(repr(s))
test_gain(s)
print()

case 1, third = empty, sample = empty, no baffle
Sphere: d=60.00, r_wall= 75.0%, r_std= 80.0%, baffle= False
    Sample: d= 20.0mm, uru=  0.0%
     Third: d= 15.0mm, uru=  0.0%
  Detector: d= 10.0mm, uru= 50.0%

upper 90% ci       = 3.699
analytic gain      = 3.518 (good)
lower 90% ci       = 3.515

case 2, third = empty, sample = 50%, no baffle

print("case 2, third = empty, sample = 50%, no baffle")
s = default_sphere()
s.sample.uru = 0.5
s.third.uru = 0.0
print(repr(s))
test_gain(s)
print()

case 2, third = empty, sample = 50%, no baffle
Sphere: d=60.00, r_wall= 75.0%, r_std= 80.0%, baffle= False
    Sample: d= 20.0mm, uru= 50.0%
     Third: d= 15.0mm, uru=  0.0%
  Detector: d= 10.0mm, uru= 50.0%

upper 90% ci       = 3.766
analytic gain      = 3.698 (good)
lower 90% ci       = 3.591

case 3, third = 95%, sample = 50%, no baffle

print("case 3, third = 95%, sample = 50%, no baffle")
s = default_sphere()
s.sample.uru = 0.5
s.third.uru = 0.95
print(repr(s))
test_gain(s)

case 3, third = 95%, sample = 50%, no baffle
Sphere: d=60.00, r_wall= 75.0%, r_std= 80.0%, baffle= False
    Sample: d= 20.0mm, uru= 50.0%
     Third: d= 15.0mm, uru= 95.0%
  Detector: d= 10.0mm, uru= 50.0%

upper 90% ci       = 4.085
analytic gain      = 3.913 (good)
lower 90% ci       = 3.901

Check baffle calculations

case 1, third = empty, sample = empty, baffle

print("case 1, third = empty, sample = empty, baffle")
s = default_sphere()
s.baffle = True
s.sample.uru = 0.0
s.third.uru = 0.00
print(repr(s))
test_gain(s)
print()

case 1, third = empty, sample = empty, baffle
Sphere: d=60.00, r_wall= 75.0%, r_std= 80.0%, baffle= True
    Sample: d= 20.0mm, uru=  0.0%
     Third: d= 15.0mm, uru=  0.0%
  Detector: d= 10.0mm, uru= 50.0%

upper 90% ci       = 3.661
analytic gain      = 3.506 (bad)
lower 90% ci       = 3.525

case 2, third = empty, sample = 50%, baffle

print("case 2, third = empty, sample = 50%, baffle")
s = default_sphere()
s.baffle = True
s.sample.uru = 0.5
s.third.uru = 0.0
print(repr(s))
test_gain(s)
print()

case 2, third = empty, sample = 50%, baffle
Sphere: d=60.00, r_wall= 75.0%, r_std= 80.0%, baffle= True
    Sample: d= 20.0mm, uru= 50.0%
     Third: d= 15.0mm, uru=  0.0%
  Detector: d= 10.0mm, uru= 50.0%

upper 90% ci       = 3.751
analytic gain      = 3.637 (good)
lower 90% ci       = 3.597

case 3, third = 95%, sample = 50%, baffle

print("case 3, third = 95%, sample = 50%, baffle")
s = default_sphere()
s.baffle = True
s.sample.uru = 0.5
s.third.uru = 0.95
print(repr(s))
test_gain(s)

case 3, third = 95%, sample = 50%, baffle
Sphere: d=60.00, r_wall= 75.0%, r_std= 80.0%, baffle= True
    Sample: d= 20.0mm, uru= 50.0%
     Third: d= 15.0mm, uru= 95.0%
  Detector: d= 10.0mm, uru= 50.0%

upper 90% ci       = 3.969
analytic gain      = 3.849 (good)
lower 90% ci       = 3.777

case 4, third = 25, sample = 95%, baffle

print("case 4, third = 25%, sample = 95%, baffle")
s = default_sphere()
s.baffle = True
s.sample.uru = 0.95
s.third.uru = 0.25
print(repr(s))
test_gain(s)

case 4, third = 25%, sample = 95%, baffle
Sphere: d=60.00, r_wall= 75.0%, r_std= 80.0%, baffle= True
    Sample: d= 20.0mm, uru= 95.0%
     Third: d= 15.0mm, uru= 25.0%
  Detector: d= 10.0mm, uru= 50.0%

upper 90% ci       = 3.882
analytic gain      = 3.821 (good)
lower 90% ci       = 3.705