Lifetime of Muons Generated by Cosmic Rays — Week 2

2.1 Prelab

After completing the prelab, you will be able to:

1. Perform a weighted exponential fit to synthetic decay data using scipy.optimize.curve_fit.
2. Assign Poisson uncertainties (\(\sqrt{N}\)) to binned count data and explain why this is appropriate.
3. Predict whether the measured muon lifetime in scintillator will be above or below the vacuum value, and explain why.

2.2 In-Lab

After completing the in-lab work, you will be able to:

1. Export MAESTRO data as an SPE file and load it into Python using the provided SPE reader.
2. Apply the timing calibration from Week 1 to convert channel numbers to microseconds.
3. Rebin data to improve statistics and explain the trade-off with time resolution.
4. Perform a weighted exponential fit to extract the muon lifetime, with proper uncertainty estimation.
5. Create a two-panel plot (data + fit and residuals) and evaluate goodness of fit using reduced chi-squared.
6. Compare the measured lifetime to the vacuum value and quantify the discrepancy.
7. Estimate the nuclear capture rate for \(\mu^-\) on carbon.

4.1 Generate Synthetic Data

Use the synthetic data generator to create a test dataset:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Generate synthetic exponential decay data
np.random.seed(42)  # For reproducibility

# Parameters
tau_true = 2.05      # "True" lifetime in microseconds
N0_true = 500        # Amplitude
bg_true = 2.0        # Background counts per bin
n_bins = 100         # Number of time bins
t_max = 10.0         # Maximum time in microseconds

# Create time array and expected counts
time = np.linspace(0.2, t_max, n_bins)  # Start above threshold
expected = N0_true * np.exp(-time / tau_true) + bg_true

# Generate Poisson-distributed counts
counts = np.random.poisson(expected)

4.2 Assign Uncertainties

Each bin contains a count \(N\) drawn from a Poisson distribution. The standard deviation of a Poisson distribution with mean \(\lambda\) is \(\sqrt{\lambda}\). Since our best estimate of \(\lambda\) is the observed count \(N\), we assign:

# Assign Poisson uncertainties
counts_err = np.sqrt(counts.astype(float))

# Handle any bins with zero counts (avoid division by zero)
counts_err = np.where(counts_err > 0, counts_err, 1.0)

4.3 Fit the Exponential Model

The model for muon decay in the detector is:

where \(N_0\) is the amplitude, \(\tau\) is the lifetime, and \(B\) is a constant background from accidental coincidences.

# Define the model
def exponential_decay(t, N0, tau, background):
    return N0 * np.exp(-t / tau) + background

# Initial parameter guesses
p0 = [np.max(counts), 2.0, np.mean(counts[-10:])]

# Weighted fit with Poisson errors
popt, pcov = curve_fit(
    exponential_decay, time, counts,
    p0=p0,
    sigma=counts_err,
    absolute_sigma=True  # Important: tells curve_fit that sigma is in data units
)
perr = np.sqrt(np.diag(pcov))

print(f"N0:         {popt[0]:.1f} +/- {perr[0]:.1f}")
print(f"tau:        {popt[1]:.3f} +/- {perr[1]:.3f} us")
print(f"background: {popt[2]:.2f} +/- {perr[2]:.2f}")

4.4 Plot Data, Fit, and Residuals

Create a two-panel figure following the same pattern as the Gaussian Beams beam profile analysis:

# Generate smooth fit curve
t_smooth = np.linspace(time.min(), time.max(), 500)
fit_smooth = exponential_decay(t_smooth, *popt)

# Residuals at data points
fit_at_data = exponential_decay(time, *popt)
residuals = counts - fit_at_data

# Two-panel plot
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8),
    gridspec_kw={'height_ratios': [3, 1]}, sharex=True)

# Top: data and fit
ax1.errorbar(time, counts, yerr=counts_err, fmt='o', markersize=3,
             capsize=2, label='Synthetic data', alpha=0.8)
ax1.plot(t_smooth, fit_smooth, 'r-', linewidth=2,
         label=f'Fit: τ = {popt[1]:.3f} μs')
ax1.set_ylabel('Counts per bin')
ax1.set_title('Prelab: Synthetic Muon Decay Fit')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Bottom: residuals
ax2.errorbar(time, residuals, yerr=counts_err, fmt='o', markersize=3,
             capsize=2, alpha=0.8)
ax2.axhline(y=0, color='r', linestyle='--')
ax2.set_xlabel('Time (μs)')
ax2.set_ylabel('Residuals')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

4.5 Calculate Reduced Chi-Squared

The chi-squared statistic measures how well the model describes the data:

The reduced chi-squared is \(\chi^2_\nu = \chi^2 / \nu\), where \(\nu = N_\text{data} - N_\text{params}\) is the number of degrees of freedom.

