Simulation Tutorial¶

[1]:

import PyPWA as pwa
import numpy as npy
import pandas
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

Define (import) amplitude (function) to simulate¶

The function will be use by the rejection method to “carve” a new distribution into the input simulated data read in next lines.

[2]:

#
import AmplitudeJPACsim
amp = AmplitudeJPACsim.NewAmplitude()

Read input (flat) simulated data (in condense format)¶

[3]:

data = pwa.read("etapiHEL2_flat.txt")

Read data full (from gamp files)

[4]:

datag = pwa.read("../TUTORIAL_FILES/raw_events.gamp")

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

[5]:

data

[5]:

	EventN	theta	phi	alpha	pol	tM	mass
0	0.0	1.902160	3.800470	-1.074900	0.4	-0.026354	0.850234
1	1.0	0.916137	1.660230	-1.438370	0.4	-0.671785	2.619530
2	2.0	1.811330	5.778960	-2.341440	0.4	-0.077582	1.268280
3	3.0	2.000940	5.250000	2.709120	0.4	-0.726730	1.253960
4	4.0	1.871130	2.697150	0.239302	0.4	-0.156830	0.927206
...	...	...	...	...	...	...	...
9855141	9855141.0	1.354170	5.931830	0.784408	0.4	-0.054668	1.083350
9855142	9855142.0	2.335000	3.347330	-0.748025	0.4	-0.844509	1.483240
9855143	9855143.0	1.328470	1.948730	-2.102880	0.4	-0.138508	0.940464
9855144	9855144.0	2.184290	2.491810	-1.821200	0.4	-0.314037	0.723398
9855145	9855145.0	1.570570	0.590933	-1.335300	0.4	-0.240994	1.187410

9855146 rows × 7 columns

Produce simulation mask¶

A boolean file of (False and True)

[6]:

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

Check on the waves and resonances that will be produced > This is a matrix with the weigths of each wave on each resonance

[7]:

amp.setup(data)
table =[]
tabler=[]

for r in range(amp.resonance.resonance_index):
    tabler.append(amp.resonance.W0[r])
    tabler.append(amp.resonance.Cr[r])
    for w in range(amp.resonance.wave_index):

        tabler.append(amp.resonance.Wave[r][w])
    table.append(tabler)
    tabler=[]

from tabulate import tabulate
headers=["Resonance","Res-weight"]
for w in range(amp.resonance.wave_index):
    headers.append(amp.resonance.wave_data[w])

print(tabulate(table,headers))

  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.98           0.65            1          0               0          0           0                0         0                0         0
      1.306          0.35            0          0               0          0           0.33             0         0.33             0         0.33
      1.584          0.08            0          0.5             0          0.5         0                0         0                0         0
      1.722          0.1             0          0               0          0           0.33             0         0.33             0         0.33

Check on how many events will be kept through the masking

[8]:

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

Removed 9454013 events from flat data.
Kept 401133 events.
4.070289775514234% of events remain.

Apply mask to input data¶

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

[9]:

new_data = data[rejection]

[10]:

new_data

[10]:

	EventN	theta	phi	alpha	pol	tM	mass
137	137.0	2.206940	1.953870	-0.247381	0.4	-0.120428	0.995635
171	171.0	1.775060	1.264290	2.600410	0.4	-0.345176	1.167540
200	200.0	0.128817	2.396510	-0.063606	0.4	-0.229984	2.269060
205	205.0	0.696665	1.981680	3.079090	0.4	-0.184068	1.885680
237	237.0	1.051710	1.753500	-2.496390	0.4	-1.500240	1.334010
...	...	...	...	...	...	...	...
9855028	9855028.0	0.858080	2.779910	-0.051326	0.4	-0.086823	0.981714
9855041	9855041.0	1.620860	0.822745	1.592530	0.4	-0.451578	1.311310
9855064	9855064.0	0.734074	5.887190	-1.768290	0.4	-0.099859	1.098970
9855075	9855075.0	1.502800	1.512360	0.860545	0.4	-0.162814	1.710660
9855113	9855113.0	1.324710	4.971240	-1.140430	0.4	-0.051079	1.154320

401133 rows × 7 columns

Mask full data format

[11]:

new_datag = datag[rejection]

Write new_data into a Pandas Dateframe

[12]:

amp.setup(new_data)
new_data = pandas.DataFrame(new_data)

Plot simulated data intensity versus mass¶

[13]:

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

# Add yerr to argment list when we have errors
yerr = npy.empty(100)
yerr = npy.sqrt(counts)
plt.errorbar(centers,counts, yerr, fmt="o")
#plt.yscale("log")
plt.xlim(.6, 2.)

[13]:

(0.6, 2.0)

../_images/examples_demo_JPAC_sim_24_1.png

Calculate Phase difference between two waves¶

In this example first and 3erd waves in amplitude list

[14]:

#if amp.Vs[1]15].all() != 0 and amp.Vs[0][0].all() != 0:
phasediff = npy.arctan(npy.imag(amp.Vs[2][3]*amp.Vs[1][6].conjugate())/npy.real(amp.Vs[2][3]*amp.Vs[1][6].conjugate()))
#phasediff = npy.arctan(npy.imag(amp.Vs[2][1]*amp.Vs[1][2].conjugate())/npy.real(amp.Vs[2][1]*amp.Vs[1][2].conjugate()))

Plot PhaseMotion

[15]:

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

# Add yerr to argment list when we have errors
yerr = npy.empty(100)
yerr = npy.sqrt(counts)
plt.errorbar(centers,counts, yerr, fmt="o")
plt.xlim(0.6, 2.)

[15]:

(0.6, 2.0)

../_images/examples_demo_JPAC_sim_28_1.png

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

