PDF.py

import lhapdf
import numpy as np
import matplotlib.pyplot as plt

# --- PDF Set Configuration ---
PDF_SET_NAME = "CT18NLO"
# Note: LHAPDF uses Q, not Q^2. So we use Q = 10 GeV.
Q_FIXED = 10.0  

# Load the PDF set. The "0" means we are using the central member of the set.
try:
    pdf = lhapdf.mkPDF(PDF_SET_NAME, 0)
except Exception as e:
    print(f"Error: Could not load PDF set '{PDF_SET_NAME}'.")
    print(f"Please make sure you have installed it by running: lhapdf install {PDF_SET_NAME}")
    exit()

# --- Setup for Plotting ---
# We use a log scale for x, as that's where the interesting physics is.
x_values = np.logspace(-4, -0.01, 200) # x from 0.0001 to ~0.977

# Quark charges (e_i^2)
# We only need up, down, strange, charm for a good approximation
e_u_sq = (2/3)**2
e_d_sq = (-1/3)**2
e_s_sq = (-1/3)**2
e_c_sq = (2/3)**2

# --- Simulation A: Plot Individual PDFs ---

print(f"Calculating individual PDFs at Q = {Q_FIXED} GeV...")

# Get PDF values for each flavor using their "PDG ID"
# 2 = up, 1 = down, 3 = strange, 4 = charm, 21 = gluon
# Negative IDs are for anti-quarks
u_pdf = np.array([pdf.xfxQ(2, x, Q_FIXED) for x in x_values])
d_pdf = np.array([pdf.xfxQ(1, x, Q_FIXED) for x in x_values])
g_pdf = np.array([pdf.xfxQ(21, x, Q_FIXED) for x in x_values])
u_bar_pdf = np.array([pdf.xfxQ(-2, x, Q_FIXED) for x in x_values])
d_bar_pdf = np.array([pdf.xfxQ(-1, x, Q_FIXED) for x in x_values])

# Calculate valence quarks: u_v = u - u_bar
u_valence = u_pdf - u_bar_pdf
d_valence = d_pdf - d_bar_pdf

plt.figure(1, figsize=(10, 7))
plt.plot(x_values, u_valence, label=r'$x \cdot u_v(x)$ (up valence)', linewidth=2)
plt.plot(x_values, d_valence, label=r'$x \cdot d_v(x)$ (down valence)', linewidth=2)
plt.plot(x_values, g_pdf, label=r'$x \cdot g(x)$ (gluon)', linestyle='--', linewidth=2)
plt.plot(x_values, 2 * u_bar_pdf, label=r'$2x \cdot \bar{u}(x)$ (up sea)', linestyle=':')
plt.plot(x_values, 2 * d_bar_pdf, label=r'$2x \cdot \bar{d}(x)$ (down sea)', linestyle=':')

plt.title(f'Parton Distribution Functions at $Q = {Q_FIXED}$ GeV (Set: {PDF_SET_NAME})')
plt.xlabel('$x$ (Momentum Fraction)')
plt.ylabel('$x \cdot f(x, Q^2)$')
plt.xscale('log') # Use log scale for x-axis
plt.ylim(0, 0.8)  # Set y-axis limits
plt.xlim(1e-4, 1.0) # Set x-axis limits
plt.legend()
plt.grid(True, which="both", ls="--")
print("Plot 1 (Individual PDFs) created.")


# --- Simulation B: Plot F2(x, Q^2) and Show Scaling Violations ---

def calculate_F2(x_val, Q_val):
    """
    Calculates the F2 structure function for a proton at a given x and Q.
    F2 = x * sum[ e_i^2 * (f_i(x,Q) + f_bar_i(x,Q)) ]
    """
    
    # Get PDF values (f_i)
    # We use Q, not Q^2, for lhapdf.xfxQ
    u = pdf.xfxQ(2, x_val, Q_val)
    d = pdf.xfxQ(1, x_val, Q_val)
    s = pdf.xfxQ(3, x_val, Q_val)
    c = pdf.xfxQ(4, x_val, Q_val)
    
    # Get anti-PDF values (f_bar_i)
    u_bar = pdf.xfxQ(-2, x_val, Q_val)
    d_bar = pdf.xfxQ(-1, x_val, Q_val)
    s_bar = pdf.xfxQ(-3, x_val, Q_val)
    c_bar = pdf.xfxQ(-4, x_val, Q_val)
    
    # Apply the formula.
    # Note: pdf.xfxQ already returns x*f(x,Q), so we don't multiply by x again.
    term_u = e_u_sq * (u + u_bar)
    term_d = e_d_sq * (d + d_bar)
    term_s = e_s_sq * (s + s_bar)
    term_c = e_c_sq * (c + c_bar)
    
    return term_u + term_d + term_s + term_c


# Define the Q^2 values (in GeV^2) you want to test
Q2_values = [10, 100, 10000] 

print("\nCalculating F2 for scaling violation plot...")

plt.figure(2, figsize=(10, 7))

for Q2 in Q2_values:
    # IMPORTANT: LHAPDF functions take Q, not Q^2
    Q = np.sqrt(Q2)
    
    # Calculate F2 for all our x_values at this Q
    F2_results = [calculate_F2(x, Q) for x in x_values]
    
    # Plot the result
    plt.plot(x_values, F2_results, label=f'$Q^2 = {Q2}$ GeV$^2$', linewidth=2)

plt.title(f'Scaling Violations in Proton $F_2(x, Q^2)$ (Set: {PDF_SET_NAME})')
plt.xlabel('$x$ (Momentum Fraction)')
plt.ylabel('$F_2(x, Q^2)$')
plt.xscale('log') # Use log scale for x-axis
plt.xlim(1e-4, 1.0) # Set x-axis limits
plt.legend()
plt.grid(True, which="both", ls="--")
print("Plot 2 (F2 Scaling) created.")

# Show both plots
print("\nDisplaying plots...")
plt.show()