chi2 = np.sum((residuals / counts_err) ** 2)
dof = len(counts) - 3  # 3 fit parameters
chi2_red = chi2 / dof

print(f"Chi-squared:         {chi2:.1f}")
print(f"Degrees of freedom:  {dof}")
print(f"Reduced chi-squared: {chi2_red:.2f}")

Interpretation: \(\chi^2_\nu \approx 1\) indicates a good fit. \(\chi^2_\nu \gg 1\) suggests the model doesn’t describe the data (or uncertainties are underestimated). \(\chi^2_\nu \ll 1\) suggests the uncertainties are overestimated.

4.6 Verify Parameter Recovery

Check that your fit recovers the known input lifetime:

print(f"\nTrue lifetime:     {tau_true:.3f} us")
print(f"Fitted lifetime:   {popt[1]:.3f} +/- {perr[1]:.3f} us")
discrepancy = abs(popt[1] - tau_true) / perr[1]
print(f"Discrepancy:       {discrepancy:.1f} sigma")

The fitted value should agree with the true value within ~1–2\(\sigma\).

4.7 Prelab Prediction

Before coming to lab, consider this question:

The free muon lifetime in vacuum is 2.197 \(\mu\)s. In your experiment, you are measuring the lifetime in plastic scintillator (essentially carbon).

Will the lifetime you measure be:

• (a) Equal to the vacuum value?
• (b) Shorter than the vacuum value?
• (c) Longer than the vacuum value?

Hint: Negative muons (\(\mu^-\)) can undergo nuclear capture (\(\mu^- + p \rightarrow n + \nu_\mu\)) in addition to free decay. Positive muons (\(\mu^+\)) cannot. About 45% of cosmic-ray muons at the surface are negatively charged. Think about how these two populations contribute to the overall observed decay rate.

Record your prediction and reasoning. You will test this quantitatively in the lab.

5.1 Export from MAESTRO

1. In MAESTRO, note the current live time and real time (displayed at the bottom of the window).
2. Stop the acquisition: Acquire → Stop.
3. Export the spectrum: File → Export. Save as an ASCII SPE file (.Spe).
4. Transfer the file to your analysis computer.

5.2 Load with Python

Use the read_spe_file() function from 01_spe_reader.py to load the data:

from importlib import import_module
spe_reader = import_module("01_spe_reader")

# Load your data
data = spe_reader.read_spe_file("your_data_file.Spe")

channels = np.arange(data['num_channels'])
counts = data['counts']

print(f"Channels:     {data['num_channels']}")
print(f"Live time:    {data['live_time']:.0f} s")
print(f"Total counts: {np.sum(counts)}")

5.3 Plot Raw Spectrum

Plot the raw spectrum to see the overall shape:

fig = spe_reader.plot_raw_spectrum(channels, counts,
    title="Muon Lifetime - Raw MCA Spectrum")
plt.show()

What to look for: You should see a peak at low channels (corresponding to short decay times, where the exponential decay rate is highest), falling off toward higher channels. There may be a flat region at high channels from accidental coincidences (background). The first few channels may have zero counts due to the discriminator threshold.

7.1 Question 7: Understanding Uncertainties

(a) For a bin with 100 counts, what is the relative uncertainty (\(\sigma/N\))? What about for a bin with 10 counts? With 1 count?

(b) How does rebinning by a factor of 10 change the counts per bin and the relative uncertainty? Explain why this is beneficial.

(c) Why would it be problematic to include bins with zero counts in the fit? (Think about the Poisson uncertainty formula and what happens in the chi-squared calculation.)

8.1 Question 8: Muon Lifetime Fit

(a) Plot your rebinned data with error bars on both linear and semilog scales:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Linear scale
ax1.errorbar(time, counts_data, yerr=counts_err, fmt='o', markersize=3,
             capsize=2, alpha=0.8)
ax1.set_xlabel('Time (μs)')
ax1.set_ylabel('Counts per bin')
ax1.set_title('Linear Scale')
ax1.grid(True, alpha=0.3)

# Semilog scale
ax2.errorbar(time, counts_data, yerr=counts_err, fmt='o', markersize=3,
             capsize=2, alpha=0.8)
ax2.set_yscale('log')
ax2.set_xlabel('Time (μs)')
ax2.set_ylabel('Counts per bin')
ax2.set_title('Semilog Scale')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

On the semilog plot, a pure exponential decay appears as a straight line. Does your data look approximately linear on the semilog plot? Where does it deviate, and why? (Think about background.)

(b) Fit the exponential decay model \(N(t) = N_0 \, e^{-t/\tau} + B\) to your data. Use the technique from the prelab exercise:

def exponential_decay(t, N0, tau, background):
    return N0 * np.exp(-t / tau) + background