[16]:

import matplotlib.colors as mcolors
from numpy.random import multivariate_normal

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"], npy.cos(new_data["theta"]), bins=100)
#axes[0, 0].hist2d(cut_list["phi"], npy.cos(cut_list["theta"]), 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["phi"], npy.cos(new_data["theta"]),
        bins=100, norm=mcolors.PowerNorm(gamma))

fig.tight_layout()

plt.show()

../_images/examples_demo_JPAC_sim_30_0.png

Plot cos(theta_HEL) vs mass for simulated data (with 4 different contrasts)

[17]:

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["mass"], npy.cos(new_data["theta"]),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["mass"], npy.cos(new_data["theta"]),
              bins=100, norm=mcolors.PowerNorm(gamma))


fig.tight_layout()
plt.xlim(.6, 2.)
plt.show()

../_images/examples_demo_JPAC_sim_32_0.png

Plot phiHEL vs mass for simulated data (with 4 different contrasts)

[18]:

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["mass"], 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["mass"], new_data["phi"],
              bins=100, norm=mcolors.PowerNorm(gamma))


fig.tight_layout()
plt.xlim(.6, 2.)
plt.show()

../_images/examples_demo_JPAC_sim_34_0.png

Histogram of alpha/Phi > alpha/Phi is the polarization angle

[19]:

plt.hist(new_data["alpha"],50)
plt.show()

../_images/examples_demo_JPAC_sim_36_0.png

Plot mass versus alpha/Phi

[20]:

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["mass"], new_data["alpha"],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["mass"], new_data["alpha"],
              bins=100, norm=mcolors.PowerNorm(gamma))


fig.tight_layout()
plt.xlim(.6, 2.)
plt.show()

../_images/examples_demo_JPAC_sim_38_0.png

Plot phi_HEL versus alpha/Phi

[21]:

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"],new_data["alpha"],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["phi"], new_data["alpha"],
              bins=100, norm=mcolors.PowerNorm(gamma))


fig.tight_layout()

plt.show()

../_images/examples_demo_JPAC_sim_40_0.png

Write simulated data to disk¶

[23]:

new_data.to_csv("simdata_JPAC.csv", index=False)

Write gamp data after maked

[24]:

pwa.write("raw_simulated_JPAC.gamp",new_datag)

Calculate Moments (for JPAC or Std) and Asymmetries¶

[25]:

H000,H010,H011,H020,H021,H022,H100,H110,H111,H120,H121,H122,sigma4,sigmay = amp.calculate_moments_JPAC()
#H00,H11,H10,H20,H21,H22 = amp.calculate_moments_STD()

Plot (all) Moments versus mass

[26]:

plt.scatter(new_data["mass"],H000,LABEL="H000")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H010,LABEL="H010")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H011,LABEL="H011")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H020,LABEL="H020")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H021,LABEL="H021")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H022,LABEL="H022")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H100,LABEL="H100")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H110,LABEL="H110")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H111,LABEL="H111")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H120,LABEL="H120")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H121,LABEL="H121")
plt.legend(loc='upper right')
plt.scatter(new_data["mass"],H122,LABEL="H122")
plt.legend(loc='upper right')
plt.xlim(0.6, 2.)

[26]:

(0.6, 2.0)

../_images/examples_demo_JPAC_sim_48_1.png

PLot each moment vs mass

[27]:

plt.xlim(0.6, 2.)
plt.scatter(new_data["mass"],H000,LABEL="H000")
plt.scatter(new_data["mass"],H100,LABEL="H100")
plt.legend(loc='upper right')
plt.show()
plt.xlim(0.6, 2.)
plt.scatter(new_data["mass"],H010,LABEL="H010")
plt.scatter(new_data["mass"],H110,LABEL="H110")
plt.legend(loc='upper right')
plt.show()
plt.xlim(0.6, 2.)
plt.scatter(new_data["mass"],H011,LABEL="H011")
plt.scatter(new_data["mass"],H111,LABEL="H111")
plt.legend(loc='upper right')
plt.show()
plt.xlim(0.6, 2.)
plt.scatter(new_data["mass"],H020,LABEL="H020")
plt.scatter(new_data["mass"],H120,LABEL="H120")
plt.legend(loc='upper right')
plt.show()
plt.xlim(0.6, 2.)
plt.scatter(new_data["mass"],H021,LABEL="H021")
plt.scatter(new_data["mass"],H121,LABEL="H121")
plt.legend(loc='upper right')
plt.show()
plt.xlim(0.6, 2.)
plt.scatter(new_data["mass"],H022,LABEL="H022")
plt.scatter(new_data["mass"],H122,LABEL="H122")
plt.legend(loc='upper right')
plt.show()

../_images/examples_demo_JPAC_sim_50_0.png

../_images/examples_demo_JPAC_sim_50_1.png

../_images/examples_demo_JPAC_sim_50_2.png

../_images/examples_demo_JPAC_sim_50_3.png

../_images/examples_demo_JPAC_sim_50_4.png

../_images/examples_demo_JPAC_sim_50_5.png

Plot asymmetry Sigma_4PI

[28]:

plt.xlim(0.6, 2.)
plt.scatter(new_data["mass"],sigma4,LABEL="sigma4pi")
plt.legend(loc='upper right')
plt.show()

../_images/examples_demo_JPAC_sim_52_0.png

Plot asymmetry Sigma_y

[29]:

plt.xlim(0.6, 2.)
plt.scatter(new_data["mass"],sigmay,LABEL="sigmay")
plt.legend(loc='upper right')

[29]:

<matplotlib.legend.Legend at 0x7f0c54a93390>

../_images/examples_demo_JPAC_sim_54_1.png

[ ]: