Fitting Tutorial

In [1]:
import PyPWA as pwa
import pandas
import numpy as npy
from IPython.display import display
import warnings
warnings.filterwarnings('ignore')

Define Waves for Fit (and input initial values of minuit and fitted parameters).*

In this example (as expected by the defined amplitude)\ each wave is defined by (epsilon.l.m) and each parameter has a real and imaginary part.\ i.e a epsilon=-1, l=1 (P wave), m=1 will produce Vs(r.-1.1.1) and Vs(i.-1.1.1) names.\ (In this example the imaginary part of the P-wave is kept fixed at 0 value\ in the fit)

In [2]:
Vs = {"errordef": 1}
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")
#    initial.append(f"{param_type}.1.3.1")
#    initial.append(f"{param_type}.1.4.1")
    
    Vs[f"{param_type}.1.0.0"] = 10
    Vs[f"limit_{param_type}.1.0.0"] = [-500, 1500.]
    Vs[f"error_{param_type}.1.0.0"] = .1
#    Vs[f"fix_i.-1.0.0"] = True

    Vs[f"{param_type}.1.1.0"] = 10
    Vs[f"limit_{param_type}.1.1.0"] = [-500, 1500.]
    Vs[f"error_{param_type}.1.1.0"] = .1

    Vs[f"{param_type}.1.1.1"] = 10
    Vs[f"limit_{param_type}.1.1.1"] = [-500, 1500.]
    Vs[f"error_{param_type}.1.1.1"] = .1
    Vs[f"fix_i.1.1.1"] = True
        
    Vs[f"{param_type}.1.2.0"] = 10
    Vs[f"limit_{param_type}.1.2.0"] = [-500, 1500.]
    Vs[f"error_{param_type}.1.2.0"] = .1
    
    Vs[f"{param_type}.1.2.1"] = 10
    Vs[f"limit_{param_type}.1.2.1"] = [-500, 1500.]
    Vs[f"error_{param_type}.1.2.1"] = .1
       
    Vs[f"{param_type}.1.2.2"] = 10
    Vs[f"limit_{param_type}.1.2.2"] = [-500, 1500.]
    Vs[f"error_{param_type}.1.2.2"] = .1 
    
#    Vs[f"{param_type}.1.3.1"] = 10
#    Vs[f"limit_{param_type}.1.3.1"] = [-500, 1500.]
#    Vs[f"error_{param_type}.1.3.1"] = .1
    
#    Vs[f"{param_type}.1.4.1"] = 10
#    Vs[f"limit_{param_type}.1.4.1"] = [-500, 1500.]
#    Vs[f"error_{param_type}.1.4.1"] = .1
    
Vs[f"i.1.1.1"] = 0.1

Read data and montecarlos (accepted and generated) samples

In [3]:
datasample =  pandas.read_csv("simdata_JPAC.csv")
#accmcsample = pandas.read_csv("simdata5.csv")
#rawmcsample = pandas.read_csv("simdata5.csv")
#datasample = pwa.read("etapi_data_data.txt")
#accmcsample = pwa.read("etapi_acc.txt")
accmcsample = pwa.read("../TUTORIAL_FILES/etapi_acc.txt")
rawmcsample = pwa.read("../TUTORIAL_FILES/etapi_raw.txt")

Binning of the data/monte-carlo and define amplitude (function) to fit*

Here the user difine number of bins, variable to be binned and range

In [4]:
#import AmplitudeOLDfit
#amplitude = AmplitudeOLDfit.FitAmplitude(initial)
import AmplitudeJPACfit
amplitude = AmplitudeJPACfit.FitAmplitude(initial)
#Define number of bins 
nbins = 20
binsda = pwa.bin_by_range(datasample, "mass", nbins, .6, 2.0)
binsma = pwa.bin_by_range(accmcsample, accmcsample["mass"], nbins, .6, 2.0)
binsmr = pwa.bin_by_range(rawmcsample, rawmcsample["mass"], nbins, .6, 2.0)

Check that bins have enough number of events for fit

In [5]:
for bin in binsda:
    print(len(bin))
