tools provided by tutors

Gammaspektroskopie/tools/Data/compton_example_data.hst

+# simulated compton spectrum
+#  counts per channel
+2243
+2013
+2041
+2137
+1847
+1924
+2106
+2076
+2035
+1969
+2099
+1999
+2058
+1916
+1839
+2009
+2072
+1869
+1990
+1828
+1755
+1867
+1891
+1913
+1768
+1981
+1872
+1874
+1708
+1695
+1788
+1812
+1900
+1871
+1881
+1933
+1845
+1708
+1921
+1641
+1745
+1712
+1653
+1481
+1708
+1764
+1522
+1780
+1725
+1771
+1594
+1676
+1768
+1703
+1739
+1644
+1687
+1649
+1663
+1710
+1624
+1610
+1600
+1436
+1743
+1482
+1365
+1676
+1685
+1564
+1545
+1611
+1485
+1460
+1466
+1720
+1538
+1598
+1524
+1596
+1487
+1546
+1463
+1617
+1493
+1517
+1537
+1527
+1587
+1481
+1520
+1481
+1516
+1340
+1584
+1415
+1424
+1403
+1492
+1622
+1529
+1469
+1598
+1339
+1428
+1561
+1359
+1561
+1488
+1456
+1336
+1503
+1467
+1531
+1309
+1321
+1318
+1353
+1284
+1502
+1292
+1385
+1311
+1379
+1403
+1332
+1331
+1419
+1399
+1426
+1320
+1371
+1242
+1347
+1452
+1450
+1337
+1419
+1363
+1423
+1256
+1286
+1210
+1349
+1371
+1546
+1376
+1477
+1381
+1440
+1351
+1322
+1318
+1403
+1266
+1432
+1436
+1403
+1289
+1386
+1323
+1366
+1239
+1391
+1392
+1439
+1492
+1315
+1404
+1462
+1523
+1491
+1266
+1366
+1399
+1375
+1334
+1257
+1381
+1381
+1307
+1449
+1388
+1512
+1475
+1441
+1459
+1343
+1458
+1328
+1332
+1466
+1445
+1405
+1319
+1321
+1285
+1466
+1614
+1423
+1466
+1388
+1393
+1348
+1453
+1443
+1288
+1206
+1423
+1473
+1419
+1593
+1419
+1343
+1355
+1507
+1639
+1489
+1438
+1555
+1380
+1503
+1413
+1409
+1478
+1442
+1394
+1592
+1391
+1293
+1330
+1372
+1640
+1542
+1513
+1514
+1533
+1561
+1559
+1287
+1525
+1519
+1519
+1683
+1540
+1580
+1567
+1412
+1391
+1535
+1540
+1517
+1338
+1426
+1499
+1574
+1577
+1454
+1564
+1557
+1335
+1575
+1292
+1420
+1797
+1498
+1529
+1609
+1484
+1589
+1432
+1642
+1473
+1615
+1633
+1507
+1506
+1608
+1551
+1471
+1567
+1644
+1716
+1499
+1437
+1629
+1507
+1433
+1602
+1567
+1425
+1489
+1456
+1650
+1428
+1503
+1530
+1428
+1430
+1435
+1555
+1434
+1289
+1358
+1437
+1510
+1466
+1395
+1441
+1344
+1204
+1273
+1277
+1269
+1194
+1316
+1098
+1123
+1150
+1218
+1197
+1051
+1216
+1117
+1152
+1170
+944
+931
+1071
+950
+1011
+1044
+992
+966
+916
+954
+763
+731
+792
+876
+915
+750
+721
+837
+742
+659
+865
+708
+615
+627
+667
+613
+592
+547
+588
+553
+609
+547
+475
+525
+451
+547
+558
+387
+466
+444
+532
+442
+422
+291
+375
+409
+373
+331
+389
+398
+351
+365
+284
+292
+276
+340
+280
+309
+337
+306
+191
+290
+225
+252
+275
+330
+173
+253
+218
+276
+227
+291
+224
+153
+248
+254
+207
+265
+242
+193
+276
+125
+179
+260
+209
+199
+202
+187
+110
+224
+142
+108
+152
+195
+240
+198
+213
+204
+193
+100
+191
+158
+327
+210
+237
+188
+184
+94
+254
+214
+170
+176
+150
+96
+141
+222
+223
+179
+184
+231
+188
+158
+186
+174
+135
+181
+185
+206
+77
+254
+242
+234
+248
+319
+188
+216
+190
+194
+179
+167
+284
+261
+145
+70
+184
+136
+130
+161
+221
+200
+207
+182
+91
+169
+249
+145
+185
+152
+239
+216
+299
+214
+131
+228
+124
+110
+100
+160
+248
+291
+163
+139
+185
+112

Gammaspektroskopie/tools/Data/peak_example_data.hst

+# simulated peak spectrum
+#  counts per channel
+219
+200
+85
+196
+165
+185
+87
+272
+117
+67
+125
+175
+216
+186
+204
+156
+88
+183
+101
+68
+176
+85
+167
+84
+129
+160
+147
+171
+182
+183
+205
+156
+167
+175
+125
+152
+263
+133
+224
+226
+144
+187
+220
+148
+239
+231
+269
+194
+200
+204
+184
+269
+231
+261
+146
+243
+194
+262
+228
+247
+263
+236
+254
+326
+321
+255
+334
+273
+236
+299
+305
+246
+311
+385
+299
+303
+334
+283
+340
+307
+445
+396
+307
+307
+447
+368
+420
+386
+498
+403
+442
+404
+491
+402
+449
+508
+522
+473
+520
+513
+551
+562
+530
+576
+539
+722
+646
+521
+683
+674
+679
+707
+691
+734
+746
+776
+760
+843
+918
+899
+843
+710
+856
+1079
+907
+975
+1129
+1170
+1053
+965
+1034
+1156
+1020
+1049
+1044
+1039
+1334
+1188
+1189
+1153
+1308
+1282
+1388
+1411
+1529
+1348
+1423
+1549
+1615
+1648
+1541
+1457
+1557
+1569
+1671
+1581
+1690
+1773
+1672
+1755
+1641
+1926
+1753
+1735
+1861
+1919
+1970
+1734
+2055
+1859
+1964
+1941
+2004
+2106
+2231
+2027
+2111
+2158
+2208
+2187
+2373
+2257
+2198
+2093
+2327
+2153
+2320
+2391
+2349
+2463
+2239
+2460
+2244
+2395
+2407
+2340
+2353
+2468
+2534
+2432
+2454
+2469
+2449
+2387
+2171
+2272
+2442
+2084
+2345
+2358
+2417
+2410
+2337
+2310
+2301
+2236
+2364
+2320
+2285
+2360
+2463
+2348
+2122
+2523
+2194
+2361
+2179
+2183
+2090
+2160
+2222
+2336
+2253
+2119
+2178
+2236
+2165
+2139
+2139
+2318
+2166
+2046
+2199
+2047
+1963
+1998
+2044
+1860
+2060
+1939
+1844
+1717
+1688
+1747
+1817
+1835
+1995
+1810
+1762
+1579
+1668
+1669
+1588
+1528
+1447
+1485
+1635
+1396
+1411
+1376
+1478
+1396
+1471
+1296
+1209
+1394
+1236
+1302
+1215
+1159
+1021
+987
+1097
+1066
+980
+1132
+963
+1120
+988
+977
+895
+900
+918
+833
+994
+776
+737
+700
+720
+712
+811
+675
+715
+691
+743
+653
+681
+665
+668
+655
+640
+703
+580
+512
+509
+506
+597
+541
+371
+552
+429
+449
+415
+420
+380
+488
+418
+370
+398
+385
+415
+421
+396
+239
+364
+297
+371
+285
+362
+297
+259
+290
+316
+234
+276
+207
+314
+209
+272
+96
+279
+245
+239
+190
+286
+329
+215
+199
+215
+101
+265
+240
+164
+254
+221
+185
+199
+193
+168
+217
+180
+171
+204
+227
+167
+148
+158
+168
+177
+225
+234
+98
+249
+140
+210
+293
+153
+126
+197
+144
+263
+205
+140
+171
+159
+50
+144
+180
+218
+95
+219
+79
+104
+170
+189
+135
+207
+113
+210
+176
+224
+135
+32
+202
+126
+122
+166
+107
+147
+161
+180
+147
+120
+146
+171
+172
+220
+79
+186
+128
+81
+158
+58
+41
+101
+170
+151
+185
+137
+151
+180
+143
+240
+152
+122
+226
+114
+150
+184
+135
+185
+218
+175
+171
+106
+123
+187
+136
+171
+213
+111
+156
+134
+122
+178
+119
+181
+231
+126
+147
+102
+228
+190
+126
+84
+147
+191
+276
+165
+126
+260
+170
+123
+216
+189
+72
+169
+149
+202
+191
+152
+163
+123
+136
+114
+145
+185
+166
+96
+154