# Estimate initial guesses from the data
p0 = [np.max(counts_data), 2.0, np.mean(counts_data[-10:])]

popt, pcov = curve_fit(exponential_decay, time, counts_data,
                       p0=p0, sigma=counts_err, absolute_sigma=True)
perr = np.sqrt(np.diag(pcov))

print(f"N0:         {popt[0]:.1f} +/- {perr[0]:.1f}")
print(f"tau:        {popt[1]:.3f} +/- {perr[1]:.3f} us")
print(f"background: {popt[2]:.2f} +/- {perr[2]:.2f}")

(c) Create a two-panel plot showing the data with the fit curve (top) and residuals (bottom). Use the code structure from your prelab exercise.

(d) Calculate the reduced chi-squared. Is the fit good? If \(\chi^2_\nu\) is significantly different from 1, discuss possible reasons.

(e) Compare your measured lifetime to the vacuum value of 2.197 \(\mu\)s. How many standard deviations away is your measurement?

9.1 Question 9: Nuclear Capture of Negative Muons

Your measured lifetime in scintillator should be shorter than the vacuum value. This is because negative muons (\(\mu^-\)) can be captured by carbon nuclei in addition to decaying freely:

When a \(\mu^-\) captures, the products (a neutron and a neutrino) typically deposit little or no visible energy in the scintillator, so capture events produce no stop signal. They are simply missing from your data. However, the possibility of capture means that \(\mu^-\) disappear faster than \(\mu^+\), so the decay electrons you do observe from \(\mu^-\) follow a shorter effective lifetime.

(a) Why can \(\mu^-\) undergo nuclear capture but \(\mu^+\) cannot? (Think about charge conservation and what nuclei contain.)

(b) When a \(\mu^-\) can either decay freely or be captured, the rates for the two processes add:

Explain in your own words why the rates add rather than the lifetimes. (Hint: Think about probability per unit time.)

(c) Your data is actually a sum of two exponentials: \(\mu^+\) decays with lifetime \(\tau_\text{free} = 2.197\;\mu\)s, and \(\mu^-\) decays with a shorter effective lifetime \(\tau_{\mu^-}\). The single-exponential fit you performed in Q8 returns a value dominated by the longer-lived \(\mu^+\) component, so it cannot directly give you \(\tau_{\mu^-}\).

To extract the capture rate, fit a two-component model with the \(\mu^+\) lifetime fixed to the known vacuum value:

def two_component_decay(t, N_plus, N_minus, tau_mu_minus, background):
    """Two-exponential model: mu+ (fixed lifetime) + mu- (free lifetime)."""
    tau_free = 2.197  # Fixed to vacuum value
    return (N_plus * np.exp(-t / tau_free)
            + N_minus * np.exp(-t / tau_mu_minus)
            + background)

# Initial guesses
p0 = [0.55 * np.max(counts_data), 0.45 * np.max(counts_data), 1.0,
      np.mean(counts_data[-10:])]

# Bounds: all amplitudes positive, tau_mu_minus between 0.1 and 2.197 us
bounds = ([0, 0, 0.1, 0], [np.inf, np.inf, 2.197, np.inf])

popt2, pcov2 = curve_fit(two_component_decay, time, counts_data,
                         p0=p0, sigma=counts_err, absolute_sigma=True,
                         bounds=bounds)
perr2 = np.sqrt(np.diag(pcov2))

print(f"N+:          {popt2[0]:.1f} +/- {perr2[0]:.1f}")
print(f"N-:          {popt2[1]:.1f} +/- {perr2[1]:.1f}")
print(f"tau_mu-:     {popt2[2]:.3f} +/- {perr2[2]:.3f} us")
print(f"background:  {popt2[3]:.2f} +/- {perr2[3]:.2f}")

(d) From your fitted \(\tau_{\mu^-}\), extract the capture lifetime using Equation 6:

\[\frac{1}{\tau_\text{capture}} = \frac{1}{\tau_{\mu^-}} - \frac{1}{\tau_\text{free}}\]

tau_mu_minus = popt2[2]
tau_free = 2.197
tau_capture = 1.0 / (1.0/tau_mu_minus - 1.0/tau_free)
print(f"Estimated capture lifetime: {tau_capture:.1f} us")
print(f"Accepted value for carbon:  ~2.0 us")

The accepted value for the \(\mu^-\) capture lifetime in carbon is approximately 2.0 \(\mu\)s. How does your estimate compare? Given the limited statistics (~5000 events split between two components), do not be surprised if your uncertainty is large.

(e) Compare the reduced chi-squared of the two-component fit to your single-exponential fit from Q8. Does the two-component model fit the data better? With only ~5000 events, the improvement may be small — discuss why resolving two exponentials with similar lifetimes is inherently difficult.