4997
6461
10341
25358
66333
27471
14968
17621
32136
45246
19835
10604
10281
12465
15419
18554
17562
12538
8345
5550

Fitting with Minuit and Extended LogLikelihood

Look at other possibilities through pypwa (use the ?pwa command\ or see https://pypwa.jlab.org or https://github.com/JeffersonLab/PyPWA)

In [6]:
results = []
for index, dbin in enumerate(binsda):
    with pwa.LogLikelihood(
        amplitude, dbin, binsma[index], generated_length=len(binsmr[index])) as Likelihood:
        results.append(
            pwa.minuit(initial, Vs, Likelihood, 1, 2, 2000)
        )

Looking at the results of fitting

In [7]:
failed = 0
for index, result in enumerate(results):
    display(f"Bin #{index}:")
    display(result.get_fmin())
    try:
        display(result.get_param_states())
    except:
        failed += 1
        pass

display(f"{(failed / len(results)) * 100}% ({failed})Failed. ")
'Bin #0:'
FCN = -4.132E+04 Ncalls = 1116 (1116 total)
EDM = 1.55E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 101 7 -500 1.5E+03
1 r.1.1.0 65 7 -500 1.5E+03
2 r.1.1.1 61 4 -500 1.5E+03
3 r.1.2.0 55 6 -500 1.5E+03
4 r.1.2.1 66 4 -500 1.5E+03
5 r.1.2.2 64.0 3.1 -500 1.5E+03
6 i.1.0.0 26 24 -500 1.5E+03
7 i.1.1.0 2.3 17.0 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 32 14 -500 1.5E+03
10 i.1.2.1 -2.2 11.3 -500 1.5E+03
11 i.1.2.2 10 10 -500 1.5E+03
'Bin #1:'
FCN = -5.47E+04 Ncalls = 1432 (1432 total)
EDM = 0.0002 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 140 4 -500 1.5E+03
1 r.1.1.0 62.0 2.6 -500 1.5E+03
2 r.1.1.1 62.2 2.2 -500 1.5E+03
3 r.1.2.0 64 5 -500 1.5E+03
4 r.1.2.1 67.6 2.2 -500 1.5E+03
5 r.1.2.2 68.7 2.9 -500 1.5E+03
6 i.1.0.0 7 26 -500 1.5E+03
7 i.1.1.0 3.0 15.8 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 17 25 -500 1.5E+03
10 i.1.2.1 2.1 10.7 -500 1.5E+03
11 i.1.2.2 -6 16 -500 1.5E+03
'Bin #2:'
FCN = -9.091E+04 Ncalls = 1957 (1957 total)
EDM = 0.000167 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 187 17 -500 1.5E+03
1 r.1.1.0 70 9 -500 1.5E+03
2 r.1.1.1 70 5 -500 1.5E+03
3 r.1.2.0 62 16 -500 1.5E+03
4 r.1.2.1 58 7 -500 1.5E+03
5 r.1.2.2 69 6 -500 1.5E+03
6 i.1.0.0 -60 50 -500 1.5E+03
7 i.1.1.0 -14 28 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -49 27 -500 1.5E+03
10 i.1.2.1 -40 9 -500 1.5E+03
11 i.1.2.2 -25 13 -500 1.5E+03
'Bin #3:'
FCN = -2.385E+05 Ncalls = 1924 (1924 total)
EDM = 0.000104 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 323 17 -500 1.5E+03
1 r.1.1.0 76 6 -500 1.5E+03
2 r.1.1.1 70 4 -500 1.5E+03
3 r.1.2.0 80 5 -500 1.5E+03
4 r.1.2.1 89 6 -500 1.5E+03
5 r.1.2.2 85 6 -500 1.5E+03
6 i.1.0.0 140 30 -500 1.5E+03
7 i.1.1.0 -11 12 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -12 14 -500 1.5E+03
10 i.1.2.1 -31 16 -500 1.5E+03
11 i.1.2.2 -8 12 -500 1.5E+03
'Bin #4:'
FCN = -6.737E+05 Ncalls = 1804 (1804 total)
EDM = 0.000157 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 270 50 -500 1.5E+03
1 r.1.1.0 47 9 -500 1.5E+03
2 r.1.1.1 86 12 -500 1.5E+03
3 r.1.2.0 95 10 -500 1.5E+03
4 r.1.2.1 97 11 -500 1.5E+03
5 r.1.2.2 90 14 -500 1.5E+03
6 i.1.0.0 554 20 -500 1.5E+03
7 i.1.1.0 13 6 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -26 10 -500 1.5E+03
10 i.1.2.1 -30 11 -500 1.5E+03
11 i.1.2.2 -22 9 -500 1.5E+03
'Bin #5:'
FCN = -2.561E+05 Ncalls = 1466 (1466 total)
EDM = 0.000174 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 -220 120 -500 1.5E+03
1 r.1.1.0 80 11 -500 1.5E+03
2 r.1.1.1 59 31 -500 1.5E+03
3 r.1.2.0 99 6 -500 1.5E+03
4 r.1.2.1 108 18 -500 1.5E+03
5 r.1.2.2 111 11 -500 1.5E+03
6 i.1.0.0 280 100 -500 1.5E+03
7 i.1.1.0 20 41 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -9 38 -500 1.5E+03
10 i.1.2.1 -7 57 -500 1.5E+03
11 i.1.2.2 4 56 -500 1.5E+03
'Bin #6:'
FCN = -1.328E+05 Ncalls = 1748 (1748 total)
EDM = 9.32E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 -150 6 -500 1.5E+03
1 r.1.1.0 85 7 -500 1.5E+03
2 r.1.1.1 69 5 -500 1.5E+03
3 r.1.2.0 129 3 -500 1.5E+03
4 r.1.2.1 129 4 -500 1.5E+03
5 r.1.2.2 127 3 -500 1.5E+03
6 i.1.0.0 -77 15 -500 1.5E+03
7 i.1.1.0 -57 10 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 18 12 -500 1.5E+03
10 i.1.2.1 -25 15 -500 1.5E+03
11 i.1.2.2 -17 14 -500 1.5E+03
'Bin #7:'
FCN = -1.604E+05 Ncalls = 1468 (1468 total)
EDM = 6.4E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 -74 10 -500 1.5E+03
1 r.1.1.0 77 7 -500 1.5E+03
2 r.1.1.1 104 5 -500 1.5E+03
3 r.1.2.0 139 8 -500 1.5E+03
4 r.1.2.1 156 8 -500 1.5E+03
5 r.1.2.2 137 11 -500 1.5E+03
6 i.1.0.0 -64 14 -500 1.5E+03
7 i.1.1.0 68 8 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 88 12 -500 1.5E+03
10 i.1.2.1 53 18 -500 1.5E+03
11 i.1.2.2 84 18 -500 1.5E+03
'Bin #8:'
FCN = -3.129E+05 Ncalls = 1605 (1605 total)
EDM = 2.36E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 -58 11 -500 1.5E+03
1 r.1.1.0 86 5 -500 1.5E+03
2 r.1.1.1 121 5 -500 1.5E+03
3 r.1.2.0 216 9 -500 1.5E+03
4 r.1.2.1 208 9 -500 1.5E+03
5 r.1.2.2 222 11 -500 1.5E+03
6 i.1.0.0 1.5 24.6 -500 1.5E+03
7 i.1.1.0 61 8 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 103 18 -500 1.5E+03
10 i.1.2.1 118 14 -500 1.5E+03
11 i.1.2.2 94 27 -500 1.5E+03
'Bin #9:'
FCN = -4.551E+05 Ncalls = 1702 (1702 total)
EDM = 0.000401 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 25 7 -500 1.5E+03
1 r.1.1.0 -70 5 -500 1.5E+03
2 r.1.1.1 -102 9 -500 1.5E+03
3 r.1.2.0 -279 13 -500 1.5E+03
4 r.1.2.1 -219 18 -500 1.5E+03
5 r.1.2.2 -261 11 -500 1.5E+03
6 i.1.0.0 65 15 -500 1.5E+03
7 i.1.1.0 34 9 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 107 29 -500 1.5E+03
10 i.1.2.1 208 17 -500 1.5E+03
11 i.1.2.2 121 25 -500 1.5E+03
'Bin #10:'
FCN = -1.82E+05 Ncalls = 2015 (2015 total)
EDM = 1.77 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
False True False True
Hesse failed Has cov. Accurate Pos. def. Forced
False True False True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 -6 7 -500 1.5E+03
1 r.1.1.0 58 5 -500 1.5E+03
2 r.1.1.1 76 6 -500 1.5E+03
3 r.1.2.0 179.9 2.9 -500 1.5E+03
4 r.1.2.1 133 6 -500 1.5E+03
5 r.1.2.2 148 3 -500 1.5E+03
6 i.1.0.0 -145.7 2.9 -500 1.5E+03
7 i.1.1.0 49.9 3.1 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -35 8 -500 1.5E+03
10 i.1.2.1 -132 6 -500 1.5E+03
11 i.1.2.2 -39 7 -500 1.5E+03
'Bin #11:'
FCN = -9.038E+04 Ncalls = 2018 (2018 total)
EDM = 3.18 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
False True False True
Hesse failed Has cov. Accurate Pos. def. Forced
False True False True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 -54 4 -500 1.5E+03
1 r.1.1.0 53.2 3.0 -500 1.5E+03
2 r.1.1.1 59 5 -500 1.5E+03
3 r.1.2.0 53 3 -500 1.5E+03
4 r.1.2.1 46 4 -500 1.5E+03
5 r.1.2.2 93 5 -500 1.5E+03
6 i.1.0.0 161.8 3.1 -500 1.5E+03
7 i.1.1.0 27.6 2.6 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 61 4 -500 1.5E+03
10 i.1.2.1 95.4 2.7 -500 1.5E+03
11 i.1.2.2 58 3 -500 1.5E+03
'Bin #12:'
FCN = -8.802E+04 Ncalls = 1019 (1019 total)
EDM = 2.91E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 -2.4 4.7 -500 1.5E+03
1 r.1.1.0 154 6 -500 1.5E+03
2 r.1.1.1 165 4 -500 1.5E+03
3 r.1.2.0 59 5 -500 1.5E+03
4 r.1.2.1 59 4 -500 1.5E+03
5 r.1.2.2 55 3 -500 1.5E+03
6 i.1.0.0 -39 18 -500 1.5E+03
7 i.1.1.0 21 16 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -29 10 -500 1.5E+03
10 i.1.2.1 -25 10 -500 1.5E+03
11 i.1.2.2 -15 15 -500 1.5E+03
'Bin #13:'
FCN = -1.1E+05 Ncalls = 1396 (1396 total)
EDM = 1.16E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 3.2 4.1 -500 1.5E+03
1 r.1.1.0 171 4 -500 1.5E+03
2 r.1.1.1 163 7 -500 1.5E+03
3 r.1.2.0 74.6 2.6 -500 1.5E+03
4 r.1.2.1 77 4 -500 1.5E+03
5 r.1.2.2 77 5 -500 1.5E+03
6 i.1.0.0 -64 22 -500 1.5E+03
7 i.1.1.0 -8 13 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 28 17 -500 1.5E+03
10 i.1.2.1 -12 24 -500 1.5E+03
11 i.1.2.2 15 10 -500 1.5E+03
'Bin #14:'
FCN = -1.403E+05 Ncalls = 1260 (1260 total)
EDM = 0.000115 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 9 6 -500 1.5E+03
1 r.1.1.0 143 4 -500 1.5E+03
2 r.1.1.1 156 5 -500 1.5E+03
3 r.1.2.0 119 3 -500 1.5E+03
4 r.1.2.1 118 4 -500 1.5E+03
5 r.1.2.2 118 5 -500 1.5E+03
6 i.1.0.0 -93 8 -500 1.5E+03
7 i.1.1.0 6 11 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 0.7 9.8 -500 1.5E+03
10 i.1.2.1 -30 6 -500 1.5E+03
11 i.1.2.2 18 10 -500 1.5E+03
'Bin #15:'
FCN = -1.733E+05 Ncalls = 1766 (1766 total)
EDM = 5.28E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 27 6 -500 1.5E+03
1 r.1.1.0 135.7 3.0 -500 1.5E+03
2 r.1.1.1 149 4 -500 1.5E+03
3 r.1.2.0 151 4 -500 1.5E+03
4 r.1.2.1 142 5 -500 1.5E+03
5 r.1.2.2 140 7 -500 1.5E+03
6 i.1.0.0 36 10 -500 1.5E+03
7 i.1.1.0 -17 8 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -32 13 -500 1.5E+03
10 i.1.2.1 -63 9 -500 1.5E+03
11 i.1.2.2 -75 13 -500 1.5E+03
'Bin #16:'
FCN = -1.633E+05 Ncalls = 1452 (1452 total)
EDM = 0.000148 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 5.0 2.8 -500 1.5E+03
1 r.1.1.0 125 4 -500 1.5E+03
2 r.1.1.1 131 4 -500 1.5E+03
3 r.1.2.0 155 4 -500 1.5E+03
4 r.1.2.1 151 4 -500 1.5E+03
5 r.1.2.2 159 3 -500 1.5E+03
6 i.1.0.0 29 23 -500 1.5E+03
7 i.1.1.0 33 11 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 26 24 -500 1.5E+03
10 i.1.2.1 39 11 -500 1.5E+03
11 i.1.2.2 13 27 -500 1.5E+03
'Bin #17:'
FCN = -1.129E+05 Ncalls = 1474 (1474 total)
EDM = 0.000107 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 18 3 -500 1.5E+03
1 r.1.1.0 114 4 -500 1.5E+03
2 r.1.1.1 121 4 -500 1.5E+03
3 r.1.2.0 128 3 -500 1.5E+03
4 r.1.2.1 127 4 -500 1.5E+03
5 r.1.2.2 128.9 3.1 -500 1.5E+03
6 i.1.0.0 14 19 -500 1.5E+03
7 i.1.1.0 6 18 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 14 20 -500 1.5E+03
10 i.1.2.1 35 9 -500 1.5E+03
11 i.1.2.2 10 17 -500 1.5E+03
'Bin #18:'
FCN = -7.181E+04 Ncalls = 1275 (1275 total)
EDM = 2.9E-06 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 15.9 2.9 -500 1.5E+03
1 r.1.1.0 99 4 -500 1.5E+03
2 r.1.1.1 98.0 3.1 -500 1.5E+03
3 r.1.2.0 105.1 2.5 -500 1.5E+03
4 r.1.2.1 104.7 3.1 -500 1.5E+03
5 r.1.2.2 104.2 2.7 -500 1.5E+03
6 i.1.0.0 9 13 -500 1.5E+03
7 i.1.1.0 -20 14 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -6 19 -500 1.5E+03
10 i.1.2.1 5 12 -500 1.5E+03
11 i.1.2.2 2.4 19.9 -500 1.5E+03
'Bin #19:'
FCN = -4.572E+04 Ncalls = 1529 (1529 total)
EDM = 7.67E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 17.4 3.1 -500 1.5E+03
1 r.1.1.0 84.9 2.7 -500 1.5E+03
2 r.1.1.1 82 3 -500 1.5E+03
3 r.1.2.0 77 3 -500 1.5E+03
4 r.1.2.1 86 3 -500 1.5E+03
5 r.1.2.2 83 4 -500 1.5E+03
6 i.1.0.0 0.22 12.14 -500 1.5E+03
7 i.1.1.0 -1.3 12.1 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -12 13 -500 1.5E+03
10 i.1.2.1 -12 10 -500 1.5E+03
11 i.1.2.2 -22 15 -500 1.5E+03
'0.0% (0)Failed. '

Checking the waves used (and filling waves variable for later use)

In [18]:
waves = amplitude.make_elm(result.values)
#waves = AmplitudeJPACfit.FitAmplitude.make_elm(result.values)
waves
Out[18]:
e l m
0 1 0 0
1 1 1 0
2 1 1 1
3 1 2 0
4 1 2 1
5 1 2 2

Filling strings with wave names (for plotting)

In [19]:
wave = npy.empty(len(waves),dtype="int")
string = npy.empty(len(waves),dtype="U5")
for index, w in waves.iterrows():
    #print(w["l"])
    string[index]= "{}{}{}".format(w["e"],w["l"],w["m"])
    wave[index]=w.name

Getting the bin mass values and number of events in datasample for those bins

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

Calculating the expected number of events in a a mass bin

In [21]:
total_nExp = npy.empty(len(binsma))
for index, the_bin, result in zip(range(len(binsma)), binsma, results):
#    amp = testAmplitudeFit.FitAmplitude(initial)
#    amp = AmplitudeJPACfit.FitAmplitude(initial)
    amp=amplitude
    amp.setup(the_bin)
    total_nExp[index] = npy.average(amp.calculate(result.values))

Plot expected number of events vs mass and data vs mass (both should agree if fitting worked)

In [22]:
import matplotlib.pyplot as plt
mni = npy.empty(len(total_nExp), dtype=[("mass", float), ("int", float)])
mda = npy.empty(len(mcounts), dtype=[("mass", float), ("intd", float)])
mni["mass"] = bmass
mni["int"] = total_nExp
mni = pandas.DataFrame(mni)
counts, bin_edges = npy.histogram(mni["mass"], nbins, weights=mni["int"])
mda["mass"] = bmass
mda["intd"] = mcounts
mda = pandas.DataFrame(mda)
dcounts, bin_edges = npy.histogram(mda["mass"], nbins, weights=mda["intd"])
#dcounts, bin_edges = npy.histogram(bmass, nbins, weights=total_nExp)
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)
myerr = npy.empty(nbins)
myerr = npy.sqrt(dcounts)
#plt.errorbar(centers,counts, yerr, fmt="s",linestyle="dashed",markersize='10',label="nExp")
plt.errorbar(centers,total_nExp, yerr, fmt="s",linestyle="dashed",markersize='10',label="nExp")
plt.errorbar(centers,mcounts, myerr, fmt="o",label="Data")
#plt.xlim(.6, 3.2)
#plt.ylim(0.,65000)
plt.legend(loc='upper right')
plt.show()