Gammaspektroskopie/tools/gamma_tools.py

+import numpy as np
+from matplotlib import pyplot as plt
+from kafe2 import HistContainer, HistFit, Plot
+from scipy import special
+
+
+
+########################################################################################################################################
+# TABLE OF CONTENTS #
+#####################
+#
+# - Gaussian model function for photopeaks
+# - Model function for compton edges using the complementary error function 'erfc' (Convolution of Gaussian and Heaviside-Step-Function) with a quadratic form correction factor 
+# - Fit a photopeak using a Gaussian model function
+# - Fit a compton edge using a complementary error function
+# - Custom wrapper for fitting a custom function to a histogram using Kafe2, combining all relevant steps
+#
+########################################################################################################################################
+
+
+
+##########################################
+# Gaussian model function for photopeaks #
+##########################################
+
+def gaussian (x, mu, sigma = 10, baseline = 0, A = 100):
+    return  A * np.exp(-(x - mu)**2 / (2 * (sigma**2))) + baseline
+
+
+
+################################################################################################################################################################################
+# Model function for compton edges using the complementary error function 'erfc' (Convolution of Gaussian and Heaviside-Step-Function) with a quadratic form correction factor #
+################################################################################################################################################################################
+
+def compton_edge (x, mu, sigma = 10, baseline = 0, A = 100, B = 1):
+    return A * special.erfc((x-mu)/(np.sqrt(2) * sigma)) * (1 + B * (x/mu - 0.5)**2) + baseline
+
+
+
+###################################################
+# Fit a photopeak using a Gaussian model function #
+###################################################
+
+def fit_peak (datax, datay, x1 = None, x2 = None, baseline = 0, showResult = True, label = None, **initparams):
+    # Check if a lower x-limit for the fit is given, otherwise use the minimum of the provided dataset
+    if x1 == None:
+        x1 = min(datax)
+    # Check if an upper x-limit for the fit is given, otherwise use the maximum of the provided dataset
+    if x2 == None:
+        x2 = max(datax)
+    # Remove datapoints outside the specified fit range
+    datax_clipped = np.delete(datax, np.where((x1 > datax) | (datax > x2)))
+    datay_clipped = np.delete(datay, np.where((x1 > datax) | (datax > x2)))
+    # Set axis labels
+    axislabels = ['Channel', 'Count']
+    # Set model LaTeX expression
+    modelexpression = 'A \\cdot e^{{\\frac{{-(x-\\mu)^2}}{{2 \\sigma^2}}}} + \\mathrm{{baseline}}'
+    # Perform a Gaussian fit on the remaining data, choosing the position of the maximum value in the data as an initial value for the position of the peak and it's value for the height of the peak
+    fit = hist_fit_data(datax_clipped, datay_clipped, model = gaussian, limitedparams = ['mu', 'sigma'], fixedparams = [["baseline", baseline]], mu = datax_clipped[np.argmax(datay_clipped)], A = np.max(datay_clipped), label = label, showResult = showResult, axislabels = axislabels, modelexpression = modelexpression, **initparams)
+    # Return fit object
+    return fit
+
+
+
+###########################################################
+# Fit a compton edge using a complementary error function #
+# ##########################################################
+
+def fit_compton (datax, datay, x1 = None, x2 = None, baseline = 0, showResult = True, label = None, **initparams):
+    # Check if a lower x-limit for the fit is given, otherwise use the minimum of the provided dataset
+    if x1 == None:
+        x1 = min(datax)
+    # Check if an upper x-limit for the fit is given, otherwise use the maximum of the provided dataset
+    if x2 == None:
+        x2 = max(datax)
+    # Remove datapoints outside the specified fit range
+    datax_clipped = np.delete(datax, np.where((x1 > datax) | (datax > x2)))
+    datay_clipped = np.delete(datay, np.where((x1 > datax) | (datax > x2)))
+    # Determine an initial value for the position of the compton edge by estimating where the slope is maximal
+    dif = np.diff(datay_clipped)
+    diffs = []
+    for i in range(len(datay_clipped-25)):
+        diffs.append(np.abs(np.sum(dif[i:i+25])))
+    mu0 = datax_clipped[np.argmax(diffs)]
+    # Set axis labels
+    axislabels = ['Channel', 'Count']
+    # Set model LaTeX expression
+    modelexpression = 'A \\cdot \\mathrm{{erfc}}(\\frac{{x-mu}}{{\\sqrt{{2}} \\sigma}}) \\cdot (1 + B (\\frac{{x}}{{mu}} - \\frac{{1}}{{2}})^2) + \\mathrm{{baseline}}'
+    # Perform the fit on the remaining data using the estimated position as an initial value
+    fit = hist_fit_data(datax_clipped, datay_clipped, model = compton_edge, limitedparams = ['mu', 'sigma', 'B'], fixedparams = [["baseline", baseline]], mu = mu0, A = np.max(datay_clipped), label = label, axislabels = axislabels, modelexpression = modelexpression, showResult = showResult, **initparams)
+    # Return fit object
+    return fit
+
+
+
+#########################################################################################################
+# Custom wrapper for fitting a custom function to a histogram using Kafe2, combining all relevant steps #
+# ########################################################################################################
+
+def hist_fit_data (x, y, model = None, label = 'data', constraints = [], fixedparams = [], limitedparams = [], parameternames = None, modellabel = None, modelname = None, axislabels = [None, None], modelexpression = None, report = False, showResult = False, **initialvalues):
+    # Organise data histogram container
+    step = x[1]-x[0]
+    data = HistContainer(bin_edges = np.arange(x[0]-0.5*step, x[-1]+step, step = step))
+    data.set_bins(y)
+    data.label = label  
+    data.axis_labels = axislabels
+    # Assign data and model to fit
+    fit = HistFit(data, model_function = model, density = False)  
+    # Assign parameter constraints and set as initial values
+    for a in constraints:
+        fit.add_parameter_constraint(name = a[0], value = a[1], uncertainty = a[2])
+        fit.set_parameter_values(**{a[0]:a[1]})  
+    # Assign fixed parameter values (Causes errors concerning non-positive uncertainties)
+    for b in fixedparams:
+        fit.fix_parameter(name = b[0], value = b[1])   
+    # Assign additional initial values
+    fit.set_parameter_values(**initialvalues)
+    # Limit parameters to be either zero or positive, not negative
+    for a in limitedparams:
+        fit.limit_parameter(a, lower = 0)
+    #Assign parameter names as well as model label, name and expression
+    if parameternames != None:
+        fit.assign_parameter_latex_names(**parameternames)
+    if modellabel != None:
+        fit.model_label = modellabel
+    if modelname != None:
+        fit.assign_model_function_latex_name(modelname)
+    if modelexpression != None:
+        fit.assign_model_function_latex_expression(modelexpression)
+    # Perform fit
+    fit.do_fit()  
+    # Plot fit results if desired
+    if showResult:
+        p = Plot(fit)
+        p.plot()
+        plt.show()
+    # Report fit results if desired
+    if report:
+        fit.report()
+    # Return fit
+    return fit
\ No newline at end of file