Note
Go to the end to download the full example code
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.
22 config_str = """
23
24 decay:
25 A:
26 - [R1, B]
27 - [R2, C]
28 - [R3, D]
29 R1: [C, D]
30 R2: [B, D]
31 R3: [B, C]
32
33 particle:
34 $top:
35 A: { mass: 1.86, J: 0, P: -1}
36 $finals:
37 B: { mass: 0.494, J: 0, P: -1}
38 C: { mass: 0.139, J: 0, P: -1}
39 D: { mass: 0.139, J: 0, P: -1}
40 R1: [ R1_a, R1_b ]
41 R1_a: { mass: 0.7, width: 0.05, J: 1, P: -1}
42 R1_b: { mass: 0.5, width: 0.05, J: 0, P: +1}
43 R2: { mass: 0.824, width: 0.05, J: 0, P: +1}
44 R3: { mass: 0.824, width: 0.05, J: 0, P: +1}
45
46 """
The config file can be loaded by yaml
library.
52 import matplotlib.pyplot as plt
53 import yaml
54
55 from tf_pwa.config_loader import ConfigLoader
56 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.
64 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())
.
71 input_params = {
72 "A->R1_a.BR1_a->C.D_total_0r": 6.0,
73 "A->R1_b.BR1_b->C.D_total_0r": 1.0,
74 "A->R2.CR2->B.D_total_0r": 2.0,
75 "A->R3.DR3->B.C_total_0r": 1.0,
76 }
77 config.set_params(input_params)
True
Here we generate some toy data and phsp mc to show the model
82 data = config.generate_toy(1000)
83 phsp = config.generate_phsp(10000)
4.8%[▓▓>-----------------------------------------------] 0.50/10.44s eff: 90.000000%
115.7%[▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓>] 1.87/1.61s eff: 4.006541%
100.0%[▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓] 1.87/1.87s eff: 4.232456%
You can also fit the data fit to the data We can omit the args when writen in config.yml
90 fit_result = config.fit([data], [phsp])
91 # After the fit, you can get the uncertaities as
92 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.2627993916195042,
"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": -3.136050125543079,
"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": -3.0643983124211265,
"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(-1028.7814128788687, shape=(), dtype=float64)
nll_grad cost time: 0.27855825424194336
nll_grad cost time: 0.27373337745666504
nll_grad cost time: 0.27289628982543945
tf.Tensor(-1029.9584106684333, shape=(), dtype=float64)
nll_grad cost time: 0.27263331413269043
tf.Tensor(-1031.1411301462558, shape=(), dtype=float64)
nll_grad cost time: 0.2725028991699219
tf.Tensor(-1032.5227180749043, shape=(), dtype=float64)
nll_grad cost time: 0.27374267578125
nll_grad cost time: 0.2745954990386963
tf.Tensor(-1032.672551646443, shape=(), dtype=float64)
nll_grad cost time: 0.2795393466949463
tf.Tensor(-1032.9239577790458, shape=(), dtype=float64)
nll_grad cost time: 0.2728538513183594
nll_grad cost time: 0.27350568771362305
tf.Tensor(-1034.8304560904971, shape=(), dtype=float64)
nll_grad cost time: 0.27144527435302734
nll_grad cost time: 0.2762782573699951
tf.Tensor(-1035.0773404217034, shape=(), dtype=float64)
nll_grad cost time: 0.2789640426635742
tf.Tensor(-1035.4878576544343, shape=(), dtype=float64)
nll_grad cost time: 0.28232526779174805
tf.Tensor(-1035.7922544633993, shape=(), dtype=float64)
nll_grad cost time: 0.28502321243286133
tf.Tensor(-1035.792857235063, shape=(), dtype=float64)
nll_grad cost time: 0.2757582664489746
tf.Tensor(-1035.7928696328354, shape=(), dtype=float64)
nll_grad cost time: 0.27629780769348145
tf.Tensor(-1035.792869796167, shape=(), dtype=float64)
Optimization terminated successfully.
Current function value: -1035.792870
Iterations: 12
Function evaluations: 17
Gradient evaluations: 17
message: Optimization terminated successfully.
success: True
status: 0
fun: -1035.792869796167
x: [ 8.291e-01 1.186e+00 1.868e+00 -3.125e+00 8.889e-01
-3.145e+00]
nit: 12
jac: [-3.099e-04 5.205e-04 7.961e-05 8.374e-05 -3.309e-04
-6.489e-04]
hess_inv: [[ 2.288e-03 1.490e-04 ... 1.545e-03 -2.191e-04]
[ 1.490e-04 1.633e-02 ... -5.497e-04 2.498e-03]
...
[ 1.545e-03 -5.497e-04 ... 3.938e-03 6.264e-04]
[-2.191e-04 2.498e-03 ... 6.264e-04 2.474e-02]]
nfev: 17
njev: 17
fit cost time: 4.996109485626221
Using Model
Time for calculating errors: 0.6261942386627197
hesse_error: [0.0479527781751716, 0.1331269396665368, 0.10405612923910544, 0.11184361809385386, 0.06114407675528138, 0.17638873219372456]
we can see that thre fit results consistant with inputs, the first one is fixed.
97 for var in input_params:
98 print(
99 f"in: {input_params[var]} => out: {fit_result.params[var]} +/- {err.get(var, 0.)}"
100 )
in: 6.0 => out: 6.0 +/- 0.0
in: 1.0 => out: 0.8291479160521567 +/- 0.0479527781751716
in: 2.0 => out: 1.868446499951094 +/- 0.10405612923910544
in: 1.0 => out: 0.8888865050046114 +/- 0.06114407675528138
We can use the amplitude to plot the fit results
105 amp = config.get_amplitude()
This is the total \(|A|^2\)
110 weight = amp(phsp)
This is the \(|A_i|^2\) for each decay chain
115 partial_weight = amp.partial_weight(phsp)
We can plot the data, Hist1D include some plot method base on matplotlib.
120 data_hist = Hist1D.histogram(
121 data.get_mass("(C, D)"), bins=60, range=(0.25, 1.45)
122 )
Read the mass var from phsp.
127 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}
136 phsp_hist = Hist1D.histogram(
137 mass_phsp, weights=weight, bins=60, range=(0.25, 1.45)
138 )
139 scale = phsp_hist.scale_to(data_hist)
140
141 pw_hist = []
142 for w in partial_weight:
143 # here we used more bins for a smooth plot
144 hist = Hist1D.histogram(
145 mass_phsp, weights=w, bins=60 * 2, range=(0.25, 1.45)
146 )
147 pw_hist.append(scale * hist * 2)
Then we can plot the histogram into matplotlib
152 for hist, dec in zip(pw_hist, config.get_decay()):
153 hist.draw_kde(label=str(dec))
154 phsp_hist.draw(label="total amplitude")
155 data_hist.draw_error(label="toy data", color="black")
156
157 plt.legend()
158 plt.ylim((0, None))
159 plt.xlabel("M(CD)/ GeV")
160 plt.ylabel("Events/ 20 MeV")
161 plt.show()
Total running time of the script: (0 minutes 8.763 seconds)