Calculate initial intensities (in case we need to check them)

In [13]:
intensities = []
for the_bin in binsda:
#    amp = AmplitudeJPACfit.FitAmplitude(initial)
    amp.setup(the_bin)
    intensities.append(amp.calculate(Vs))
intesities=pandas.DataFrame(intensities)

Calculate the expected number of events for each wave (vs mass)

In [14]:
wave_nExp = npy.empty([len(waves),len(binsma)],"f16")
for i in range(len(waves)):
    for index, the_bin, result in zip(range(len(binsma)), binsma, results):
#        amp = testAmplitudeFit.FitAmplitude(initial)
        amp.setup(the_bin)
        wave_nExp[i][index] = npy.average(amp.calculate_wave(result.values,wave[i]))
In [ ]:
results.

Plot expected number of events vs mass for each wave

In [15]:
for i in range(len(waves)):
    mni0 = npy.empty(len(wave_nExp[i]), dtype=[("mass", float), ("int0", float)])
    mni0["mass"] = bmass
    mni0["int0"] = wave_nExp[i]
    mni0 = pandas.DataFrame(mni0)
    counts0, bin_edges = npy.histogram(mni0["mass"], nbins, weights=mni0["int0"])
    centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# Add yerr to argment list when we have errors
    yerr = npy.empty(nbins)
    yerr = npy.sqrt(counts0)
    plt.errorbar(centers,wave_nExp[i], yerr, linestyle="dashed", fmt="o",label=string[i])
   
    plt.legend(loc='upper right', title="ε L M")
    plt.show()

