Fitting Tutorial

[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) ==========================================================================================================================================

[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

[3]:
datasample =  pandas.read_csv("simdata_JPAC-np.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

[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

[5]:
for bin in binsda:
    print(len(bin))
6577
9364
16674
45953
127024
48129
21406
23014
41661
58948
20078
5179
3404
5434
11026
18505
19504
13486
8425
5481

Fitting with Minuit and Extended LogLikelihood

Look at other possibilities through pypwa (use the ?pwa command
[6]:
from IPython.display import display
intensities = []
for the_bin in binsda:
#    amp = AmplitudeJPACfit.FitAmplitude(initial)
    amplitude.setup(the_bin)
#    intensities.append(amplitude.calculate(Vs))
    display(pandas.DataFrame(amplitude.calculate(Vs)))
0
0 524.818443
1 112.544543
2 75.976211
3 64.049006
4 143.808291
... ...
6572 105.370331
6573 78.846342
6574 120.701116
6575 1477.089631
6576 123.955803

6577 rows × 1 columns

0
0 86.996618
1 447.912162
2 161.396588
3 759.107477
4 471.011055
... ...
9359 161.296093
9360 68.394902
9361 120.358836
9362 146.710380
9363 1196.737830

9364 rows × 1 columns

0
0 102.392089
1 595.649766
2 65.344248
3 280.705037
4 602.855631
... ...
16669 67.959747
16670 124.297777
16671 301.015157
16672 456.641653
16673 28.929693

16674 rows × 1 columns

0
0 4.379040
1 29.951455
2 7.905595
3 7.029131
4 710.321421
... ...
45948 10.878720
45949 945.802674
45950 206.911768
45951 179.314430
45952 9.352145

45953 rows × 1 columns

0
0 45.201774
1 43.788671
2 68.660084
3 9.205558
4 114.399723
... ...
127019 126.680693
127020 1049.419573
127021 51.796865
127022 157.035089
127023 38.986098

127024 rows × 1 columns

0
0 274.367106
1 35.553065
2 522.821937
3 92.515519
4 33.446217
... ...
48124 25.919379
48125 75.859785
48126 48.402383
48127 293.509456
48128 421.146119

48129 rows × 1 columns

0
0 45.795303
1 15.878769
2 21.845367
3 14.488416
4 344.731653
... ...
21401 59.818207
21402 13.964827
21403 128.051507
21404 176.612802
21405 245.816445

21406 rows × 1 columns

0
0 118.713985
1 67.395399
2 27.073464
3 112.184943
4 67.211745
... ...
23009 20.237335
23010 1213.036052
23011 53.378732
23012 99.544672
23013 97.962964

23014 rows × 1 columns

0
0 148.878988
1 88.325662
2 53.639603
3 58.828367
4 1037.404460
... ...
41656 507.615053
41657 317.901180
41658 30.504177
41659 62.302882
41660 728.571526

41661 rows × 1 columns

0
0 110.002848
1 100.119518
2 35.670939
3 104.887662
4 73.708474
... ...
58943 135.970979
58944 99.755129
58945 155.988587
58946 80.796297
58947 116.174652

58948 rows × 1 columns

0
0 18.156982
1 770.564549
2 120.300152
3 62.457731
4 65.750478
... ...
20073 84.231031
20074 914.184940
20075 60.042079
20076 412.368861
20077 19.055842

20078 rows × 1 columns

0
0 398.611511
1 55.839262
2 77.243593
3 1304.790980
4 72.391068
... ...
5174 36.506975
5175 1115.010985
5176 143.689987
5177 48.646845
5178 464.702500

5179 rows × 1 columns

0
0 60.668913
1 34.574231
2 135.356656
3 1201.104334
4 60.986209
... ...
3399 1428.391920
3400 1399.261054
3401 66.499178
3402 98.375089
3403 73.720043

3404 rows × 1 columns

0
0 103.818553
1 113.029484
2 105.535405
3 956.658340
4 37.009373
... ...
5429 109.659727
5430 245.627433
5431 89.075206
5432 46.218475
5433 68.453003

5434 rows × 1 columns

0
0 117.917716
1 60.368189
2 106.026393
3 18.873867
4 867.596778
... ...
11021 49.616098
11022 1108.300986
11023 121.330039
11024 625.404605
11025 103.814676

11026 rows × 1 columns

0
0 1110.502242
1 345.911006
2 39.504982
3 99.835906
4 75.525220
... ...
18500 529.065292
18501 119.301695
18502 106.702676
18503 276.974839
18504 538.951930

18505 rows × 1 columns

0
0 68.286312
1 39.953749
2 60.416889
3 50.038970
4 64.739444
... ...
19499 77.684112
19500 129.915471
19501 106.508003
19502 77.068989
19503 35.953535

19504 rows × 1 columns

0
0 339.181092
1 122.302602
2 634.641964
3 151.480484
4 95.661702
... ...
13481 720.673161
13482 235.111937
13483 729.677239
13484 48.733459
13485 92.768180

13486 rows × 1 columns

0
0 171.258809
1 540.359474
2 73.843787
3 70.820574
4 55.000461
... ...
8420 177.519620
8421 435.200540
8422 1199.026642
8423 1107.307499
8424 118.770827

8425 rows × 1 columns

0
0 83.857029
1 114.315472
2 208.307961
3 188.421035
4 17.385281
... ...
5476 697.596041
5477 94.189151
5478 114.175654
5479 71.607560
5480 60.128946

5481 rows × 1 columns

[7]:
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

[8]:
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 = -5.432E+04 Ncalls = 767 (767 total)
EDM = 0.000149 (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 133 18 -500 1.5E+03
1 r.1.1.0 6 9 -500 1.5E+03
2 r.1.1.1 -1.4 2.2 -500 1.5E+03
3 r.1.2.0 55 14 -500 1.5E+03
4 r.1.2.1 67 11 -500 1.5E+03
5 r.1.2.2 66 9 -500 1.5E+03
6 i.1.0.0 79 30 -500 1.5E+03
7 i.1.1.0 -7 12 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 57 17 -500 1.5E+03
10 i.1.2.1 36 19 -500 1.5E+03
11 i.1.2.2 30 20 -500 1.5E+03
'Bin #1:'
FCN = -8.009E+04 Ncalls = 1359 (1359 total)
EDM = 0.00172 (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 False False True
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 -50 60 -500 1.5E+03
1 r.1.1.0 -39 22 -500 1.5E+03
2 r.1.1.1 2.5 6.5 -500 1.5E+03
3 r.1.2.0 -4 33 -500 1.5E+03
4 r.1.2.1 -33 27 -500 1.5E+03
5 r.1.2.2 -23 28 -500 1.5E+03
6 i.1.0.0 186 18 -500 1.5E+03
7 i.1.1.0 -7 17 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 82 5 -500 1.5E+03
10 i.1.2.1 76 11 -500 1.5E+03
11 i.1.2.2 81 8 -500 1.5E+03
'Bin #2:'
FCN = -1.504E+05 Ncalls = 2607 (2607 total)
EDM = 0.00186 (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 100 160 -500 1.5E+03
1 r.1.1.0 -6 10 -500 1.5E+03
2 r.1.1.1 -4 7 -500 1.5E+03
3 r.1.2.0 50 50 -500 1.5E+03
4 r.1.2.1 70 40 -500 1.5E+03
5 r.1.2.2 60 50 -500 1.5E+03
6 i.1.0.0 260 60 -500 1.5E+03
7 i.1.1.0 2.7 2.9 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 78 28 -500 1.5E+03
10 i.1.2.1 60 40 -500 1.5E+03
11 i.1.2.2 70 40 -500 1.5E+03
'Bin #3:'
FCN = -4.538E+05 Ncalls = 1108 (1108 total)
EDM = 0.00788 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
False True True 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 395 28 -500 1.5E+03
1 r.1.1.0 -8 9 -500 1.5E+03
2 r.1.1.1 -5 4 -500 1.5E+03
3 r.1.2.0 45 11 -500 1.5E+03
4 r.1.2.1 4 14 -500 1.5E+03
5 r.1.2.2 27 12 -500 1.5E+03
6 i.1.0.0 310 40 -500 1.5E+03
7 i.1.1.0 8 11 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 82 11 -500 1.5E+03
10 i.1.2.1 125 10 -500 1.5E+03
11 i.1.2.2 110 11 -500 1.5E+03
'Bin #4:'
FCN = -1.371E+06 Ncalls = 2001 (2001 total)
EDM = 0.00707 (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 580 60 -500 1.5E+03
1 r.1.1.0 10 7 -500 1.5E+03
2 r.1.1.1 5 5 -500 1.5E+03
3 r.1.2.0 -77 13 -500 1.5E+03
4 r.1.2.1 -82 11 -500 1.5E+03
5 r.1.2.2 -76 10 -500 1.5E+03
6 i.1.0.0 640 50 -500 1.5E+03
7 i.1.1.0 -9 4 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 96 6 -500 1.5E+03
10 i.1.2.1 95 9 -500 1.5E+03
11 i.1.2.2 100 11 -500 1.5E+03
'Bin #5:'
FCN = -4.759E+05 Ncalls = 1560 (1560 total)
EDM = 0.000158 (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 16 -500 1.5E+03
1 r.1.1.0 165 11 -500 1.5E+03
2 r.1.1.1 125 13 -500 1.5E+03
3 r.1.2.0 -2.9 5.5 -500 1.5E+03
4 r.1.2.1 -1.8 7.3 -500 1.5E+03
5 r.1.2.2 2.1 11.2 -500 1.5E+03
6 i.1.0.0 475 7 -500 1.5E+03
7 i.1.1.0 6 5 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -106 7 -500 1.5E+03
10 i.1.2.1 -118 7 -500 1.5E+03
11 i.1.2.2 -95.4 2.8 -500 1.5E+03
'Bin #6:'
FCN = -1.975E+05 Ncalls = 1431 (1431 total)
EDM = 6.85E-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 -1.1 14.4 -500 1.5E+03
1 r.1.1.0 94 4 -500 1.5E+03
2 r.1.1.1 74 7 -500 1.5E+03
3 r.1.2.0 6 8 -500 1.5E+03
4 r.1.2.1 -0.4 7.8 -500 1.5E+03
5 r.1.2.2 6 13 -500 1.5E+03
6 i.1.0.0 237.5 3.0 -500 1.5E+03
7 i.1.1.0 3 4 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -150 3 -500 1.5E+03
10 i.1.2.1 -146 4 -500 1.5E+03
11 i.1.2.2 -139.8 2.8 -500 1.5E+03
'Bin #7:'
FCN = -2.144E+05 Ncalls = 1645 (1645 total)
EDM = 2.65E-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 -80 19 -500 1.5E+03
1 r.1.1.0 -5 6 -500 1.5E+03
2 r.1.1.1 -12 7 -500 1.5E+03
3 r.1.2.0 0.7 27.3 -500 1.5E+03
4 r.1.2.1 4 31 -500 1.5E+03
5 r.1.2.2 -8 25 -500 1.5E+03
6 i.1.0.0 -139 11 -500 1.5E+03
7 i.1.1.0 -3.5 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 203.4 2.6 -500 1.5E+03
10 i.1.2.1 198.9 3.0 -500 1.5E+03
11 i.1.2.2 195 4 -500 1.5E+03
'Bin #8:'
FCN = -4.129E+05 Ncalls = 1510 (1510 total)
EDM = 1.59E-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 False False True
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 -3.1 16.3 -500 1.5E+03
1 r.1.1.0 2.1 5.3 -500 1.5E+03
2 r.1.1.1 0.06 3.26 -500 1.5E+03
3 r.1.2.0 156 29 -500 1.5E+03
4 r.1.2.1 243 20 -500 1.5E+03
5 r.1.2.2 167 27 -500 1.5E+03
6 i.1.0.0 105 8 -500 1.5E+03
7 i.1.1.0 1.0 4.9 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -248 19 -500 1.5E+03
10 i.1.2.1 -169 29 -500 1.5E+03
11 i.1.2.2 -223 20 -500 1.5E+03
'Bin #9:'
FCN = -6.071E+05 Ncalls = 1829 (1829 total)
EDM = 0.00232 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
False True True False
Hesse failed Has cov. Accurate Pos. def. Forced
False True False False True
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 50 21 -500 1.5E+03
1 r.1.1.0 -9 9 -500 1.5E+03
2 r.1.1.1 -0.8 10.2 -500 1.5E+03
3 r.1.2.0 -80 80 -500 1.5E+03
4 r.1.2.1 25 89 -500 1.5E+03
5 r.1.2.2 -30 80 -500 1.5E+03
6 i.1.0.0 0.16 10.69 -500 1.5E+03
7 i.1.1.0 6 3 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -344 18 -500 1.5E+03
10 i.1.2.1 -357 7 -500 1.5E+03
11 i.1.2.2 -340 7 -500 1.5E+03
'Bin #10:'
FCN = -1.854E+05 Ncalls = 1098 (1098 total)
EDM = 0.00301 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
False True True False
Hesse failed Has cov. Accurate Pos. def. Forced
False True False False True
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 36 11 -500 1.5E+03
1 r.1.1.0 -8 16 -500 1.5E+03
2 r.1.1.1 1.6 3.4 -500 1.5E+03
3 r.1.2.0 202 12 -500 1.5E+03
4 r.1.2.1 203 6 -500 1.5E+03
5 r.1.2.2 195 4 -500 1.5E+03
6 i.1.0.0 31 15 -500 1.5E+03
7 i.1.1.0 -48 7 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -30 70 -500 1.5E+03
10 i.1.2.1 17 66 -500 1.5E+03
11 i.1.2.2 5 69 -500 1.5E+03
'Bin #11:'
FCN = -4.068E+04 Ncalls = 1721 (1721 total)
EDM = 0.000158 (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 16 59 -500 1.5E+03
1 r.1.1.0 6 9 -500 1.5E+03
2 r.1.1.1 6 8 -500 1.5E+03
3 r.1.2.0 -70 60 -500 1.5E+03
4 r.1.2.1 -70 60 -500 1.5E+03
5 r.1.2.2 -80 60 -500 1.5E+03
6 i.1.0.0 66 19 -500 1.5E+03
7 i.1.1.0 2.8 4.8 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 70 70 -500 1.5E+03
10 i.1.2.1 70 60 -500 1.5E+03
11 i.1.2.2 70 70 -500 1.5E+03
'Bin #12:'
FCN = -2.508E+04 Ncalls = 1372 (1372 total)
EDM = 7.06E-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 -36 15 -500 1.5E+03
1 r.1.1.0 -9 5 -500 1.5E+03
2 r.1.1.1 -5 5 -500 1.5E+03
3 r.1.2.0 60 23 -500 1.5E+03
4 r.1.2.1 57 23 -500 1.5E+03
5 r.1.2.2 51 22 -500 1.5E+03
6 i.1.0.0 36 16 -500 1.5E+03
7 i.1.1.0 1.2 5.0 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 56 24 -500 1.5E+03
10 i.1.2.1 58 22 -500 1.5E+03
11 i.1.2.2 54 21 -500 1.5E+03
'Bin #13:'
FCN = -4.279E+04 Ncalls = 1079 (1079 total)
EDM = 0.000103 (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 10 -500 1.5E+03
1 r.1.1.0 43 7 -500 1.5E+03
2 r.1.1.1 31 9 -500 1.5E+03
3 r.1.2.0 -11 15 -500 1.5E+03
4 r.1.2.1 -12 12 -500 1.5E+03
5 r.1.2.2 -18 16 -500 1.5E+03
6 i.1.0.0 -20 4 -500 1.5E+03
7 i.1.1.0 4 6 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 102.4 3.1 -500 1.5E+03
10 i.1.2.1 102 3 -500 1.5E+03
11 i.1.2.2 97 4 -500 1.5E+03
'Bin #14:'
FCN = -9.5E+04 Ncalls = 1294 (1294 total)
EDM = 0.000535 (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 -11 14 -500 1.5E+03
1 r.1.1.0 2.4 6.6 -500 1.5E+03
2 r.1.1.1 2.8 4.6 -500 1.5E+03
3 r.1.2.0 115 23 -500 1.5E+03
4 r.1.2.1 90 26 -500 1.5E+03
5 r.1.2.2 97 28 -500 1.5E+03
6 i.1.0.0 -4 14 -500 1.5E+03
7 i.1.1.0 -2.9 6.1 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 103 25 -500 1.5E+03
10 i.1.2.1 123 19 -500 1.5E+03
11 i.1.2.2 112 25 -500 1.5E+03
'Bin #15:'
FCN = -1.692E+05 Ncalls = 1111 (1111 total)
EDM = 0.00244 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
False True True False
Hesse failed Has cov. Accurate Pos. def. Forced
False True False False True
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 8 19 -500 1.5E+03
1 r.1.1.0 -2.5 9.3 -500 1.5E+03
2 r.1.1.1 3 4 -500 1.5E+03
3 r.1.2.0 193 9 -500 1.5E+03
4 r.1.2.1 198 3 -500 1.5E+03
5 r.1.2.2 190 15 -500 1.5E+03
6 i.1.0.0 40 40 -500 1.5E+03
7 i.1.1.0 -28 13 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -22 75 -500 1.5E+03
10 i.1.2.1 -2.7 59.9 -500 1.5E+03
11 i.1.2.2 -25 87 -500 1.5E+03
'Bin #16:'
FCN = -1.798E+05 Ncalls = 1159 (1159 total)
EDM = 0.00754 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
False True True False
Hesse failed Has cov. Accurate Pos. def. Forced
False True False False True
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 9 13 -500 1.5E+03
1 r.1.1.0 -2.1 5.4 -500 1.5E+03
2 r.1.1.1 2.3 3.1 -500 1.5E+03
3 r.1.2.0 201 7 -500 1.5E+03
4 r.1.2.1 196 17 -500 1.5E+03
5 r.1.2.2 197 7 -500 1.5E+03
6 i.1.0.0 40 23 -500 1.5E+03
7 i.1.1.0 -14 13 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 17 77 -500 1.5E+03
10 i.1.2.1 50 70 -500 1.5E+03
11 i.1.2.2 19 81 -500 1.5E+03
'Bin #17:'
FCN = -1.193E+05 Ncalls = 1486 (1486 total)
EDM = 0.000675 (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 False False True
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 r.1.0.0 21 13 -500 1.5E+03
1 r.1.1.0 0.5 7.1 -500 1.5E+03
2 r.1.1.1 -0.013 3.887 -500 1.5E+03
3 r.1.2.0 165 7 -500 1.5E+03
4 r.1.2.1 164 16 -500 1.5E+03
5 r.1.2.2 166 6 -500 1.5E+03
6 i.1.0.0 28 25 -500 1.5E+03
7 i.1.1.0 15 17 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 11 96 -500 1.5E+03
10 i.1.2.1 40 80 -500 1.5E+03
11 i.1.2.2 13 97 -500 1.5E+03
'Bin #18:'
FCN = -7.069E+04 Ncalls = 1547 (1547 total)
EDM = 0.000113 (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 12 49 -500 1.5E+03
1 r.1.1.0 -2.7 9.7 -500 1.5E+03
2 r.1.1.1 0.5 4.2 -500 1.5E+03
3 r.1.2.0 110 120 -500 1.5E+03
4 r.1.2.1 90 150 -500 1.5E+03
5 r.1.2.2 90 140 -500 1.5E+03
6 i.1.0.0 -20 41 -500 1.5E+03
7 i.1.1.0 -1.8 13.4 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 -80 170 -500 1.5E+03
10 i.1.2.1 -90 150 -500 1.5E+03
11 i.1.2.2 -90 140 -500 1.5E+03
'Bin #19:'
FCN = -4.365E+04 Ncalls = 1778 (1778 total)
EDM = 0.000126 (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 4 23 -500 1.5E+03
1 r.1.1.0 5 14 -500 1.5E+03
2 r.1.1.1 -10 13 -500 1.5E+03
3 r.1.2.0 -24 29 -500 1.5E+03
4 r.1.2.1 -5 31 -500 1.5E+03
5 r.1.2.2 -22 28 -500 1.5E+03
6 i.1.0.0 26 5 -500 1.5E+03
7 i.1.1.0 3 3 -500 1.5E+03
8 i.1.1.1 0.10 0.10 -500 1.5E+03 yes
9 i.1.2.0 103 7 -500 1.5E+03
10 i.1.2.1 106.7 2.7 -500 1.5E+03
11 i.1.2.2 103 6 -500 1.5E+03
'0.0% (0)Failed. '

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

[9]:
waves = amplitude.make_elm(result.values)
#waves = AmplitudeJPACfit.FitAmplitude.make_elm(result.values)
waves
[9]:
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)

[10]:
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

[11]:
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

[12]:
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)

[13]:
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()
../_images/examples_demo_JPAC_fit_24_0.png

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

[14]:
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)

[15]:
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]))

Plot expected number of events vs mass for each wave

[16]:
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()
../_images/examples_demo_JPAC_fit_30_0.png
../_images/examples_demo_JPAC_fit_30_1.png
../_images/examples_demo_JPAC_fit_30_2.png
../_images/examples_demo_JPAC_fit_30_3.png
../_images/examples_demo_JPAC_fit_30_4.png
../_images/examples_demo_JPAC_fit_30_5.png

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

[17]:
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()
../_images/examples_demo_JPAC_fit_32_0.png

Calculate Moments

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

Calculate PhaseMotion between two waves

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

[21]:
plt.xlim(0.6, 2.)
plt.scatter(mni["mass"],H000,LABEL="H000")
plt.scatter(mni["mass"],H100,LABEL="H100")
plt.legend(loc='upper right')
plt.show()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-21-d1e42b4db3a2> in <module>
      1 plt.xlim(0.6, 2.)
----> 2 plt.scatter(mni["mass"],H000,LABEL="H000")
      3 plt.scatter(mni["mass"],H100,LABEL="H100")
      4 plt.legend(loc='upper right')
      5 plt.show()

~/.conda/envs/pypwa-development/lib/python3.7/site-packages/matplotlib/pyplot.py in scatter(x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, edgecolors, plotnonfinite, data, **kwargs)
   2846         verts=verts, edgecolors=edgecolors,
   2847         plotnonfinite=plotnonfinite, **({"data": data} if data is not
-> 2848         None else {}), **kwargs)
   2849     sci(__ret)
   2850     return __ret

~/.conda/envs/pypwa-development/lib/python3.7/site-packages/matplotlib/__init__.py in inner(ax, data, *args, **kwargs)
   1597     def inner(ax, *args, data=None, **kwargs):
   1598         if data is None:
-> 1599             return func(ax, *map(sanitize_sequence, args), **kwargs)
   1600
   1601         bound = new_sig.bind(ax, *args, **kwargs)

~/.conda/envs/pypwa-development/lib/python3.7/site-packages/matplotlib/axes/_axes.py in scatter(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, edgecolors, plotnonfinite, **kwargs)
   4441         y = np.ma.ravel(y)
   4442         if x.size != y.size:
-> 4443             raise ValueError("x and y must be the same size")
   4444
   4445         if s is None:

ValueError: x and y must be the same size
../_images/examples_demo_JPAC_fit_36_1.png
[18]:
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

[19]:
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)
[19]:
<ErrorbarContainer object of 3 artists>
../_images/examples_demo_JPAC_fit_39_1.png

Write fitted values (of production amplitudes) to disk

[20]:
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)
[21]:
pandas.DataFrame([res.np_errors() for res in results], columns=results[0].parameters)
[21]:
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 2.852249 2.454338 2.233434 2.357102 2.383402 2.626812 30.355296 23.058165 0.1 29.117763 19.515329 21.461233
1 4.143564 6.957555 2.850610 2.725406 2.577218 3.350895 24.972032 16.099636 0.1 20.407775 15.720090 14.143814
2 5.692306 4.785399 2.383390 8.044186 5.157728 3.613297 35.220094 15.107411 0.1 14.095800 13.392294 16.332394
3 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.1 0.000000 0.000000 0.000000
4 44.106133 13.751319 19.829211 10.964590 13.588421 15.587226 12.307185 5.567356 0.1 10.009993 9.033432 6.326170
5 34.436933 9.159502 11.498057 5.682668 8.133921 7.363222 28.735040 13.816089 0.1 11.345555 17.985841 17.540164
6 7.198252 8.374227 4.763529 4.189970 4.520014 3.353074 17.183687 11.151904 0.1 13.529245 17.287970 16.054413
7 11.671223 5.996570 4.279334 7.625037 5.538663 12.454837 19.612032 6.767201 0.1 15.038266 14.425965 25.855571
8 6.475163 5.382701 5.416978 8.161982 8.444590 8.678492 17.661079 8.321693 0.1 17.488448 10.915935 22.474796
9 8.721956 8.181113 8.761297 13.376127 16.434883 13.947150 17.630661 12.990305 0.1 25.651122 14.814478 25.210855
10 8.493679 6.084142 8.505738 14.310373 10.479354 6.749547 8.246809 5.742596 0.1 10.345708 5.497219 9.727813
11 3.855473 5.012773 4.741625 7.609596 3.837674 5.775665 7.119417 12.499101 0.1 5.919945 5.846466 7.079995
12 3.037836 4.442392 5.428389 3.954342 3.764583 3.667913 13.051691 11.682621 0.1 9.414705 9.861358 12.532426
13 5.321230 7.210871 6.045126 4.710401 3.813843 4.735120 17.598563 16.636521 0.1 15.900800 14.740404 13.281171
14 3.567710 3.193201 3.337286 4.090491 3.627887 4.866442 8.841624 9.924517 0.1 8.551742 7.499289 9.294340
15 6.727505 3.681380 4.365790 3.104282 5.219031 4.103600 12.789609 16.679240 0.1 19.571438 13.145244 15.012328
16 6.177552 2.676033 3.828072 3.833456 3.289779 5.726847 18.829520 13.113563 0.1 20.012938 17.978662 23.420750
17 2.933203 3.324188 3.654907 3.609211 4.131418 2.916676 19.305438 13.300829 0.1 20.013559 11.407195 22.344017
18 3.051973 3.116716 3.546187 5.042575 3.908125 5.411048 18.860401 11.611010 0.1 21.188980 12.235497 23.397019
19 2.687607 3.314947 4.027697 2.654148 4.601623 2.709061 23.221308 25.126377 0.1 25.907302 15.642661 23.220975
[22]:
final_values
[22]:
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 111.704412 60.453623 59.805884 61.078017 61.115139 62.403759 1.320825 5.792675 0.1 2.262332 4.399611 -5.139527
1 137.562715 64.190450 62.573028 65.067900 65.866904 70.420397 -2.994290 29.600163 0.1 6.868300 2.703924 -7.967649
2 198.184937 64.268233 66.615670 64.595816 64.606409 69.656890 30.953681 20.593808 0.1 45.946604 29.400564 10.318831
3 337.491982 72.959054 66.363838 71.371487 88.351834 80.233547 110.288478 -5.319404 0.1 -1.097682 -44.011232 -12.085719
4 191.948115 72.071364 103.066780 91.422853 83.604409 72.393911 585.360143 5.557063 0.1 -13.656104 -10.491047 -6.151311
5 -229.743178 82.435298 69.243004 91.917776 101.022541 107.038100 275.792660 29.126976 0.1 -10.282243 -3.671119 -0.001849
6 -151.579375 86.849678 71.505687 129.835388 127.271735 119.766565 79.726460 59.105635 0.1 -22.104475 31.613361 10.078527
7 -76.655901 75.693118 103.524504 146.645818 156.565194 152.959609 -57.516859 68.733135 0.1 75.957406 47.150633 61.484341
8 -65.816781 86.599772 124.420103 216.389109 200.722624 222.109393 13.190943 59.468474 0.1 99.732686 131.442304 84.706465
9 -26.800948 67.763534 105.828129 269.445703 216.992426 244.480875 42.960919 39.327445 0.1 130.453626 215.905357 153.900684
10 -67.342887 58.718273 74.548323 111.526792 83.759368 129.187741 -127.001798 -1.308107 0.1 -132.286247 -157.557282 -123.256984
11 -13.410233 144.691054 144.451525 39.366624 43.680159 36.585875 37.882643 -3.078285 0.1 76.413219 73.946001 83.273424
12 -5.720551 162.913340 153.937952 58.902362 60.435890 57.743208 -28.108761 11.247377 0.1 -20.433062 -25.311471 -50.581979
13 1.810090 165.935186 165.998990 78.621711 77.283796 79.933338 -63.243685 0.625528 0.1 4.086042 -17.541890 17.548998
14 1.808067 165.754963 161.907458 105.701858 106.483733 107.956861 -24.282072 -11.740948 0.1 50.413737 43.696063 41.968264
15 4.371681 144.679254 145.952516 149.573552 143.805461 152.500957 71.885659 14.356363 0.1 13.678425 53.303994 -21.211547
16 14.178999 130.629862 126.308624 157.756320 155.465998 155.974918 39.733446 -0.988772 0.1 -25.118227 -10.691878 -22.150941
17 15.816585 114.488989 117.105987 129.755851 126.990058 125.500639 23.952670 6.442825 0.1 20.606003 31.539561 6.359807
18 12.896210 97.348436 99.699072 101.426950 103.016910 102.753527 11.472671 6.887124 0.1 19.385896 17.023952 23.056549
19 17.379010 80.181185 83.618056 83.784640 88.080221 85.754140 -10.540275 -8.638290 0.1 5.881675 -16.019906 1.178200
[ ]: