Prediction Tutorial¶

Simulates “true data” corrected by acceptance, according to

the fitted values for the amplitude

[1]:

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

Define waves and amplitude (function) to simulate¶

[2]:

initial=[]
for param_type in ["r", "i"]:
    initial.append(f"{param_type}.1.0.0")
    initial.append(f"{param_type}.1.1.0")
    initial.append(f"{param_type}.1.1.1")
    initial.append(f"{param_type}.1.2.0")
    initial.append(f"{param_type}.1.2.1")
    initial.append(f"{param_type}.1.2.2")
#import AmplitudeOLDfit
#amp = AmplitudeOLDfit.FitAmplitude(initial)
#
import AmplitudeJPACsim
amp = AmplitudeJPACsim.NewAmplitude()
#
#import AmplitudeJPACfit
#amp = AmplitudeJPACfit.FitAmplitude(initial)

Read input (flat) simulated (generated) data¶

[3]:

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

Read flat data in gamp format (full information)

[4]:

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

Define number of bins (MUST be the same as the fitted parameters)

[5]:

nbins=20
bins = pwa.bin_by_range(data, "mass", nbins, .6, 2.0)

Calculate the mass value and number of events in each bin

[6]:

bmass=[]
mcounts=[]
for index, bin in enumerate(bins):
    if len(bin)==0:
        bmass.append(0.)
        mcounts.append(0.)
    else:
        bmass.append(npy.average(bin["mass"]))
        mcounts.append(len(bin))

Read parameters from fit¶

[7]:

par = pwa.read("final_values_JPAC.csv")

Prepare for a binned simulation > Find intensities in each bin and max intensity for each bin.

[8]:

int_values = []
max_values = []
params={}
for index, bin in enumerate(bins):
    for param_type in ["r", "i"]:
        params.update({f"{param_type}.1.0.0":par[f"{param_type}.1.0.0"][index]})
        params.update({f"{param_type}.1.1.0":par[f"{param_type}.1.1.0"][index]})
        params.update({f"{param_type}.1.1.1":par[f"{param_type}.1.1.1"][index]})
        params.update({f"{param_type}.1.2.0":par[f"{param_type}.1.2.0"][index]})
        params.update({f"{param_type}.1.2.1":par[f"{param_type}.1.2.1"][index]})
        params.update({f"{param_type}.1.2.2":par[f"{param_type}.1.2.2"][index]})


    [int,intmax] = pwa.simulate.process_user_function(amp,bin,params,16)
    int_values.append(int)
    max_values.append(intmax)

Simulate events for each bin (produce mask and mask them)¶

[9]:

rejected_bins = []
masked_final_values = []
for int_value, bin in zip(int_values, bins):
    rejection = pwa.simulate.make_rejection_list(int_value, max_values)
    rejected_bins.append(bin[rejection])
    masked_final_values.append(int_value[rejection])

Check how many events simulated by bin

[10]:

for index, simulated_bin in enumerate(rejected_bins):
    print(
        f"Bin {index+1}'s length is {len(simulated_bin)}, "
        f"{(len(simulated_bin) / len(bin)) * 100:.2f}% events were kept"
    )

Bin 1's length is 4850, 1.79% events were kept
Bin 2's length is 6545, 2.42% events were kept
Bin 3's length is 10388, 3.84% events were kept
Bin 4's length is 25492, 9.42% events were kept
Bin 5's length is 65913, 24.34% events were kept
Bin 6's length is 27495, 10.15% events were kept
Bin 7's length is 15419, 5.69% events were kept
Bin 8's length is 17872, 6.60% events were kept
Bin 9's length is 31916, 11.79% events were kept
Bin 10's length is 45793, 16.91% events were kept
Bin 11's length is 20121, 7.43% events were kept
Bin 12's length is 10427, 3.85% events were kept
Bin 13's length is 10301, 3.80% events were kept
Bin 14's length is 12541, 4.63% events were kept
Bin 15's length is 15215, 5.62% events were kept
Bin 16's length is 18517, 6.84% events were kept
Bin 17's length is 17288, 6.39% events were kept
Bin 18's length is 12478, 4.61% events were kept
Bin 19's length is 8226, 3.04% events were kept
Bin 20's length is 5718, 2.11% events were kept

Stak data for all bins in one file (new_data)

[11]:

for index, the_bin in zip(range(len(rejected_bins)), rejected_bins):
    if index ==0:
        new_data=pandas.DataFrame(the_bin)

    new_data =  new_data.append(the_bin,ignore_index=True)

#print(new_data)

[12]:

new_data

[12]:

	EventN	theta	phi	alpha	pol	tM	mass