Expected number of events vs mass for each wave in a same plot

In [16]:
for i in range(len(waves)):
    plt.errorbar(centers,wave_nExp[i], yerr, linestyle="dashed",fmt="o",label=string[i])
    plt.legend(loc='upper right', title="ε L M")
plt.show()

Calculate PhaseMotion between two waves

In this example between 2nd and 3rd waves in amp list

In [17]:
phase = npy.empty(len(binsda))
for index, the_bin, result in zip(range(len(binsda)), binsda, results):
#    amp = testAmplitudeFit.FitAmplitude()
#    amp = AmplitudeJPACfit.FitAmplitude(initial)
    amp.setup(the_bin)
    phase[index] = amp.Phasediff(result.values,wave[4],wave[2])

Plot PhaseMotion

In [18]:
mnip = npy.empty(len(bmass), dtype=[("mass", float), ("intp", float)])
mnip["mass"] = bmass
mnip["intp"] = phase
mnip = pandas.DataFrame(mnip)
countsp, bin_edges = npy.histogram(mnip["mass"], nbins, weights=mnip["intp"])
#countsp, bin_edges = npy.histogram(mnip["mass"], 25)

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

# Add yerr to argment list when we have errors
yerr = npy.empty(nbins)
yerr = npy.sqrt(npy.abs(countsp))
plt.errorbar(centers,countsp, yerr, fmt="o")
#plt.xlim(.6, 2.5)
Out[18]:
<ErrorbarContainer object of 3 artists>

Write fitted values (of production amplitudes) to disk

In [48]:
final_values = pandas.DataFrame([res.np_values() for res in results], columns=results[0].parameters)
final_values.to_csv("final_values_JPAC.csv", index=False)
In [49]:
pandas.DataFrame([res.np_errors() for res in results], columns=results[0].parameters)
Out[49]:
r.1.0.0 r.1.1.0 r.1.1.1 r.1.2.0 r.1.2.1 r.1.2.2 i.1.0.0 i.1.1.0 i.1.1.1 i.1.2.0 i.1.2.1 i.1.2.2
0 6.954829 6.827580 3.554664 6.018165 3.976862 3.121498 23.742340 17.046423 0.1 13.907862 11.264237 10.455112
1 3.802987 2.606027 2.155161 4.931763 2.201047 2.864100 25.674425 15.797999 0.1 24.915356 10.674988 15.745680
2 17.040772 9.297779 4.566329 16.217777 7.233109 6.223982 51.559825 27.907861 0.1 26.606110 9.219745 12.651604
3 16.663525 5.503232 3.681815 5.228695 6.205652 6.286110 34.349073 12.080370 0.1 14.208845 16.395569 12.242415
4 45.232868 9.477409 11.742453 9.604031 10.595964 14.102345 20.092418 5.964615 0.1 10.151890 10.604341 8.624453
5 117.579745 11.084983 30.532282 5.623476 17.876051 10.854289 95.148666 41.011584 0.1 37.507708 56.539039 56.141556
6 6.388357 6.741500 4.536766 3.334645 4.100674 3.270934 14.559920 9.978258 0.1 11.849006 15.382238 14.071663
7 9.674936 6.591454 4.655869 7.911636 7.514837 10.782131 13.518799 8.166621 0.1 12.071510 17.608551 17.748484
8 10.535564 5.193630 5.196341 8.722004 9.092910 11.320792 24.556608 8.433569 0.1 18.478807 13.644900 27.333611
9 7.006934 5.265910 9.043399 12.551561 18.423677 10.526510 15.301628 9.072606 0.1 29.455452 16.838985 24.953206
10 7.030241 5.426042 5.847165 2.896373 6.227190 3.494455 2.915719 3.129739 0.1 8.321822 6.009766 6.825844
11 4.401767 3.046399 4.518873 3.243339 3.723411 4.818405 3.069889 2.562528 0.1 3.653762 2.717735 3.232582
12 4.702111 6.437627 4.022266 5.199858 3.658756 3.321823 18.390600 16.388502 0.1 10.462186 9.566562 15.331334
13 4.110845 4.309469 6.863718 2.640409 3.599147 4.509742 22.200358 13.452693 0.1 16.926294 23.687215 10.220994
14 6.242526 3.706594 5.056283 3.352333 4.192237 4.599247 8.376552 11.073956 0.1 9.842125 5.801195 10.487145
15 6.010469 2.994718 3.866779 3.843395 5.182108 6.924436 9.784828 8.492263 0.1 13.240668 9.010339 12.520236
16 2.827927 3.697508 4.226044 4.084553 4.389821 3.454387 23.199488 11.077514 0.1 24.427168 11.324746 27.486141
17 3.174387 3.598044 3.731598 3.476954 4.313309 3.050902 18.547449 18.347815 0.1 20.004608 8.999662 17.015956
18 2.914699 4.156095 3.146181 2.450015 3.084025 2.706066 13.419362 14.321449 0.1 18.908542 12.106658 19.887643
19 3.146650 2.677247 3.362296 3.336464 3.461697 4.229378 12.136651 12.102438 0.1 12.635719 9.602312 15.353926
In [50]:
final_values
Out[50]:
r.1.0.0 r.1.1.0 r.1.1.1 r.1.2.0 r.1.2.1 r.1.2.2 i.1.0.0 i.1.1.0 i.1.1.1 i.1.2.0 i.1.2.1 i.1.2.2
0 101.082208 65.138042 61.001270 55.000546 66.110221 64.039694 26.032063 2.337883 0.1 31.674827 -2.213773 9.942500
1 139.674858 61.971084 62.188751 64.236630 67.568739 68.739345 6.852412 2.975524 0.1 17.437169 2.096012 -5.689982
2 187.348631 70.044407 70.450004 62.156223 58.328096 68.753584 -59.471006 -14.020646 0.1 -48.587415 -39.711547 -24.980981
3 323.457581 76.494688 69.871990 80.312123 88.685671 85.283301 143.034647 -10.904887 0.1 -11.542512 -30.675265 -7.934322
4 265.310200 46.773274 85.863858 94.680866 96.586948 90.222383 554.233154 13.135420 0.1 -26.423956 -30.397467 -21.938487
5 -222.228190 79.661793 59.414017 98.521007 107.690304 111.036456 281.358882 19.787339 0.1 -9.085648 -6.783025 4.398715
6 -150.295619 84.587045 69.472345 129.448897 129.034043 127.254223 -77.090868 -56.792958 0.1 18.232602 -25.402967 -17.365946
7 -73.507421 76.685208 104.469047 139.207053 156.021989 137.310467 -63.814576 67.994759 0.1 88.176964 53.289035 83.844111
8 -57.946405 85.706192 121.441780 216.396291 208.121464 222.494438 1.525251 61.228584 0.1 102.984406 118.347510 94.222542
9 25.126468 -70.063300 -101.841991 -279.395238 -219.050288 -261.424226 65.441440 34.421981 0.1 106.689361 208.475223 120.543431
10 -5.811603 57.525949 75.943609 179.925780 132.889785 147.619872 -145.684550 49.852821 0.1 -34.651359 -132.290765 -39.137939
11 -54.484340 53.239005 59.471448 52.810897 46.002772 93.003283 161.769891 27.600019 0.1 61.198545 95.424898 58.143936
12 -2.425633 154.077802 164.973803 59.395493 59.354709 55.118071 -39.446433 21.124348 0.1 -28.952056 -24.579618 -14.737104
13 3.159826 170.783543 163.281327 74.579066 77.026218 77.120156 -63.803290 -7.694733 0.1 27.872617 -12.342172 15.304048
14 8.506245 143.246131 155.557679 118.778850 117.810134 117.931485 -93.346335 5.874608 0.1 0.669072 -30.246420 17.726145
15 27.033009 135.679204 149.111848 150.647481 142.387799 139.996493 36.428653 -17.248620 0.1 -32.094729 -62.536027 -74.887932
16 5.039317 124.972595 131.434780 154.959420 150.987142 159.123139 29.221314 33.219667 0.1 26.463682 39.182434 13.314298
17 18.333928 113.537124 121.398279 128.300397 126.713200 128.937173 13.624954 6.056093 0.1 13.505022 35.103369 9.816295
18 15.935776 98.595268 98.049081 105.128338 104.687480 104.178499 9.308810 -19.591011 0.1 -6.404105 5.487549 2.444948
19 17.409842 84.880760 81.797274 77.363298 85.572313 82.639824 0.223294 -1.307494 0.1 -12.341157 -11.921014 -21.776772
In [ ]: