Simulation Tutorial

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import PyPWA as pwa
from tqdm.notebook import tqdm
import torch as tc
import jpac

/home/salgado/anaconda3/envs/PyPWA/lib/python3.10/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Set up the amplitude for Simulation

The function will be use by the rejection method to "carve" a new distribution into the input simulated data read in next lines. In this example we read: elm : waves defined by (reflectivity(e), angular momentum(l) and angular momentum z projection (m)) resonances: mass and width of resonances input

amp = jpac.Simulate().load_resonances_from_uri(
    'data/elms.csv', 'data/resonance.csv'
)

Read input (flat) generated data (in condensed format)

data = pwa.read("data/helicity.csv")
data

array([(0.000000e+00, 1.90216 , 3.80047 , -1.0749 , 0.4, -0.0263538, 0.850234),
       (1.000000e+00, 0.916137, 1.66023 , -1.43837, 0.4, -0.671785 , 2.61953 ),
       (2.000000e+00, 1.81133 , 5.77896 , -2.34144, 0.4, -0.0775825, 1.26828 ),
       ...,
       (9.855143e+06, 1.32847 , 1.94873 , -2.10288, 0.4, -0.138508 , 0.940464),
       (9.855144e+06, 2.18429 , 2.49181 , -1.8212 , 0.4, -0.314037 , 0.723398),
       (9.855145e+06, 1.57057 , 0.590933, -1.3353 , 0.4, -0.240994 , 1.18741 )],
      dtype={'names': ['EventN', 'theta', 'phi', 'alpha', 'pol', 'tM', 'mass'], 'formats': ['<f8', '<f8', '<f8', '<f8', '<f8', '<f8', '<f8'], 'offsets': [0, 8, 16, 24, 32, 40, 48], 'itemsize': 56, 'aligned': True})

The format of the input data will depend on the Amplitudes (function): In this example the standard HEL angles, polarization angle (alpha) and other information necessary are given for PWA

Produce simulation mask

A boolean file of (False and True)

rejection = pwa.monte_carlo_simulation(amp, data, dict())

Check on the waves and resonances that will be produced

This is a matrix with the weigths of each wave on each resonance

resonance_df = {
    "Resonance": amp.resonance.w0,
    "Res-weight": amp.resonance.cr,
}
for i in range(len(amp.resonance.e)):
    res = amp.resonance
    resonance_df[f"({res.e[i]}, {res.l[i]}, {res.m[i]})"] = res.wave.T[i]

resonance_df = pd.DataFrame(resonance_df)
resonance_df

	Resonance	Res-weight	(1, 0, 0)	(1, 1, 0)	(1, 1, -1)	(1, 1, 1)	(1, 2, 0)	(1, 2, -1)	(1, 2, 1)	(1, 2, -2)	(1, 2, 2)
0	0.980	0.6	1.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
1	1.306	0.4	0.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0
2	1.584	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
3	1.722	0.3	0.0	0.0	0.0	1.0	0.0	0.0	0.0	0.0	0.0

Check on how many events will be kept through the masking

print(f"Removed {len(data) - np.sum(rejection)} events from flat data.")
print(f"Kept {np.sum(rejection)} events.")
print(f"{(np.sum(rejection) / len(data)) * 100}% of events remain.")

Removed 9062608 events from flat data.
Kept 792538 events.
8.041869699342861% of events remain.

Apply mask to input data

new_data will contain the simulated data in the same format that data/datag

new_data = data[rejection]
results = amp.setup(new_data).calculate()

new_data = pd.DataFrame(new_data)
new_data

	EventN	theta	phi	alpha	pol	tM	mass
0	19.0	1.252070	1.568260	-0.902285	0.4	-0.576279	0.968759
1	25.0	1.093960	1.172510	2.855050	0.4	-0.150605	1.771060
2	33.0	1.626260	3.103820	3.026000	0.4	-0.009246	0.846339
3	37.0	2.105310	5.148610	0.079403	0.4	-0.084543	1.458540
4	50.0	2.304940	0.220618	-2.450040	0.4	-0.697125	1.244430
...	...	...	...	...	...	...	...
792533	9855118.0	2.642010	2.581130	-0.056481	0.4	-0.032798	0.944652
792534	9855120.0	2.215190	1.773200	-1.382420	0.4	-0.052298	1.465130
792535	9855130.0	1.103580	4.182710	-3.113100	0.4	-0.098202	1.771670
792536	9855134.0	0.667566	4.076400	-1.236310	0.4	-0.188685	1.679170
792537	9855139.0	0.930979	1.097230	-1.213530	0.4	-0.730447	1.716200

792538 rows × 7 columns

Plot simulated data intensity versus mass

mni = np.empty(len(new_data), dtype=[("mass", float), ("intensity", float)])
mni["mass"] = new_data["mass"]
mni["intensity"] = results.real
mni = pd.DataFrame(mni)
counts, bin_edges = np.histogram(mni["mass"], 200, weights=mni["intensity"])
centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Add y_err to argument list when we have errors
y_err = np.sqrt(counts)
plt.errorbar(centers, counts, y_err, fmt="o")
plt.xlabel(r'$mass$')
plt.ylabel("INtensity")
plt.xlim(.6, 2.)

(0.6, 2.0)

Calculate Phase difference between two waves

In this example first and 3erd waves in amplitude list

phase_diff = np.arctan(
    (amp.vs[0][0] * amp.vs[3][3].conj()).imag
    / (amp.vs[0][0] * amp.vs[3][3].conj()).real
)

Plot PhaseMotion

mnip = np.empty(len(new_data), dtype=[("mass", float), ("phase", float)])
mnip["mass"] = new_data["mass"]
mnip["phase"] = phase_diff
mnip = pd.DataFrame(mnip)
counts, bin_edges = np.histogram(mnip["mass"], 100, weights=mnip["phase"])
centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Add y_err to argument list when we have errors
y_err = np.sqrt(counts)
plt.errorbar(centers, counts, y_err, fmt="o")
plt.xlim(0.6, 2.)

/tmp/ipykernel_203453/414601304.py:9: RuntimeWarning: invalid value encountered in sqrt
  y_err = np.sqrt(counts)

(0.6, 2.0)

Plot phi_HEL vs cosHEL) of simulated data (with 4 different contracts(gamma))

import matplotlib.colors as mcolors

gammas = [0.8, 0.5, 0.3]

fig, axes = plt.subplots(nrows=2, ncols=2)

axes[0, 0].set_title('Linear normalization')
axes[0, 0].hist2d(new_data["phi"], np.cos(new_data["theta"]), bins=100)
plt.xlabel(r'$\phi$')
plt.ylabel("cos($\theta$)")
for ax, gamma in zip(axes.flat[1:], gammas):
    ax.set_title(r'Power law $(\gamma=%1.1f)$' % gamma)
    ax.hist2d(new_data["phi"], np.cos(new_data["theta"]),
        bins=100, norm=mcolors.PowerNorm(gamma))
    
fig.tight_layout()


plt.show()

import matplotlib.colors as mcolors

gammas = [0.8, 0.5, 0.3]

fig, axes = plt.subplots(nrows=2, ncols=2)

axes[0, 0].set_title('Linear normalization')
axes[0, 0].hist2d(new_data["theta"], new_data["phi"], bins=100)

for ax, gamma in zip(axes.flat[1:], gammas):
    ax.set_title(r'Power law $(\gamma=%1.1f)$' % gamma)
    ax.hist2d(new_data["theta"], new_data["phi"],
    bins=100, norm=mcolors.PowerNorm(gamma))
    
fig.tight_layout()

plt.show()