Examples for config.yml file

Configuration file config.yml use YAML (https://yaml.org) format to describe decay process.

The main parts of config.yml is decay and particle.

The decay part describe the particle (or an id of a list of particle) decay into which particles, it can be a list of a list of list. A list means that there is ony one decay mode, A list of list is the list of possible decay mode. The list item can be the particle name (or a dict to describe the decay parameters). All name should appear in particle part.

The particle part describe the parameters of particles. There are two special parts $top and $finals describe the top and finals particles. The other parts are lists of particle name or dicts of particle parameters. The list is the same type particle in decay part. The dict is the parameters of the particle name.

config_str = """

decay:
    A:
       - [R1, B]
       - [R2, C]
       - [R3, D]
    R1: [C, D]
    R2: [B, D]
    R3: [B, C]

particle:
    $top:
       A: { mass: 1.86, J: 0, P: -1}
    $finals:
       B: { mass: 0.494, J: 0, P: -1}
       C: { mass: 0.139, J: 0, P: -1}
       D: { mass: 0.139, J: 0, P: -1}
    R1: [ R1_a, R1_b ]
    R1_a: { mass: 0.7, width: 0.05, J: 1, P: -1}
    R1_b: { mass: 0.5, width: 0.05, J: 0, P: +1}
    R2: { mass: 0.824, width: 0.05, J: 0, P: +1}
    R3: { mass: 0.824, width: 0.05, J: 0, P: +1}

"""

The config file can be loaded by yaml library.

import matplotlib.pyplot as plt
import yaml

from tf_pwa.config_loader import ConfigLoader
from tf_pwa.histogram import Hist1D

The simple way to create config is write it to config.yml file and then you can load it as config = ConfigLoader("config.yml"). Here we used config_str directly.

config = ConfigLoader(yaml.full_load(config_str))

We set parameters to a blance value. And we can generate some toy data and calclute the weights The full params can be found by print(config.get_params()).

input_params = {
    "A->R1_a.BR1_a->C.D_total_0r": 6.0,
    "A->R1_b.BR1_b->C.D_total_0r": 1.0,
    "A->R2.CR2->B.D_total_0r": 2.0,
    "A->R3.DR3->B.C_total_0r": 1.0,
}
config.set_params(input_params)

True

Here we generate some toy data and phsp mc to show the model

data = config.generate_toy(1000)
phsp = config.generate_phsp(10000)

8.5%[▓▓▓▓>---------------------------------------------] 0.44/5.14s eff: 90.000000%
59.3%[▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓>--------------------] 1.36/2.30s eff: 7.031889%
100.3%[▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓>] 2.19/2.18s eff: 3.823378%
100.0%[▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓] 2.19/2.19s  eff: 3.685080%

You can also fit the data fit to the data We can omit the args when writen in config.yml

fit_result = config.fit([data], [phsp])
# After the fit, you can get the uncertaities as
err = config.get_params_error(fit_result, [data], [phsp])

Using Model
decay chains included:
   [A->R1_a+B, R1_a->C+D]  ls:  ((1, 1),) ((1, 0),)
   [A->R1_b+B, R1_b->C+D]  ls:  ((0, 0),) ((0, 0),)
   [A->R2+C, R2->B+D]  ls:  ((0, 0),) ((0, 0),)
   [A->R3+D, R3->B+C]  ls:  ((0, 0),) ((0, 0),)

########### initial parameters
{
  "R1_a_mass": 0.7,
  "R1_a_width": 0.05,
  "R1_b_mass": 0.5,
  "R1_b_width": 0.05,
  "R2_mass": 0.824,
  "R2_width": 0.05,
  "R3_mass": 0.824,
  "R3_width": 0.05,
  "A->R1_a.BR1_a->C.D_total_0r": 6.0,
  "A->R1_a.BR1_a->C.D_total_0i": 0.0,
  "A->R1_a.B_g_ls_0r": 1.0,
  "A->R1_a.B_g_ls_0i": 0.0,
  "R1_a->C.D_g_ls_0r": 1.0,
  "R1_a->C.D_g_ls_0i": 0.0,
  "A->R1_b.BR1_b->C.D_total_0r": 1.0,
  "A->R1_b.BR1_b->C.D_total_0i": 1.6125033928659116,
  "A->R1_b.B_g_ls_0r": 1.0,
  "A->R1_b.B_g_ls_0i": 0.0,
  "R1_b->C.D_g_ls_0r": 1.0,
  "R1_b->C.D_g_ls_0i": 0.0,
  "A->R2.CR2->B.D_total_0r": 2.0,
  "A->R2.CR2->B.D_total_0i": 1.6836207330993584,
  "A->R2.C_g_ls_0r": 1.0,
  "A->R2.C_g_ls_0i": 0.0,
  "R2->B.D_g_ls_0r": 1.0,
  "R2->B.D_g_ls_0i": 0.0,
  "A->R3.DR3->B.C_total_0r": 1.0,
  "A->R3.DR3->B.C_total_0i": -2.2865658090068535,
  "A->R3.D_g_ls_0r": 1.0,
  "A->R3.D_g_ls_0i": 0.0,
  "R3->B.C_g_ls_0r": 1.0,
  "R3->B.C_g_ls_0i": 0.0
}
initial NLL:  tf.Tensor(-1009.0690556201171, shape=(), dtype=float64)
nll_grad  cost time: 0.24228906631469727
nll_grad  cost time: 0.23493552207946777
nll_grad  cost time: 0.23424363136291504
tf.Tensor(-1009.9362024564653, shape=(), dtype=float64)
nll_grad  cost time: 0.234511137008667
tf.Tensor(-1010.1607292561866, shape=(), dtype=float64)
nll_grad  cost time: 0.2311232089996338
nll_grad  cost time: 0.23242688179016113
tf.Tensor(-1010.2383755695182, shape=(), dtype=float64)
nll_grad  cost time: 0.2350614070892334
nll_grad  cost time: 0.23209166526794434
tf.Tensor(-1010.262177756922, shape=(), dtype=float64)
nll_grad  cost time: 0.23210716247558594
tf.Tensor(-1010.2676491830634, shape=(), dtype=float64)
nll_grad  cost time: 0.23252654075622559
nll_grad  cost time: 0.2337489128112793
tf.Tensor(-1010.3200291979856, shape=(), dtype=float64)
nll_grad  cost time: 0.23587965965270996
tf.Tensor(-1010.3397647468155, shape=(), dtype=float64)
nll_grad  cost time: 0.23921632766723633
tf.Tensor(-1010.3775503993093, shape=(), dtype=float64)
nll_grad  cost time: 0.23255085945129395
tf.Tensor(-1010.44486771213, shape=(), dtype=float64)
nll_grad  cost time: 0.23274683952331543
tf.Tensor(-1010.5319860448744, shape=(), dtype=float64)
nll_grad  cost time: 0.23270416259765625
tf.Tensor(-1010.5395282283571, shape=(), dtype=float64)
nll_grad  cost time: 0.23105502128601074
tf.Tensor(-1010.5395284442584, shape=(), dtype=float64)
Optimization terminated successfully.
         Current function value: -1010.539528
         Iterations: 12
         Function evaluations: 17
         Gradient evaluations: 17
  message: Optimization terminated successfully.
  success: True
   status: 0
      fun: -1010.5395284442584
        x: [ 9.781e-01  1.531e+00  2.052e+00  1.629e+00  1.027e+00
            -2.223e+00]
      nit: 12
      jac: [ 1.300e-04 -8.077e-06  6.332e-05 -1.643e-05 -8.139e-05
             3.357e-05]
 hess_inv: [[ 3.869e-03 -2.088e-04 ...  2.808e-03  3.087e-04]
            [-2.088e-04  1.279e-02 ... -1.001e-04  5.875e-03]
            ...
            [ 2.808e-03 -1.001e-04 ...  5.526e-03  1.515e-03]
            [ 3.087e-04  5.875e-03 ...  1.515e-03  1.502e-02]]
     nfev: 17
     njev: 17
fit  cost time: 4.241307258605957
Using Model
Time for calculating errors: 0.5472018718719482
hesse_error: [0.061930489373256625, 0.11329283992581098, 0.11883471556261103, 0.09609044161580227, 0.07508951185378457, 0.12296377456412602]

we can see that thre fit results consistant with inputs, the first one is fixed.

for var in input_params:
    print(
        f"in: {input_params[var]} => out: {fit_result.params[var]} +/- {err.get(var, 0.)}"
    )

in: 6.0 => out: 6.0 +/- 0.0
in: 1.0 => out: 0.9780654317910398 +/- 0.061930489373256625
in: 2.0 => out: 2.0517922895579135 +/- 0.11883471556261103
in: 1.0 => out: 1.026998360907368 +/- 0.07508951185378457

We can use the amplitude to plot the fit results

amp = config.get_amplitude()

This is the total \(|A|^2\)

weight = amp(phsp)

This is the \(|A_i|^2\) for each decay chain

partial_weight = amp.partial_weight(phsp)

We can plot the data, Hist1D include some plot method base on matplotlib.

data_hist = Hist1D.histogram(
    data.get_mass("(C, D)"), bins=60, range=(0.25, 1.45)
)

Read the mass var from phsp.

mass_phsp = phsp.get_mass("(C, D)")

For helicity angle it can be read as phsp.get_angle("(C, D)", "C") The first arg is used to deteminate the decay chain. (decay chain include all particles) The second arg is used to deteminate the angle in decay chain. The return value is a dict as {"alpha": phi, "beta": theta, "gamma": 0}

phsp_hist = Hist1D.histogram(
    mass_phsp, weights=weight, bins=60, range=(0.25, 1.45)
)
scale = phsp_hist.scale_to(data_hist)

pw_hist = []
for w in partial_weight:
    # here we used more bins for a smooth plot
    hist = Hist1D.histogram(
        mass_phsp, weights=w, bins=60 * 2, range=(0.25, 1.45)
    )
    pw_hist.append(scale * hist * 2)

Then we can plot the histogram into matplotlib

for hist, dec in zip(pw_hist, config.get_decay()):
    hist.draw_kde(label=str(dec))
phsp_hist.draw(label="total amplitude")
data_hist.draw_error(label="toy data", color="black")

plt.legend()
plt.ylim((0, None))
plt.xlabel("M(CD)/ GeV")
plt.ylabel("Events/ 20 MeV")
plt.show()