0	432.0	1.075030	3.11080	-2.279760	0.4	-0.036087	0.730692
1	3277.0	1.106290	1.62386	2.912540	0.4	-0.021695	0.706537
2	5034.0	2.092900	1.77176	-1.243810	0.4	-0.065768	0.730773
3	5547.0	0.201674	-1.60246	-1.418660	0.4	-0.028844	0.701134
4	6527.0	1.325940	3.00931	-1.637510	0.4	-0.023814	0.722924
...	...	...	...	...	...	...	...
387360	9849089.0	1.940630	-1.85137	1.307640	0.4	-0.521166	1.968390
387361	9851359.0	1.055120	-2.61521	-2.743190	0.4	-0.403025	1.960020
387362	9853817.0	0.189324	-1.00286	-1.850720	0.4	-0.370487	1.944500
387363	9854636.0	1.043920	2.97397	-0.604751	0.4	-0.325758	1.971750
387364	9854972.0	0.607401	-1.62378	-1.638410	0.4	-0.282602	1.978580

387365 rows × 7 columns

Calculate the number of expected events versus mass bins¶

[13]:

total_nExp = npy.empty(len(rejected_bins))
for index, the_bin in zip(range(len(rejected_bins)), rejected_bins):

    for param_type in ["r", "i"]:
        params.update({f"{param_type}.1.0.0":par[f"{param_type}.1.0.0"][index]})
        params.update({f"{param_type}.1.1.0":par[f"{param_type}.1.1.0"][index]})
        params.update({f"{param_type}.1.1.1":par[f"{param_type}.1.1.1"][index]})
        params.update({f"{param_type}.1.2.0":par[f"{param_type}.1.2.0"][index]})
        params.update({f"{param_type}.1.2.1":par[f"{param_type}.1.2.1"][index]})
        params.update({f"{param_type}.1.2.2":par[f"{param_type}.1.2.2"][index]})
    amp.setup(the_bin)
    total_nExp[index] = npy.average(amp.calculate(params))

Plot number of true (100% acc) events versus mass

[14]:

mni = npy.empty(len(total_nExp), dtype=[("mass", float), ("int", float)])
mni["mass"] = bmass
mni["int"] = total_nExp
mni = pandas.DataFrame(mni)
counts, bin_edges = npy.histogram(mni["mass"], nbins, weights=mni["int"])

centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Add yerr to argment list when we have errors
yerr = npy.empty(nbins)
yerr = npy.sqrt(counts)
yerr=0.
plt.errorbar(centers,counts, yerr, fmt="s",linestyle="dashed",markersize='10',label="nExp")

#plt.xlim(.6, 3.2)
#plt.ylim(0.,65000)
plt.legend(loc='upper right')
plt.show()

../_images/examples_demo_JPAC_pre_26_0.png

Calculate Phase difference between two waves¶

[15]:

V1=npy.ndarray(len(par))
V2=npy.ndarray(len(par))
V1 = par["r.1.1.1"]+par["i.1.1.1"]*(1j)
V2 = par["r.1.2.1"]+par["i.1.2.1"]*(1j)
if V1.all() != 0 and V2.all() != 0:
    phasediff = npy.arctan(npy.imag(V1*npy.conj(V2))/npy.real(V1*npy.conj(V2)))

Plot PhaseMotion

[16]:

mnip = npy.empty(len(bins), dtype=[("mass", float), ("phase", float)])
mnip["mass"] = bmass
mnip["phase"] = phasediff
mnip = pandas.DataFrame(mnip)
counts, bin_edges = npy.histogram(mnip["mass"], 25, 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(npy.abs(counts))
plt.errorbar(centers,counts, yerr, fmt="o")
plt.xlim(0.6, 3.2)

[16]:

(0.6, 3.2)

../_images/examples_demo_JPAC_pre_30_1.png

Plot phi_HEL vs cos(theta_HEL) of predicted data (with 4 different contracts(gamma))

[17]:

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_pre_32_0.png

Plot cos(theta_HEL) 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"], 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_pre_34_0.png

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

[19]:

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_pre_36_0.png

[ ]:

[20]:

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

../_images/examples_demo_JPAC_pre_38_0.png

Plot mass 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["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_pre_40_0.png

Plot mass versus alpha/Phi

[22]:

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_pre_42_0.png

Write predicted (true) data to disk¶

[23]:

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

Write gamp format predicted to data (true)

[21]:

x = (data["EventN"]).astype(int)
y = (new_data["EventN"]).astype(int)

predm = npy.isin(x,y)
pdatag = datag[predm]
pwa.write("raw_predicted.gamp",pdatag)

Write predicted accepted data > events.pf is a mask of accepted data that has been produced by Geant

[24]:

acc = pwa.read("events.pf")
accn = acc.to_numpy()

mask_acc_phy = npy.logical_and(predm,accn)
pdatag_acc = datag[mask_acc_phy]
pdata_acc = data[mask_acc_phy]

pdata_acc.to_csv("predictedJPAC_ACC.csv", index=False)
pwa.write("acc_predictedJPAC.gamp",pdatag_acc)