gwm17/spspy
1
0
Fork 0
mirror of https://github.com/gwm17/spspy.git synced 2024-11-26 12:08:51 -05:00
Code Issues Projects Releases Wiki Activity

Merge pull request #2 from gwm17/devel

Browse Source 
Make the code actual python code not just C++ pretending to be python
This commit is contained in:
Gordon McCann 2022-12-22 15:00:33 -05:00 committed by GitHub
parent 1a7a139ef4 2562156943
commit 3a3c6aa28d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
43 changed files with 1960 additions and 3945 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
.gitignore vendored
View File

@ -1,8 +1,8 @@
*.pickle .venv/
*.txt .env/
!test_config.txt __pycache__/
__pycache__ .vscode/
*.sublime-project .vs/
*.sublime-workspace *.sps
!/etc/mass.txt *.spanc
!.gitignore *.csv

887
11C_fit.eps
View File

@ -1,887 +0,0 @@
%!PS-Adobe-3.0 EPSF-3.0
%%Title: 11C_fit.eps
%%Creator: Matplotlib v3.5.1, https://matplotlib.org/
%%CreationDate: Wed Mar 2 17:08:47 2022
%%Orientation: portrait
%%BoundingBox: 90 320 522 472
%%HiResBoundingBox: 90.000000 320.400000 522.000000 471.600000
%%EndComments
%%BeginProlog
/mpldict 12 dict def
mpldict begin
/_d { bind def } bind def
/m { moveto } _d
/l { lineto } _d
/r { rlineto } _d
/c { curveto } _d
/cl { closepath } _d
/ce { closepath eofill } _d
/box {
m
1 index 0 r
0 exch r
neg 0 r
cl
} _d
/clipbox {
box
clip
newpath
} _d
/sc { setcachedevice } _d
%!PS-Adobe-3.0 Resource-Font
%%Creator: Converted from TrueType to Type 3 by Matplotlib.
10 dict begin
/FontName /DejaVuSans def
/PaintType 0 def
/FontMatrix [0.00048828125 0 0 0.00048828125 0 0] def
/FontBBox [-2090 -948 3673 2524] def
/FontType 3 def
/Encoding [/space /c /parenleft /parenright /m /zero /one /minus /two /four /six /seven /eight] def
/CharStrings 14 dict dup begin
/.notdef 0 def
/space{651 0 0 0 0 0 sc
ce} _d
/c{1126 0 113 -29 999 1147 sc
999 1077 m
999 905 l
947 934 895 955 842 969 c
790 984 737 991 684 991 c
565 991 472 953 406 877 c
340 802 307 696 307 559 c
307 422 340 316 406 240 c
472 165 565 127 684 127 c
737 127 790 134 842 148 c
895 163 947 184 999 213 c
999 43 l
948 19 894 1 839 -11 c
784 -23 726 -29 664 -29 c
495 -29 361 24 262 130 c
163 236 113 379 113 559 c
113 742 163 885 263 990 c
364 1095 501 1147 676 1147 c
733 1147 788 1141 842 1129 c
896 1118 948 1100 999 1077 c
ce} _d
/parenleft{799 0 176 -270 635 1554 sc
635 1554 m
546 1401 479 1249 436 1099 c
393 949 371 797 371 643 c
371 489 393 336 436 185 c
480 34 546 -117 635 -270 c
475 -270 l
375 -113 300 41 250 192 c
201 343 176 494 176 643 c
176 792 201 941 250 1092 c
299 1243 374 1397 475 1554 c
635 1554 l
ce} _d
/parenright{799 0 164 -270 623 1554 sc
164 1554 m
324 1554 l
424 1397 499 1243 548 1092 c
598 941 623 792 623 643 c
623 494 598 343 548 192 c
499 41 424 -113 324 -270 c
164 -270 l
253 -117 319 34 362 185 c
406 336 428 489 428 643 c
428 797 406 949 362 1099 c
319 1249 253 1401 164 1554 c
ce} _d
/m{1995 0 186 0 1821 1147 sc
1065 905 m
1111 988 1166 1049 1230 1088 c
1294 1127 1369 1147 1456 1147 c
1573 1147 1663 1106 1726 1024 c
1789 943 1821 827 1821 676 c
1821 0 l
1636 0 l
1636 670 l
1636 777 1617 857 1579 909 c
1541 961 1483 987 1405 987 c
1310 987 1234 955 1179 892 c
1124 829 1096 742 1096 633 c
1096 0 l
911 0 l
911 670 l
911 778 892 858 854 909 c
816 961 757 987 678 987 c
584 987 509 955 454 891 c
399 828 371 742 371 633 c
371 0 l
186 0 l
186 1120 l
371 1120 l
371 946 l
413 1015 463 1065 522 1098 c
581 1131 650 1147 731 1147 c
812 1147 881 1126 938 1085 c
995 1044 1038 984 1065 905 c
ce} _d
/zero{1303 0 135 -29 1167 1520 sc
651 1360 m
547 1360 469 1309 416 1206 c
364 1104 338 950 338 745 c
338 540 364 387 416 284 c
469 182 547 131 651 131 c
756 131 834 182 886 284 c
939 387 965 540 965 745 c
965 950 939 1104 886 1206 c
834 1309 756 1360 651 1360 c
651 1520 m
818 1520 946 1454 1034 1321 c
1123 1189 1167 997 1167 745 c
1167 494 1123 302 1034 169 c
946 37 818 -29 651 -29 c
484 -29 356 37 267 169 c
179 302 135 494 135 745 c
135 997 179 1189 267 1321 c
356 1454 484 1520 651 1520 c
ce} _d
/one{1303 0 225 0 1114 1493 sc
254 170 m
584 170 l
584 1309 l
225 1237 l
225 1421 l
582 1493 l
784 1493 l
784 170 l
1114 170 l
1114 0 l
254 0 l
254 170 l
ce} _d
/minus{1716 0 217 557 1499 727 sc
217 727 m
1499 727 l
1499 557 l
217 557 l
217 727 l
ce} _d
/two{1303 0 150 0 1098 1520 sc
393 170 m
1098 170 l
1098 0 l
150 0 l
150 170 l
227 249 331 356 463 489 c
596 623 679 709 713 748 c
778 821 823 882 848 932 c
874 983 887 1032 887 1081 c
887 1160 859 1225 803 1275 c
748 1325 675 1350 586 1350 c
523 1350 456 1339 385 1317 c
315 1295 240 1262 160 1217 c
160 1421 l
241 1454 317 1478 388 1495 c
459 1512 523 1520 582 1520 c
737 1520 860 1481 952 1404 c
1044 1327 1090 1223 1090 1094 c
1090 1033 1078 974 1055 919 c
1032 864 991 800 930 725 c
913 706 860 650 771 557 c
682 465 556 336 393 170 c
ce} _d
/four{1303 0 100 0 1188 1493 sc
774 1317 m
264 520 l
774 520 l
774 1317 l
721 1493 m
975 1493 l
975 520 l
1188 520 l
1188 352 l
975 352 l
975 0 l
774 0 l
774 352 l
100 352 l
100 547 l
721 1493 l
ce} _d
/six{1303 0 143 -29 1174 1520 sc
676 827 m
585 827 513 796 460 734 c
407 672 381 587 381 479 c
381 372 407 287 460 224 c
513 162 585 131 676 131 c
767 131 838 162 891 224 c
944 287 971 372 971 479 c
971 587 944 672 891 734 c
838 796 767 827 676 827 c
1077 1460 m
1077 1276 l
1026 1300 975 1318 923 1331 c
872 1344 821 1350 770 1350 c
637 1350 535 1305 464 1215 c
394 1125 354 989 344 807 c
383 865 433 909 492 940 c
551 971 617 987 688 987 c
838 987 956 941 1043 850 c
1130 759 1174 636 1174 479 c
1174 326 1129 203 1038 110 c
947 17 827 -29 676 -29 c
503 -29 371 37 280 169 c
189 302 143 494 143 745 c
143 981 199 1169 311 1309 c
423 1450 573 1520 762 1520 c
813 1520 864 1515 915 1505 c
967 1495 1021 1480 1077 1460 c
ce} _d
/seven{1303 0 168 0 1128 1493 sc
168 1493 m
1128 1493 l
1128 1407 l
586 0 l
375 0 l
885 1323 l
168 1323 l
168 1493 l
ce} _d
/eight{1303 0 139 -29 1163 1520 sc
651 709 m
555 709 479 683 424 632 c
369 581 342 510 342 420 c
342 330 369 259 424 208 c
479 157 555 131 651 131 c
747 131 823 157 878 208 c
933 260 961 331 961 420 c
961 510 933 581 878 632 c
823 683 748 709 651 709 c
449 795 m
362 816 295 857 246 916 c
198 975 174 1048 174 1133 c
174 1252 216 1347 301 1416 c
386 1485 503 1520 651 1520 c
800 1520 916 1485 1001 1416 c
1086 1347 1128 1252 1128 1133 c
1128 1048 1104 975 1055 916 c
1007 857 940 816 854 795 c
951 772 1027 728 1081 662 c
1136 596 1163 515 1163 420 c
1163 275 1119 164 1030 87 c
942 10 816 -29 651 -29 c
486 -29 360 10 271 87 c
183 164 139 275 139 420 c
139 515 166 596 221 662 c
276 728 352 772 449 795 c
375 1114 m
375 1037 399 976 447 933 c
496 890 564 868 651 868 c
738 868 805 890 854 933 c
903 976 928 1037 928 1114 c
928 1191 903 1252 854 1295 c
805 1338 738 1360 651 1360 c
564 1360 496 1338 447 1295 c
399 1252 375 1191 375 1114 c
ce} _d
end readonly def
/BuildGlyph {
exch begin
CharStrings exch
2 copy known not {pop /.notdef} if
true 3 1 roll get exec
end
} _d
/BuildChar {
1 index /Encoding get exch get
1 index /BuildGlyph get exec
} _d
FontName currentdict end definefont pop
%!PS-Adobe-3.0 Resource-Font
%%Creator: Converted from TrueType to Type 3 by Matplotlib.
10 dict begin
/FontName /DejaVuSans-Oblique def
/PaintType 0 def
/FontMatrix [0.00048828125 0 0 0.00048828125 0 0] def
/FontBBox [-2080 -717 3398 2187] def
/FontType 3 def
/Encoding [/x /rho] def
/CharStrings 3 dict dup begin
/.notdef 0 def
/x{1212 0 -53 0 1229 1120 sc
1229 1120 m
715 571 l
1030 0 l
819 0 l
582 444 l
170 0 l
-53 0 l
498 590 l
205 1120 l
416 1120 l
631 715 l
1006 1120 l
1229 1120 l
ce} _d
/rho{1300 0 33 -426 1278 1147 sc
385 920 m
438 988 521 1052 634 1112 c
678 1135 761 1147 882 1147 c
1018 1147 1118 1093 1182 985 c
1246 877 1261 735 1227 559 c
1192 383 1122 241 1016 133 c
910 25 789 -29 653 -29 c
571 -29 504 -13 451 20 c
398 52 359 101 334 168 c
218 -426 l
33 -426 l
223 549 l
253 703 307 827 385 920 c
1036 559 m
1062 694 1055 801 1014 878 c
973 955 904 993 807 993 c
710 993 626 955 555 878 c
484 801 436 694 410 559 c
383 424 390 317 431 240 c
472 163 541 125 638 125 c
735 125 819 163 890 240 c
961 317 1009 424 1036 559 c
ce} _d
end readonly def
/BuildGlyph {
exch begin
CharStrings exch
2 copy known not {pop /.notdef} if
true 3 1 roll get exec
end
} _d
/BuildChar {
1 index /Encoding get exch get
1 index /BuildGlyph get exec
} _d
FontName currentdict end definefont pop
end
%%EndProlog
mpldict begin
90 320.4 translate
432 151.2 0 0 clipbox
0.500 setlinewidth
0 setlinejoin
0 setlinecap
[] 0 setdash
0.000 setgray
gsave
0 0 m
432 0 l
432 151.2 l
0 151.2 l
cl
gsave
1.000 setgray
fill
grestore
stroke
grestore
gsave
37.71899 33.672115 m
428.99976 33.672115 l
428.99976 148.19976 l
37.71899 148.19976 l
cl
1.000 setgray
fill
grestore
0.800 setlinewidth
1 setlinejoin
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
0 0 m
0 -3.5 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
54.8816 33.6721 o
grestore
/DejaVuSans 10.000 selectfont
gsave
41.155048 19.078365 translate
0.000000 rotate
0.000000 0 m /minus glyphshow
8.378906 0 m /two glyphshow
14.741211 0 m /six glyphshow
21.103516 0 m /zero glyphshow
grestore
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
0 0 m
0 -3.5 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
113.126 33.6721 o
grestore
gsave
99.399278 19.078365 translate
0.000000 rotate
0.000000 0 m /minus glyphshow
8.378906 0 m /two glyphshow
14.741211 0 m /four glyphshow
21.103516 0 m /zero glyphshow
grestore
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
0 0 m
0 -3.5 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
171.37 33.6721 o
grestore
gsave
157.643508 19.078365 translate
0.000000 rotate
0.000000 0 m /minus glyphshow
8.378906 0 m /two glyphshow
14.741211 0 m /two glyphshow
21.103516 0 m /zero glyphshow
grestore
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
0 0 m
0 -3.5 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
229.614 33.6721 o
grestore
gsave
215.887739 19.078365 translate
0.000000 rotate
0.000000 0 m /minus glyphshow
8.378906 0 m /two glyphshow
14.741211 0 m /zero glyphshow
21.103516 0 m /zero glyphshow
grestore
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
0 0 m
0 -3.5 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
287.859 33.6721 o
grestore
gsave
274.131969 19.078365 translate
0.000000 rotate
0.000000 0 m /minus glyphshow
8.378906 0 m /one glyphshow
14.741211 0 m /eight glyphshow
21.103516 0 m /zero glyphshow
grestore
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
0 0 m
0 -3.5 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
346.103 33.6721 o
grestore
gsave
332.376199 19.078365 translate
0.000000 rotate
0.000000 0 m /minus glyphshow
8.378906 0 m /one glyphshow
14.741211 0 m /six glyphshow
21.103516 0 m /zero glyphshow
grestore
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
0 0 m
0 -3.5 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
404.347 33.6721 o
grestore
gsave
390.620429 19.078365 translate
0.000000 rotate
0.000000 0 m /minus glyphshow
8.378906 0 m /one glyphshow
14.741211 0 m /four glyphshow
21.103516 0 m /zero glyphshow
grestore
gsave
214.859375 5.078365 translate
0.000000 rotate
/DejaVuSans-Oblique 10.0 selectfont
0.000000 0.406250 moveto
/x glyphshow
/DejaVuSans 10.0 selectfont
5.917969 0.406250 moveto
/space glyphshow
9.096680 0.406250 moveto
/parenleft glyphshow
12.998047 0.406250 moveto
/m glyphshow
22.739258 0.406250 moveto
/m glyphshow
32.480469 0.406250 moveto
/parenright glyphshow
grestore
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
-0 0 m
-3.5 0 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
37.719 50.8677 o
grestore
gsave
18.000240 47.070840 translate
0.000000 rotate
0.000000 0 m /seven glyphshow
6.362305 0 m /zero glyphshow
grestore
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
-0 0 m
-3.5 0 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
37.719 96.096 o
grestore
gsave
18.000240 92.299152 translate
0.000000 rotate
0.000000 0 m /seven glyphshow
6.362305 0 m /two glyphshow
grestore
gsave
/o {
gsave
newpath
translate
0.8 setlinewidth
1 setlinejoin
0 setlinecap
-0 0 m
-3.5 0 l
gsave
0.000 setgray
fill
grestore
stroke
grestore
} bind def
37.719 141.324 o
grestore
gsave
18.000240 137.527464 translate
0.000000 rotate
0.000000 0 m /seven glyphshow
6.362305 0 m /four glyphshow
grestore
gsave
11.000240 74.435937 translate
90.000000 rotate
/DejaVuSans-Oblique 10.0 selectfont
0.000000 0.406250 moveto
/rho glyphshow
/DejaVuSans 10.0 selectfont
6.347656 0.406250 moveto
/space glyphshow
9.526367 0.406250 moveto
/parenleft glyphshow
13.427734 0.406250 moveto
/c glyphshow
18.925781 0.406250 moveto
/m glyphshow
28.666992 0.406250 moveto
/parenright glyphshow
grestore
1.500 setlinewidth
0.122 0.467 0.706 setrgbcolor
gsave
391.281 114.528 37.719 33.672 clipbox
55.50448 38.936602 m
57.261815 38.936602 l
stroke
grestore
gsave
391.281 114.528 37.719 33.672 clipbox
115.79737 56.684726 m
117.570009 56.684726 l
stroke
grestore
gsave
391.281 114.528 37.719 33.672 clipbox
134.046325 62.41111 m
135.850671 62.41111 l
stroke
grestore
gsave
391.281 114.528 37.719 33.672 clipbox
343.718019 123.926997 m
345.473365 123.926997 l
stroke
grestore
gsave
391.281 114.528 37.719 33.672 clipbox
409.454726 142.767063 m
411.21427 142.767063 l
stroke
grestore
gsave
391.281 114.528 37.719 33.672 clipbox
56.383147 38.877917 m
56.383147 38.995286 l
stroke
grestore
gsave
391.281 114.528 37.719 33.672 clipbox
116.68369 56.631041 m
116.68369 56.738412 l
stroke
grestore
gsave
391.281 114.528 37.719 33.672 clipbox
134.948498 62.35357 m
134.948498 62.46865 l
stroke
grestore
gsave
391.281 114.528 37.719 33.672 clipbox
344.595692 123.880074 m
344.595692 123.97392 l
stroke
grestore
gsave
391.281 114.528 37.719 33.672 clipbox
410.334498 142.720827 m
410.334498 142.813299 l
stroke
grestore
2 setlinecap
1.000 0.498 0.055 setrgbcolor
gsave
391.281 114.528 37.719 33.672 clipbox
56.383147 39.086887 m
410.334498 142.993958 l
410.334498 142.993958 l
stroke
grestore
1.000 setlinewidth
0 setlinecap
0.122 0.467 0.706 setrgbcolor
gsave
391.281 114.528 37.719 33.672 clipbox
/o {
gsave
newpath
translate
1.0 setlinewidth
1 setlinejoin
0 setlinecap
0 -3 m
0.795609 -3 1.55874 -2.683901 2.12132 -2.12132 c
2.683901 -1.55874 3 -0.795609 3 0 c
3 0.795609 2.683901 1.55874 2.12132 2.12132 c
1.55874 2.683901 0.795609 3 0 3 c
-0.795609 3 -1.55874 2.683901 -2.12132 2.12132 c
-2.683901 1.55874 -3 0.795609 -3 0 c
-3 -0.795609 -2.683901 -1.55874 -2.12132 -2.12132 c
-1.55874 -2.683901 -0.795609 -3 0 -3 c
cl
gsave
0.122 0.467 0.706 setrgbcolor
fill
grestore
stroke
grestore
} bind def
56.3831 38.9366 o
116.684 56.6847 o
134.948 62.4111 o
344.596 123.927 o
410.334 142.767 o
grestore
0.800 setlinewidth
0 setlinejoin
2 setlinecap
0.000 setgray
gsave
37.71899 33.672115 m
37.71899 148.19976 l
stroke
grestore
gsave
428.99976 33.672115 m
428.99976 148.19976 l
stroke
grestore
gsave
37.71899 33.672115 m
428.99976 33.672115 l
stroke
grestore
end
showpage

22
README.md
View File

@ -1,21 +1,23 @@
# SPSPy # SPSPy
SPSPy is a Python based package of tools for use with the Super-Enge Split-Pole Spectrograph at FSU. Much of the code here is based on Java programs originally written at Yale University by D.W. Visser, C.M. Deibel, and others. Currently the package contains spsplot, a tool aimed at informing users which states should appear at the focal plane of the SESPS, and spanc, a tool for calibrating the position spectra from the focal plane. SPSPy is a Python based package of tools for use with the Super-Enge Split-Pole Spectrograph at FSU. Much of the code here is based on Java programs originally written at Yale University by D.W. Visser, C.M. Deibel, and others. Currently the package contains SPSPlot, a tool aimed at informing users which states should appear at the focal plane of the SESPS, and SPANC, a tool for calibrating the position spectra from the focal plane.
# Depencencies and Requirements ## Depencencies and Requirements
SPSPy requires the Python packages qtpy (along with a functional Qt distriubtion for python), matplotlib, numpy, lxml, and scipy to guarantee full functionality. The most straightforward way to intstall all dependencies is by using either pip or having a distribution such as Anaconda. The requirements for running SPSPy are outlined in the requirements.txt file located in the repository. It is recommended to install these to a local virtual environment using `pip install -r requirements.txt`. For conda use the environments.yml file to create a conda environment for SPSPy. Simply run `conda env create -f environment.yml` from the SPSPy directory. conda will make a new virtual environment named spsenv with the dependencies outlined in environments.yml. If you already have an environment named spsenv or would like to change the name simply edit the first line of the enviornments.yml.
Spsplot also requires that the user have an internet connection. The recommended install for SPSPy dependencies is via pip.
## spsplot ### Creating a virtual environment with pip
This tool is intended to be used for guiding the settings of the SPS to show specific states on the focal plane detector. The user gives the program reaction information, and the program runs through the kinematics to calculate the energies of ejecta into the the SESPS. To evaluate different states, the program scrapes a list of levels from www.nndc.bnl.gov, and these levels are then passed on to the reaction handler. These levels are then shown on the screen with labels. The labels can be modified to show either the excitation energy of the state, or the kinetic energy of the ejectile. To create a virtual environment with pip in the terminal for MacOS or Linux use `python3 -m venv env` to create a local virtual environment named `env` (or whatever name you'd like), or on Windows use `py -m venv env` to do the same. To activate your new environment run `source env/bin/activate` in MacOS or Linux, or `.\env\Scripts\activate`. Now you can run the above `pip` command to install all dependencies to the virtual environment. To leave the virtual environment use the command `deactivate` in your terminal.
## spanc ## SPSPlot
This tool is intended to be used for guiding the settings of the SPS to show specific states on the focal plane detector. The user gives the program reaction information, and the program runs through the kinematics to calculate the energies of ejecta into the the SESPS. To evaluate different states, the program scrapes a list of levels from [NNDC](https://www.nndc.bnl.gov/), and these levels are then passed on to the reaction handler. These levels are then shown on the screen with labels. The labels can be modified to show either the excitation energy of the state, the kinetic energy of the ejectile, or the focal plane z-offset for a state. Note that since levels are obtained from NNDC, SPSPlot requires an internet connection. SPSPlot can also export the calculated reaction information to a csv file.
## SPANC
SPANC is the program used to calibrate SESPS focal plane spectra. It works by the user specifying a target, reaction, calibration peaks, and output peaks. The target is a description of the physical target foil used in the SPS, which is used to calculate energy loss effects. The target must contain the isotope used as the target in the reaction description. The reaction indicates to the program what type of ejecta are expected, as well as the settings of the spectrograph. Calibration data is given as centroids from a spectrum with correspoding excitation energies, as well as associated uncertainties. The calibration peaks are then fit using the scipy ODR package (see scipy ODR for more documentation). The fit is plotted, and the results are shown in a table. Additionally, residuals are plotted and shown in a table. The user can then feed the program an output peak, or a peak for which the user would like to calculate the excitation energy of a state using the calibration fit. The peak excitation energy will then be reported, with uncertainty. The user can also give a FWHM to be converted from focal plane position to energy. SPANC is the program used to calibrate SESPS focal plane spectra. It works by the user specifying a target, reaction, calibration peaks, and output peaks. The target is a description of the physical target foil used in the SPS, which is used to calculate energy loss effects. The target must contain the isotope used as the target in the reaction description. The reaction indicates to the program what type of ejecta are expected, as well as the settings of the spectrograph. Calibration data is given as centroids from a spectrum with correspoding excitation energies, as well as associated uncertainties. The calibration peaks are then fit using the scipy ODR package (see scipy ODR for more documentation). The fit is plotted, and the results are shown in a table. Additionally, residuals are plotted and shown in a table. The user can then feed the program an output peak, or a peak for which the user would like to calculate the excitation energy of a state using the calibration fit. The peak excitation energy will then be reported, with uncertainty. The user can also give a FWHM to be converted from focal plane position to energy.
# Running the tools ## Running the tools
Use `./bin/spanc` or `./bin/spsplot`. Note that they should be run from the SPSPy directory, as there are some file paths which need to be maintained. Activate the environment containing the requirements and then simply `python main.py` from the repository. This will bring up a launcher from which you can select a tool.
### Known issues ### Known issues
1. NNDC sometimes puts annoying characters in the ENDSF list; each of these "special characters" needs to be added to a list of exclusions 1. NNDC sometimes puts annoying characters in the ENDSF list; each of these "special characters" needs to be added to a list of exclusions
2. Not really an issue but with high level density reactions, spsplot becomes quite crowded. Working on implementing level removal. 2. Not really an issue but with high level density reactions, spsplot becomes quite crowded. Working on implementing level removal.
3. Debian -- mostly relevant to SPS DAQ machine, but Qt on debian running with a TightVNC instance causes a crash of the VNC server. Current fix is to implement a virtual env for a ssh into the DAQ machine from which the tools will be used, while leaving the VNC window free from the Qt related code.

3
bin/fpcheck
View File

@ -1,3 +0,0 @@
#!/bin/bash
./fpcheck/FPCheckGUI.py

3
bin/spanc
View File

@ -1,3 +0,0 @@
#!/bin/bash
./spanc/SpancGUI.py

3
bin/spsplot
View File

@ -1,3 +0,0 @@
#!/bin/bash
./spsplot/SPSPlotGUI.py

51
environment.yml Normal file
View File

@ -0,0 +1,51 @@
name: spsenv
channels:
- defaults
dependencies:
- _libgcc_mutex=0.1=main
- _openmp_mutex=5.1=1_gnu
- bzip2=1.0.8=h7b6447c_0
- ca-certificates=2022.10.11=h06a4308_0
- certifi=2022.9.24=py310h06a4308_0
- ld_impl_linux-64=2.38=h1181459_1
- libffi=3.4.2=h6a678d5_6
- libgcc-ng=11.2.0=h1234567_1
- libgomp=11.2.0=h1234567_1
- libstdcxx-ng=11.2.0=h1234567_1
- libuuid=1.41.5=h5eee18b_0
- ncurses=6.3=h5eee18b_3
- openssl=1.1.1s=h7f8727e_0
- pip=22.2.2=py310h06a4308_0
- python=3.10.8=h7a1cb2a_1
- readline=8.2=h5eee18b_0
- setuptools=65.5.0=py310h06a4308_0
- sqlite=3.40.0=h5082296_0
- tk=8.6.12=h1ccaba5_0
- tzdata=2022f=h04d1e81_0
- wheel=0.37.1=pyhd3eb1b0_0
- xz=5.2.8=h5eee18b_0
- zlib=1.2.13=h5eee18b_0
- pip:
- charset-normalizer==2.1.1
- contourpy==1.0.6
- cycler==0.11.0
- fonttools==4.38.0
- idna==3.4
- kiwisolver==1.4.4
- lxml==4.9.1
- matplotlib==3.6.2
- numpy==1.23.5
- packaging==21.3
- pillow==9.3.0
- pycatima==1.71
- pyparsing==3.0.9
- pyqtdarktheme==1.2.1
- pyside6==6.4.1
- pyside6-addons==6.4.1
- pyside6-essentials==6.4.1
- python-dateutil==2.8.2
- requests==2.28.1
- scipy==1.9.3
- shiboken6==6.4.1
- six==1.16.0
- urllib3==1.26.13

0
etc/mass.txt → etc/amdc2016_mass.txt
View File

179
fpcheck/FPCheckGUI.py
View File

@ -1,179 +0,0 @@
#!/usr/bin/env python3
from Reaction import Reaction
from qtpy.QtWidgets import QApplication, QWidget, QMainWindow
from qtpy.QtWidgets import QLabel, QMenuBar, QAction
from qtpy.QtWidgets import QHBoxLayout, QVBoxLayout, QGridLayout, QGroupBox
from qtpy.QtWidgets import QPushButton, QButtonGroup, QRadioButton
from qtpy.QtWidgets import QSpinBox, QDoubleSpinBox, QComboBox
from qtpy.QtWidgets import QDialog, QFileDialog, QDialogButtonBox
from qtpy.QtWidgets import QTableWidget, QTableWidgetItem
from qtpy.QtWidgets import QLineEdit, QTabWidget, QFormLayout
from qtpy.QtCore import Signal
import sys
class ReactionDialog(QDialog):
new_reaction = Signal(Reaction)
update_reaction = Signal(float, float, float, str)
def __init__(self, parent=None, rxn=None):
super().__init__(parent)
self.setWindowTitle("Add A Reaction")
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
if rxn is not None:
self.buttonBox.accepted.connect(self.SendReactionUpdate)
else:
self.buttonBox.accepted.connect(self.SendReaction)
self.buttonBox.rejected.connect(self.reject)
self.layout = QVBoxLayout()
self.setLayout(self.layout)
self.CreateReactionInputs()
if rxn is not None:
self.SetInitialValues(rxn)
self.layout.addWidget(self.buttonBox)
def SendReaction(self) :
rxn = Reaction(self.ztInput.value(), self.atInput.value(), self.zpInput.value(),self.apInput.value(),self.zeInput.value(),self.aeInput.value(), self.bkeInput.value(), self.thetaInput.value(), self.bfieldInput.value())
self.new_reaction.emit(rxn)
def SendReactionUpdate(self):
self.update_reaction.emit(self.bkeInput.value(), self.thetaInput.value(), self.bfieldInput.value(), self.rxnKey)
def CreateReactionInputs(self) :
self.nucleiGroupBox = QGroupBox("Reaction Nuclei",self)
inputLayout = QFormLayout()
self.ztInput = QSpinBox(self.nucleiGroupBox)
self.ztInput.setRange(1, 110)
self.atInput = QSpinBox(self.nucleiGroupBox)
self.atInput.setRange(1,270)
self.zpInput = QSpinBox(self.nucleiGroupBox)
self.zpInput.setRange(1, 110)
self.apInput = QSpinBox(self.nucleiGroupBox)
self.apInput.setRange(1,270)
self.zeInput = QSpinBox(self.nucleiGroupBox)
self.zeInput.setRange(1, 110)
self.aeInput = QSpinBox(self.nucleiGroupBox)
self.aeInput.setRange(1,270)
inputLayout.addRow("ZT",self.ztInput)
inputLayout.addRow("AT",self.atInput)
inputLayout.addRow("ZP",self.zpInput)
inputLayout.addRow("AP",self.apInput)
inputLayout.addRow("ZE",self.zeInput)
inputLayout.addRow("AE",self.aeInput)
self.parameterGroupBox = QGroupBox("Reaction Parameters", self)
parameterLayout = QFormLayout()
self.bkeInput = QDoubleSpinBox(self.parameterGroupBox)
self.bkeInput.setRange(0.0, 40.0)
self.bkeInput.setDecimals(4)
self.thetaInput = QDoubleSpinBox(self.parameterGroupBox)
self.thetaInput.setRange(0.0, 180.0)
self.thetaInput.setDecimals(4)
self.bfieldInput = QDoubleSpinBox(self.parameterGroupBox)
self.bfieldInput.setRange(0.0, 16.0)
self.bfieldInput.setDecimals(6)
parameterLayout.addRow("Beam KE(Mev)",self.bkeInput)
parameterLayout.addRow("Theta(deg)",self.thetaInput)
parameterLayout.addRow("Bfield(kG)",self.bfieldInput)
self.nucleiGroupBox.setLayout(inputLayout)
self.parameterGroupBox.setLayout(parameterLayout)
self.layout.addWidget(self.nucleiGroupBox)
self.layout.addWidget(self.parameterGroupBox)
def SetInitialValues(self, rxn):
self.ztInput.setValue(rxn.Target.Z)
self.ztInput.setEnabled(False)
self.atInput.setValue(rxn.Target.A)
self.atInput.setEnabled(False)
self.zpInput.setValue(rxn.Projectile.Z)
self.zpInput.setEnabled(False)
self.apInput.setValue(rxn.Projectile.A)
self.apInput.setEnabled(False)
self.zeInput.setValue(rxn.Ejectile.Z)
self.zeInput.setEnabled(False)
self.aeInput.setValue(rxn.Ejectile.A)
self.aeInput.setEnabled(False)
self.bkeInput.setValue(rxn.BKE)
self.thetaInput.setValue(rxn.Theta/rxn.DEG2RAD)
self.bfieldInput.setValue(rxn.Bfield)
class FPCheckGUI(QMainWindow):
def __init__(self, parent=None):
super().__init__(parent)
self.setWindowTitle("FPCheck")
self.reactions = [] #data
self.centralLayout = QVBoxLayout()
self.centralWidget = QWidget(self)
self.setCentralWidget(self.centralWidget)
self.centralWidget.setLayout(self.centralLayout)
self.CreateMenus()
self.CreateTable()
self.show()
def CreateMenus(self):
self.fileMenu = self.menuBar().addMenu("File")
saveAction = QAction("Save", self)
loadAction = QAction("Load", self)
exitAction = QAction("Exit", self)
self.fileMenu.addAction(saveAction)
self.fileMenu.addAction(loadAction)
self.fileMenu.addAction(exitAction)
def CreateTable(self):
self.reactionTable = QTableWidget(self)
self.reactionTable.setColumnCount(8)
self.reactionTable.setHorizontalHeaderLabels(["Target","Projectile","Ejectile","Residual","Beam KE(MeV)","BField(kG)","Angle(deg)","FP ZOffset (cm)"])
self.centralLayout.addWidget(self.reactionTable)
self.reactionTable.resizeColumnsToContents()
self.addButton = QPushButton("Add", self)
self.addButton.clicked.connect(self.HandleAddReaction)
self.centralLayout.addWidget(self.addButton)
def HandleAddReaction(self):
rxnDia = ReactionDialog(self)
rxnDia.new_reaction.connect(self.AddReaction)
rxnDia.exec()
def UpdateTable(self):
self.reactionTable.setRowCount(len(self.reactions))
curRow = 0
for rxn in self.reactions:
self.reactionTable.setItem(curRow, 0, QTableWidgetItem(rxn.Target.Symbol))
self.reactionTable.setItem(curRow, 1, QTableWidgetItem(rxn.Projectile.Symbol))
self.reactionTable.setItem(curRow, 2, QTableWidgetItem(rxn.Ejectile.Symbol))
self.reactionTable.setItem(curRow, 3, QTableWidgetItem(rxn.Residual.Symbol))
self.reactionTable.setItem(curRow, 4, QTableWidgetItem(str(rxn.BKE)))
self.reactionTable.setItem(curRow, 5, QTableWidgetItem(str(rxn.Bfield)))
self.reactionTable.setItem(curRow, 6, QTableWidgetItem(str(rxn.Theta)))
self.reactionTable.setItem(curRow, 7, QTableWidgetItem(str(rxn.GetFocalPlaneZOffset())))
curRow += 1
self.reactionTable.resizeColumnsToContents()
self.reactionTable.resizeRowsToContents()
def AddReaction(self, rxn):
self.reactions.append(rxn)
self.UpdateTable()
def main():
myapp = QApplication(sys.argv)
window = FPCheckGUI()
sys.exit(myapp.exec_())
if __name__ == '__main__':
main()

70
fpcheck/NucData.py
View File

@ -1,70 +0,0 @@
#!/usr/bin/env python3
import numpy as np
import requests
import lxml.html as xhtml
class MassTable:
def __init__(self):
file = open("./etc/mass.txt","r")
self.mtable = {}
u2mev = 931.4940954
me = 0.000548579909
self.etable = {}
file.readline()
file.readline()
for line in file:
entries = line.split()
n = entries[0]
z = entries[1]
a = entries[2]
element = entries[3]
massBig = float(entries[4])
massSmall = float(entries[5])
key = '('+z+','+a+')'
value = ((massBig+massSmall*1e-6) - float(z)*me)*u2mev
self.mtable[key] = value
self.etable[key] = element
file.close()
def GetMass(self, z, a):
key = '('+str(z)+','+str(a)+')'
if key in self.mtable:
return self.mtable[key]
else:
return 0
def GetSymbol(self, z, a):
key = '('+str(z)+','+str(a)+')'
if key in self.etable:
return str(a)+self.etable[key]
else:
return 'none'
Masses = MassTable()
def GetExcitations(symbol):
levels = np.array(np.empty(0))
text = ''
site = requests.get("https://www.nndc.bnl.gov/nudat2/getdatasetClassic.jsp?nucleus="+symbol+"&unc=nds")
contents = xhtml.fromstring(site.content)
tables = contents.xpath("//table")
rows = tables[2].xpath("./tr")
for row in rows[1:-2]:
entries = row.xpath("./td")
if len(entries) != 0:
entry = entries[0]
data = entry.xpath("./a")
if len(data) == 0:
text = entry.text
else:
text = data[0].text
text = text.replace('?', '')
text = text.replace('\xa0\xa0≈','')
levels = np.append(levels, float(text)/1000.0)
return levels

67
fpcheck/Reaction.py
View File

@ -1,67 +0,0 @@
#!/usr/bin/env python3
from NucData import Masses
import numpy as np
class Nucleus:
def __init__(self, z, a):
self.Z = z
self.A = a
self.Symbol = Masses.GetSymbol(self.Z, self.A)
self.GSMass = Masses.GetMass(self.Z, self.A)
def Minus(self, rhs):
final_Z = self.Z - rhs.Z
final_A = self.A - rhs.A
if final_A < 0 or final_Z < 0:
print("Illegal minus operation on Nuclei!")
return Nucleus(0,0)
else:
return Nucleus(final_Z, final_A)
def Plus(self, rhs):
return Nucleus(self.Z + rhs.Z, self.A + rhs.A)
class Reaction:
DEG2RAD = np.pi/180.0 #degrees to radians
C = 299792458 #speed of light m/s
QBRHO2P = 1.0E-9*C #Converts qbrho to p (kG*cm -> MeV/c)
DISP = 1.96
MAG = 0.39
def __init__(self, zt, at, zp, ap, ze, ae, beamKE, theta, bfield):
self.Target = Nucleus(zt, at)
self.Projectile = Nucleus(zp, ap)
self.Ejectile = Nucleus(ze, ae)
self.Residual = (self.Target.Plus(self.Projectile)).Minus(self.Ejectile)
self.BKE = beamKE
self.Theta = theta * self.DEG2RAD
self.Bfield = bfield
self.Name = self.Target.Symbol +"("+ self.Projectile.Symbol +","+ self.Ejectile.Symbol +")"+ self.Residual.Symbol
def GetEjectileKineticEnergy(self, Elevel) :
Q = self.Target.GSMass + self.Projectile.GSMass - (self.Ejectile.GSMass + self.Residual.GSMass + Elevel)
Ethresh = -Q*(self.Ejectile.GSMass+self.Residual.GSMass)/(self.Ejectile.GSMass + self.Residual.GSMass - self.Projectile.GSMass)
if self.BKE < Ethresh:
return 0.0
term1 = np.sqrt(self.Projectile.GSMass*self.Ejectile.GSMass*self.BKE)/(self.Ejectile.GSMass + self.Residual.GSMass)*np.cos(self.Theta)
term2 = (self.BKE*(self.Residual.GSMass - self.Projectile.GSMass) + self.Residual.GSMass*Q)/(self.Ejectile.GSMass + self.Residual.GSMass)
ke1 = term1 + np.sqrt(term1**2.0 + term2)
ke2 = term1 - np.sqrt(term1**2.0 + term2)
if ke1 > 0:
return ke1**2.0
else :
return ke2**2.0
def GetFocalPlaneZOffset(self):
ejectKE = self.GetEjectileKineticEnergy(0.0)
ejectP = np.sqrt(ejectKE**2.0 + 2.0*ejectKE*self.Ejectile.GSMass)
rho = ejectP/(self.QBRHO2P*self.Bfield*self.Ejectile.Z)
K = np.sqrt(self.Projectile.GSMass*self.Ejectile.GSMass*self.BKE/ejectKE)*np.sin(self.Theta)
K /= self.Ejectile.GSMass + self.Residual.GSMass - np.sqrt(self.Projectile.GSMass*self.Ejectile.GSMass*self.BKE/ejectKE)*np.cos(self.Theta)
return -1.0*K*rho*self.DISP*self.MAG
def ChangeReactionParameters(self, bke, theta, bfield) :
self.BKE = bke
self.Theta = theta*self.DEG2RAD
self.Bfield = bfield

9
main.py Normal file
View File

@ -0,0 +1,9 @@
import os
os.environ["QT_API"] = "pyside6"
from spspy.SPSPlotUI import run_spsplot_ui
from spspy.SpancUI import run_spanc_ui
from spspy.Launcher import run_launcher
#run_spsplot_ui()
#run_spanc_ui()
run_launcher()

8
requirements.txt Normal file
View File

@ -0,0 +1,8 @@
lxml==4.9.1
matplotlib==3.6.2
numpy==1.23.5
pycatima==1.71
pyqtdarktheme==1.2.1
PySide6==6.4.1
requests==2.28.1
scipy==1.9.3

213
spanc/EnergyLoss.py
View File

@ -1,213 +0,0 @@
#!/usr/bin/env python3
import numpy as np
import EnergyLossData as edata
from NucData import Masses
class EnergyLoss:
MAX_FRACTIONAL_STEP=0.001
MAX_H_E_PER_U=100000.0
AVOGADRO=0.60221367
MEV2U=1.0/931.4940954
def __init__(self):
self.ZP = 0
self.AP = 0
self.MP = 0
self.ZT = np.zeros(0)
self.AT = np.zeros(0)
self.Stoich = np.zeros(0)
self.illegalFlag = True
def SetTargetData(self, zt, at, stoich):
self.ZT = zt
self.AT = at
total = np.sum(stoich)
self.Stoich = stoich/total
for z in self.ZT:
if z >= edata.MaxZ:
self.illegalFlag = True
return
self.illegalFlag = False
def GetEnergyLoss(self, zp, ap, e_initial, thickness):
if self.illegalFlag:
print("Unable to get energy loss with unset target data... returning 0")
return 0.0
if self.ZP != zp:
self.ZP = zp
self.AP = ap
self.MP = Masses.GetMass(self.ZP, self.AP)*self.MEV2U
e_final = e_initial
x_traversed = 0
x_step = 0.25*thickness
e_step = self.GetTotalStoppingPower(e_final)*x_step/1000.0
if thickness == 0.0:
return 0.0
go = True
while go:
if e_step/e_final > self.MAX_FRACTIONAL_STEP:
x_step *= 0.5
e_step = self.GetTotalStoppingPower(e_final)*x_step/1000.0
elif x_step+x_traversed >= thickness:
go = False
x_step = thickness - x_traversed #get valid portion of last chunk
e_final -= self.GetTotalStoppingPower(e_final)*x_step/1000.0
if e_final <= 0.0:
return e_initial
else:
x_traversed += x_step
e_step = self.GetTotalStoppingPower(e_final)*x_step/1000.0
e_final -= e_step
if e_final <= 0.0:
return e_initial
return e_initial - e_final
def GetReverseEnergyLoss(self, zp, ap, e_final, thickness):
if self.illegalFlag:
print("Unable to get energy loss with unset target data... returning 0")
return 0.0
if self.ZP != zp:
self.ZP = zp
self.AP = ap
self.MP = Masses.GetMass(self.ZP, self.AP)*self.MEV2U
e_initial = e_final
x_traversed = 0
x_step = 0.25*thickness
e_step = self.GetTotalStoppingPower(e_initial)*x_step/1000.0
if thickness == 0.0:
return 0.0
go = True
while go:
if e_step/e_final > self.MAX_FRACTIONAL_STEP:
x_step *= 0.5
e_step = self.GetTotalStoppingPower(e_initial)*x_step/1000.0
elif x_step+x_traversed >= thickness:
go = False
x_step = thickness - x_traversed #get valid portion of last chunk
e_initial += self.GetTotalStoppingPower(e_initial)*x_step/1000.0
else:
x_traversed += x_step
e_step = self.GetTotalStoppingPower(e_initial)*x_step/1000.0
e_final += e_step
if e_final <= 0.0:
return e_initial
return abs(e_initial - e_final)
def GetTotalStoppingPower(self, energy):
return self.GetElectronicStoppingPower(energy)+self.GetNuclearStoppingPower(energy)
def GetElectronicStoppingPower(self, energy):
e_per_u = energy*1000.0/self.MP
values = np.zeros(len(self.ZT))
if e_per_u > self.MAX_H_E_PER_U:
print("Bombarding energy exceeds maximum allowed value for energy loss! Returning 0")
return 0.0
elif e_per_u > 1000.0:
for i in range(len(self.ZT)):
values[i] = self.Hydrogen_dEdx_High(e_per_u, energy, self.ZT[i])
elif e_per_u > 10.0:
for i in range(len(self.ZT)):
values[i] = self.Hydrogen_dEdx_Med(e_per_u, self.ZT[i])
elif e_per_u > 0.0:
for i in range(len(self.ZT)):
values[i] = self.Hydrogen_dEdx_Low(e_per_u, self.ZT[i])
else:
print("Negative energy encountered at EnergyLoss::GetElectronicStoppingPower! Returning 0")
return 0.0
if self.ZP > 1:
for i in range(len(values)):
values[i] *= self.CalculateEffectiveChargeRatio(e_per_u, self.ZT[i])
stopping_total = np.sum(values*self.Stoich)
conv_factor = 0.0
for i in range(len(self.ZT)):
conv_factor += self.Stoich[i]*edata.NaturalMass[self.ZT[i]]
stopping_total *= self.AVOGADRO/conv_factor
return stopping_total
def GetNuclearStoppingPower(self, energy):
e = energy*1000.0
stopping_total = 0.0
for i in range(len(self.ZT)):
zt = self.ZT[i]
mt = edata.NaturalMass[self.ZT[i]]
x = (self.MP + mt) * np.sqrt(self.ZP**(2.0/3.0) + zt**(2.0/3.0))
epsilon = 32.53*mt*e/(self.ZP*zt*x)
sn = 8.462*(0.5*np.log(1.0+epsilon)/(epsilon+0.10718*(epsilon**0.37544)))*self.ZP*zt*self.MP/x
conversion_factor = self.AVOGADRO/mt
stopping_total += sn*conversion_factor*self.Stoich[i]
return stopping_total
def Hydrogen_dEdx_Low(self, e_per_u, zt):
return np.sqrt(e_per_u)*edata.HydrogenCoeff[zt][0]
def Hydrogen_dEdx_Med(self, e_per_u, zt):
x = edata.HydrogenCoeff[zt][1]*e_per_u**0.45
y = edata.HydrogenCoeff[zt][2]/e_per_u * np.log(1.0 + edata.HydrogenCoeff[zt][3]/e_per_u + edata.HydrogenCoeff[zt][4]*e_per_u)
return x*y/(x+y)
def Hydrogen_dEdx_High(self, e_per_u, energy, zt):
beta_sq = energy * (energy + 2.0*self.MP/self.MEV2U)/((energy + self.MP/self.MEV2U)**2.0)
alpha = edata.HydrogenCoeff[zt][5]/beta_sq
epsilon = edata.HydrogenCoeff[zt][6]*beta_sq/(1.0 - beta_sq) - beta_sq - edata.HydrogenCoeff[zt][7]
for i in range(1,5):
epsilon += edata.HydrogenCoeff[zt][7+i]*(np.log(e_per_u))**float(i)
return alpha * np.log(epsilon)
def CalculateEffectiveChargeRatio(self, e_per_u, zt):
z_ratio=0
if self.ZP == 2:
ln_epu = np.log(e_per_u)
gamma = 1.0+(0.007+0.00005*zt)*np.exp(-1.0*(7.6-ln_epu)**2.0)
alpha = 0.7446 + 0.1429*ln_epu + 0.01562*ln_epu**2.0 - 0.00267*ln_epu**3.0 + 1.338e-6*ln_epu**8.0
z_ratio = gamma*(1.0-np.exp(-alpha))*2.0
elif self.ZP == 3:
ln_epu = np.log(e_per_u)
gamma = 1.0+(0.007+0.00005*zt)*np.exp(-1.0*(7.6-ln_epu)**2.0)
alpha = 0.7138+0.002797*e_per_u+1.348e-6*e_per_u**2.0
z_ratio = gamma*(1-np.exp(-alpha))*3.0
else:
B = 0.886*(e_per_u/25.0)**0.5/(self.ZP**(2.0/3.0))
A = B + 0.0378*np.sin(np.pi/2.0*B)
z_ratio = (1.0 - np.exp(-A)*(1.034-0.1777*np.exp(-0.08114*self.ZP)))*self.ZP
return z_ratio*z_ratio
def main():
targetA = np.array([12])
targetZ = np.array([6])
targetS = np.array([1])
beamKE = 16.0
thickness = 20.0
eloss = EnergyLoss()
eloss.SetTargetData(targetZ, targetA, targetS)
print("Testing various cases for energy loss. Using 12C target with 20 ug/cm^2 thickness. Compare to values given by LISE++ or SRIM")
result = eloss.GetEnergyLoss(1, 1, beamKE, 20.0)
print("Case 1: ZP = 1, AP=1, Beam energy = 16 MeV -> Resulting energy loss = ", result, " MeV")
beamKE = 1.0
result = eloss.GetEnergyLoss(1, 1, beamKE, 20.0)
print("Case 2: ZP = 1, AP=1, Beam energy = 1.0 MeV -> Resulting energy loss = ", result, " MeV")
beamKE = 0.1
result = eloss.GetEnergyLoss(1, 1, beamKE, 20.0)
print("Case 3: ZP = 1, AP=1, Beam energy = 0.1 MeV -> Resulting energy loss = ", result, " MeV")
beamKE = 0.01
result = eloss.GetEnergyLoss(1, 1, beamKE, 20.0)
print("Case 4: ZP = 1, AP=1, Beam energy = 0.01 MeV -> Resulting energy loss = ", result, " MeV")
beamKE = 24.0
result = eloss.GetEnergyLoss(2, 4, beamKE, 20.0)
print("Case 5: ZP = 2, AP=4, Beam energy = 24.0 MeV -> Resulting energy loss = ", result, " MeV")
beamKE = 24.0
result = eloss.GetEnergyLoss(3, 6, beamKE, 20.0)
print("Case 6: ZP = 3, AP=6, Beam energy = 24 MeV -> Resulting energy loss = ", result, " MeV")
print("Finished.")
if __name__ == '__main__':
main()

122
spanc/EnergyLossData.py
View File

@ -1,122 +0,0 @@
#!/usr/bin/env python3
import numpy as np
MaxZ = 93
NaturalMass = np.array([0, 1.00797, 4.0026, 6.939, 9.0122, 10.818,
12.01115, 14.0067, 15.9994, 18.99984, 20.183,
22.9898, 24.312, 26.9815, 28.086, 30.9738,
32.064, 35.453, 39.948, 39.102, 40.08,
44.956, 47.90, 50.942, 51.996, 54.938,
55.847, 58.933, 58.71, 63.54, 65.37,
69.72, 72.59, 74.922, 78.96, 79.909,
83.80, 85.47, 87.62, 88.909, 91.22,
92.906, 95.94, 98., 101.07, 102.905,
106.4, 107.87, 112.4, 114.82, 118.69,
121.75, 127.60, 126.904, 131.3, 132.905,
137.34, 138.91, 140.12, 140.907, 144.24,
146., 150.35, 151.96, 157.25, 158.924,
162.50, 164.93, 167.26, 168.934, 173.04,
174.97, 178.49, 180.948, 183.85, 186.2,
190.2, 192.2, 195.09, 196.967, 200.59,
204.37, 207.19, 208.98, 209., 210.,
222., 223., 226., 227., 232.038,
231., 238.03])
HydrogenCoeff = np.array([
[0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.],#Blank
[1.262,1.44,242.6,1.2E4,0.1159,0.0005099,5.436E4,-5.052,2.049,-0.3044,0.01966,-0.0004659],#H
[1.229,1.397,484.5,5873,0.05225,0.00102,2.451E4,-2.158,0.8278,-0.1172,0.007259,-0.000166],#He
[1.411,1.6,725.6,3013,0.04578,0.00153,2.147E4,-0.5831,0.562,-0.1183,0.009298,-0.0002498],#Li
[2.248,2.59,966,153.8,0.03475,0.002039,1.63E4,0.2779,0.1745,-0.05684,0.005155,-0.0001488],#Be
[2.474,2.815,1206,1060,0.02855,0.002549,1.345E4,-2.445,1.283,-0.2205,0.0156,-0.000393],#B
[2.631,2.989,1445,957.2,0.02819,0.003059,1.322E4,-4.38,2.044,-0.3283,0.02221,-0.0005417],#C
[2.954,3.35,1683,1900,0.02513,0.003569,1.179E4,-5.054,2.325,-0.3713,0.02506,-0.0006109],#N
[2.652,3,1920,2000,0.0223,0.004079,1.046E4,-6.734,3.019,-0.4748,0.03171,-0.0007669],#O
[2.085,2.352,2157,2634,0.01816,0.004589,8517,-5.571,2.449,-0.3781,0.02483,-0.0005919],#F
[1.951,2.199,2393,2699,0.01568,0.005099,7353,-4.408,1.879,-0.2814,0.01796,-0.0004168],#Ne
[2.542,2.869,2628,1854,0.01472,0.005609,6905,-4.959,2.073,-0.3054,0.01921,-0.0004403],#Na
[3.792,4.293,2862,1009,0.01397,0.006118,6551,-5.51,2.266,-0.3295,0.02047,-0.0004637],#Mg
[4.154,4.739,2766,164.5,0.02023,0.006628,6309,-6.061,2.46,-0.3535,0.02173,-0.0004871],#Al
[4.15,4.7,3329,550,0.01321,0.007138,6194,-6.294,2.538,-0.3628,0.0222,-0.0004956],#Si
[3.232,3.647,3561,1560,0.01267,0.007648,5942,-6.527,2.616,-0.3721,0.02267,-0.000504],#P
[3.447,3.891,3792,1219,0.01211,0.008158,5678,-6.761,2.694,-0.3814,0.02314,-0.0005125],#S
[5.047,5.714,4023,878.6,0.01178,0.008668,5524,-6.994,2.773,-0.3907,0.02361,-0.0005209],#Cl
[5.731,6.5,4253,530,0.01123,0.009178,5268,-7.227,2.851,-0.4,0.02407,-0.0005294],#Ar
[5.151,5.833,4482,545.7,0.01129,0.009687,5295,-7.44,2.923,-0.4094,0.02462,-0.0005411],#K
[5.521,6.252,4710,553.3,0.01112,0.0102,5214,-7.653,2.995,-0.4187,0.02516,-0.0005529],#Ca
[5.201,5.884,4938,560.9,0.009995,0.01071,4688,-8.012,3.123,-0.435,0.02605,-0.0005707],#Sc.....
[4.862,5.496,5165,568.5,0.009474,0.01122,4443,-8.371,3.251,-0.4513,0.02694,-0.0005886],
[4.48,5.055,5391,952.3,0.009117,0.01173,4276,-8.731,3.379,-0.4676,0.02783,-0.0006064],
[3.983,4.489,5616,1336,0.008413,0.01224,3946,-9.09,3.507,-0.4838,0.02872,-0.0006243],
[3.469,3.907,5725,1461,0.008829,0.01275,3785,-9.449,3.635,-0.5001,0.02961,-0.0006421],
[3.519,3.963,6065,1243,0.007782,0.01326,3650,-9.809,3.763,-0.5164,0.0305,-0.00066],
[3.14,3.535,6288,1372,0.007361,0.01377,3453,-10.17,3.891,-0.5327,0.03139,-0.0006779],
[3.553,4.004,6205,555.1,0.008763,0.01428,3297,-10.53,4.019,-0.549,0.03229,-0.0006957],
[3.696,4.175,4673,387.8,0.02188,0.01479,3174,-11.18,4.252,-0.5791,0.03399,-0.0007314],
[4.21,4.75,6953,295.2,0.006809,0.0153,3194,-11.57,4.394,-0.598,0.03506,-0.0007537],
[5.041,5.697,7173,202.6,0.006725,0.01581,3154,-11.95,4.537,-0.6169,0.03613,-0.0007759],
[5.554,6.3,6496,110,0.009689,0.01632,3097,-12.34,4.68,-0.6358,0.03721,-0.0007981],
[5.323,6.012,7611,292.5,0.006447,0.01683,3024,-12.72,4.823,-0.6547,0.03828,-0.0008203],
[5.874,6.656,7395,117.5,0.007684,0.01734,3006,-13.11,4.965,-0.6735,0.03935,-0.0008425],
[5.611,6.335,8046,365.2,0.006244,0.01785,2928,-13.4,5.083,-0.6906,0.04042,-0.0008675],
[6.411,7.25,8262,220,0.006087,0.01836,2855,-13.69,5.2,-0.7076,0.0415,-0.0008925],
[5.694,6.429,8478,292.9,0.006087,0.01886,2855,-13.92,5.266,-0.714,0.04173,-0.0008943],
[6.339,7.159,8693,330.3,0.006003,0.01937,2815,-14.14,5.331,-0.7205,0.04196,-0.0008962],
[6.407,7.234,8907,367.8,0.005889,0.01988,2762,-14.36,5.397,-0.7269,0.04219,-0.000898],
[6.734,7.603,9120,405.2,0.005765,0.02039,2704,-14.59,5.463,-0.7333,0.04242,-0.0008998],
[6.902,7.791,9333,442.7,0.005587,0.0209,2621,-16.22,6.094,-0.8225,0.04791,-0.001024],
[6.425,7.248,9545,480.2,0.005367,0.02141,2517,-17.85,6.725,-0.9116,0.05339,-0.001148],
[6.799,7.671,9756,517.6,0.005315,0.02192,2493,-17.96,6.752,-0.9135,0.05341,-0.001147],
[6.108,6.887,9966,555.1,0.005151,0.02243,2416,-18.07,6.779,-0.9154,0.05342,-0.001145],
[5.924,6.677,1.018E4,592.5,0.004919,0.02294,2307,-18.18,6.806,-0.9173,0.05343,-0.001143],
[5.238,5.9,1.038E4,630,0.004758,0.02345,2231,-18.28,6.833,-0.9192,0.05345,-0.001142],
[5.623,6.354,7160,337.6,0.01394,0.02396,2193,-18.39,6.86,-0.9211,0.05346,-0.00114],
[5.814,6.554,1.08E4,355.5,0.004626,0.02447,2170,-18.62,6.915,-0.9243,0.0534,-0.001134],
[6.23,7.024,1.101E4,370.9,0.00454,0.02498,2129,-18.85,6.969,-0.9275,0.05335,-0.001127],
[6.41,7.227,1.121E4,386.4,0.004474,0.02549,2099,-19.07,7.024,-0.9308,0.05329,-0.001121],
[7.5,8.48,8608,348,0.009074,0.026,2069,-19.57,7.225,-0.9603,0.05518,-0.001165],
[6.979,7.871,1.162E4,392.4,0.004402,0.02651,2065,-20.07,7.426,-0.9899,0.05707,-0.001209],
[7.725,8.716,1.183E4,394.8,0.004376,0.02702,2052,-20.56,7.627,-1.019,0.05896,-0.001254],
[8.231,9.289,1.203E4,397.3,0.004384,0.02753,2056,-21.06,7.828,-1.049,0.06085,-0.001298],
[7.287,8.218,1.223E4,399.7,0.004447,0.02804,2086,-20.4,7.54,-1.004,0.05782,-0.001224],
[7.899,8.911,1.243E4,402.1,0.004511,0.02855,2116,-19.74,7.252,-0.9588,0.05479,-0.001151],
[8.041,9.071,1.263E4,404.5,0.00454,0.02906,2129,-19.08,6.964,-0.9136,0.05176,-0.001077],
[7.489,8.444,1.283E4,406.9,0.00442,0.02957,2073,-18.43,6.677,-0.8684,0.04872,-0.001003],
[7.291,8.219,1.303E4,409.3,0.004298,0.03008,2016,-17.77,6.389,-0.8233,0.04569,-0.0009292],
[7.098,8,1.323E4,411.8,0.004182,0.03059,1962,-17.11,6.101,-0.7781,0.04266,-0.0008553],
[6.91,7.786,1.343E4,414.2,0.00405,0.0311,1903,-16.45,5.813,-0.733,0.03963,-0.0007815],
[6.728,7.58,1.362E4,416.6,0.003976,0.03161,1865,-15.79,5.526,-0.6878,0.0366,-0.0007077],
[6.551,7.38,1.382E4,419,0.003877,0.03212,1819,-15.13,5.238,-0.6426,0.03357,-0.0006339],
[6.739,7.592,1.402E4,421.4,0.003863,0.03263,1812,-14.47,4.95,-0.5975,0.03053,-0.0005601],
[6.212,6.996,1.421E4,423.9,0.003725,0.03314,1747,-14.56,4.984,-0.6022,0.03082,-0.0005668],
[5.517,6.21,1.44E4,426.3,0.003632,0.03365,1703,-14.65,5.018,-0.6069,0.03111,-0.0005734],
[5.219,5.874,1.46E4,428.7,0.003498,0.03416,1640,-14.74,5.051,-0.6117,0.03141,-0.0005801],
[5.071,5.706,1.479E4,433,0.003405,0.03467,1597,-14.83,5.085,-0.6164,0.0317,-0.0005867],
[4.926,5.542,1.498E4,433.5,0.003342,0.03518,1567,-14.91,5.119,-0.6211,0.03199,-0.0005933],
[4.787, 5.386,1.517E4,435.9,0.003292,0.03569,1544,-15,5.153,-0.6258,0.03228,-0.0006],
[4.893, 5.505,1.536E4,438.4,0.003243,0.0362,1521,-15.09,5.186,-0.6305,0.03257,-0.0006066],
[5.028, 5.657,1.555E4,440.8,0.003195,0.03671,1499,-15.18,5.22,-0.6353,0.03286,-0.0006133],
[4.738, 5.329,1.574E4,443.2,0.003186,0.03722,1494,-15.27,5.254,-0.64,0.03315,-0.0006199],
[4.574, 5.144,1.593E4,442.4,0.003144,0.03773,1475,-15.67,5.392,-0.6577,0.03418,-0.0006426],
[5.2, 5.851,1.612E4,441.6,0.003122,0.03824,1464,-16.07,5.529,-0.6755,0.03521,-0.0006654],
[5.07, 5.704,1.63E4,440.9,0.003082,0.03875,1446,-16.47,5.667,-0.6932,0.03624,-0.0006881],
[4.945, 5.563,1.649E4,440.1,0.002965,0.03926,1390,-16.88,5.804,-0.711,0.03727,-0.0007109],
[4.476, 5.034,1.667E4,439.3,0.002871,0.03977,1347,-17.28,5.942,-0.7287,0.0383,-0.0007336],
[4.856, 5.46,1.832E4,438.5,0.002542,0.04028,1354,-17.02,5.846,-0.7149,0.0374,-0.0007114],
[4.308, 4.843,1.704E4,487.8,0.002882,0.04079,1352,-17.84,6.183,-0.7659,0.04076,-0.0007925],
[4.723, 5.311,1.722E4,537,0.002913,0.0413,1366,-18.66,6.52,-0.8169,0.04411,-0.0008737],
[5.319, 5.982,1.74E4,586.3,0.002871,0.04181,1347,-19.48,6.857,-0.8678,0.04747,-0.0009548],
[5.956, 6.7,1.78E4,677,0.00266,0.04232,1336,-19.55,6.871,-0.8686,0.04748,-0.0009544],
[6.158, 6.928,1.777E4,586.3,0.002812,0.04283,1319,-19.62,6.884,-0.8694,0.04748,-0.000954],
[6.204, 6.979,1.795E4,586.3,0.002776,0.04334,1302,-19.69,6.898,-0.8702,0.04749,-0.0009536],
[6.181, 6.954,1.812E4,586.3,0.002748,0.04385,1289,-19.76,6.912,-0.871,0.04749,-0.0009532],
[6.949, 7.82,1.83E4,586.3,0.002737,0.04436,1284,-19.83,6.926,-0.8718,0.0475,-0.0009528],
[7.506, 8.448,1.848E4,586.3,0.002727,0.04487,1279,-19.9,6.94,-0.8726,0.04751,-0.0009524],
[7.649, 8.609,1.866E4,586.3,0.002697,0.04538,1265,-19.97,6.953,-0.8733,0.04751,-0.000952],
[7.71, 8.679,1.883E4,586.3,0.002641,0.04589,1239,-20.04,6.967,-0.8741,0.04752,-0.0009516],
[7.407, 8.336,1.901E4,586.3,0.002603,0.0464,1221,-20.11,6.981,-0.8749,0.04752,-0.0009512],
[7.29, 8.204,1.918E4,586.3,0.002573,0.04691,1207,-20.18,6.995,-0.8757,0.04753,-0.0009508]
])

158
spanc/Fitter.py
View File

@ -1,158 +0,0 @@
#!/usr/bin/env python3
import numpy as np
from scipy.odr import RealData, ODR, polynomial
class Fit:
def __init__(self, order):
self.poly_order = order
self.model = polynomial(order)
self.x_data = None
self.y_data = None
self.x_errors = None
self.y_errors = None
self.fit_data = None
self.fitter = None
self.output = None
self.parameters = None
def __getstate__(self):
return self.x_data, self.y_data, self.x_errors, self.y_errors, self.parameters, self.poly_order
def __setstate__(self, state):
self.x_data, self.y_data, self.x_errors, self.y_errors, self.parameters, self.poly_order = state
self.model = polynomial(self.poly_order)
def RunFit(self, xarray, yarray, y_errors, x_errors):
self.x_data = xarray
self.y_data = yarray
self.x_errors = x_errors
self.y_errors = y_errors
self.fit_data = RealData(self.x_data, y=self.y_data, sx=self.y_errors, sy=self.x_errors)
self.fitter = ODR(self.fit_data, self.model)
self.output = self.fitter.run()
self.parameters = self.output.beta
def GetParameterError(self, par_index):
return self.output.sd_beta[par_index]
def GetNDF(self):
return len(self.x_data) - self.poly_order+1
class LinearFit(Fit):
def __init__(self):
super().__init__(1)
def EvaluateFunction(self, x):
return self.parameters[0] + self.parameters[1]*x
def EvaluateFunctionDeriv(self, x):
return self.parameters[1]
def EvaluateFunctionParamDeriv(self, x, par_index):
if par_index == 0:
return 1.0
elif par_index == 1:
return x
else:
return 0.0
def ReducedChiSquare(self):
ndf = len(self.x_data)-len(self.parameters)
y_eff_errors = np.zeros(len(self.x_data))
for i in range(len(self.x_data)):
y_eff_errors[i] = np.sqrt(self.y_errors[i]**2.0 + (self.x_errors[i]*self.EvaluateFunctionDeriv(self.x_data[i]))**2.0)
chisq = np.sum(((self.y_data - self.EvaluateFunction(self.x_data))/y_eff_errors)**2.0)
if ndf > 0:
return chisq/ndf
else:
return 0
class QuadraticFit(Fit):
def __init__(self):
super().__init__(2)
def EvaluateFunction(self, x):
return self.parameters[0] + self.parameters[1]*x + self.parameters[2]*x**2.0
def EvaluateFunctionDeriv(self, x):
return self.parameters[1] + 2.0*self.parameters[2]*x
def EvaluateFunctionParamDeriv(self, x, par_index):
if par_index == 0:
return 1.0
elif par_index == 1:
return x
elif par_index == 2:
return x**2.0
else:
return 0.0
def ReducedChiSquare(self):
ndf = len(self.x_data) - len(self.parameters)
y_eff_errors = np.zeros(len(self.x_data))
for i in range(len(self.x_data)):
y_eff_errors[i] = np.sqrt(self.y_errors[i]**2.0 + (self.x_errors[i]*self.EvaluateFunctionDeriv(self.x_data[i]))**2.0)
chisq = np.sum(((self.y_data - self.EvaluateFunction(self.x_data))/y_eff_errors)**2.0)
if ndf >= 0:
return chisq/ndf
else:
return 0
class CubicFit(Fit):
def __init__(self):
super().__init__(3)
def EvaluateFunction(self, x):
return self.parameters[0] + self.parameters[1]*x + self.parameters[2]*x**2.0 + self.parameters[3]*x**3.0
def EvaluateFunctionDeriv(self, x):
return self.parameters[1] + 2.0*self.parameters[2]*x + 3.0*self.parameters[3]*x**2.0
def EvaluateFunctionParamDeriv(self, x, par_index):
if par_index == 0:
return 1.0
elif par_index == 1:
return x
elif par_index == 2:
return x**2.0
elif par_index == 3:
return x**3.0
else:
return 0.0
def ReducedChiSquare(self):
ndf = len(self.x_data) - len(self.parameters)
y_eff_errors = np.zeros(len(self.x_data))
for i in range(len(self.x_data)):
y_eff_errors[i] = np.sqrt(self.y_errors[i]**2.0 + (self.x_errors[i]*self.EvaluateFunctionDeriv(self.x_data[i]))**2.0)
chisq = np.sum(((self.y_data - self.EvaluateFunction(self.x_data))/y_eff_errors)**2.0)
if ndf >= 0:
return chisq/ndf
else:
return 0
def main():
ndata = 100
x = np.zeros(ndata)
y = np.zeros(ndata)
dy = np.zeros(ndata)
dx = np.zeros(ndata)
for i in range(ndata):
x[i] = i
y[i] = i+0.0001*i*i+7.0
dy[i] = 0.1
dx[i] = 0.1
my_fit = LinearFit()
print("Testing SPANC fitting routine using test data and linear function...")
my_fit.RunFit(x, y, dy, dx)
print("Results from fit with y-errors: param[0] =",my_fit.parameters[0],"param[1] =",my_fit.parameters[1],"Reduced chi-sq =",my_fit.ReducedChiSquare())
if __name__ == '__main__':
main()

335
spanc/IngoFit.py
View File

@ -1,335 +0,0 @@
#!/usr/bin/env python3
import numpy as np
class FitFunction :
def __init__(self, nparams, numericFlag=True):
self.nparams = nparams
self.numericFlag = numericFlag
def Evaluate(self, x, params):
return None
def EvaluateParamDeriv(self, x, params, this_param):
return None
def EvaluateDeriv(self, x, params):
return None
def NumericDeriv(self, x, y, dx, params):
return (self.Evaluate(x+dx, params) - y)/dx
def NumericParamDeriv(self, x, y, params, this_param, dpar):
par_val = params[this_param]
params[this_param] += dpar
value = (self.Evaluate(x, params) - y)/dpar
params[this_param] = par_val
return value
class SpancFit:
PRECISION = 1e-6
MAX_ITERS = 100
def __init__(self, func) :
self.function = func
self.nparams = func.nparams
self.numericFlag = func.numericFlag
self.covMatrix = np.zeros((self.nparams, self.nparams))
self.xErrorFlag = False
self.chiSq = 0
self.ndf = 0
self.redChiSq = 0
self.Lambda = 0.1
self.parameters = np.zeros(self.nparams)
self.x_data = np.zeros(0)
self.y_data = np.zeros(0)
self.y_errors = np.zeros(0)
self.x_errors = np.zeros(0)
self.x_incr = np.zeros(0)
self.param_incr = np.zeros(self.nparams)
self.ndata = 0
def EvaluateFunction(self, x):
return self.function.Evaluate(x, self.parameters)
def SetParamIncrement(self, func_vals):
vor = 0.0
nach = 0.0
zoom = 0
index = 0
for i in range(len(self.parameters)):
zoom = 0
self.param_incr[i] = abs(self.parameters[i]*1e-3)
if self.param_incr[i] == 0.0:
self.param_incr[i] = 1e-3
for j in range(5):
index = j*self.ndata/5
nach = self.function.NumericParamDeriv(self.x_data[index], func_vals[index], self.parameters, i, self.param_incr[i])
vor = 0
while zoom <= 4 and abs(nach - vor) >= vor:
zoom += 1
vor = nach
self.param_incr[i] *= 0.1
nach = self.function.NumericParamDeriv(self.x_data[index], func_vals[index], self.parameters, i, self.param_incr[i])
self.param_incr[i]*10.0
def SetXIncrement(self, func_vals):
vor = 0.0
nach = 0.0
zoom = 0
for i in range(self.ndata):
zoom = 0
self.x_incr[i] = abs(self.x_data[i]*1e-3)
if self.x_incr[i] == 0.0:
self.x_incr[i] = 1e-3
nach = self.function.NumericDeriv(self.x_data[i], func_vals[i], self.x_incr[i], self.parameters)
vor = 0
while zoom <= 4 and abs(nach - vor) >= vor:
zoom += 1
vor = nach
self.x_incr[i] *= 0.1
nach = self.function.NumericDeriv(self.x_data[index], func_vals[index], self.x_incr[i], self.parameters)
self.x_incr[i]*10.0
def Fit(self, x_array, y_array, yerr_array, xerr_array=np.empty(0)):
self.x_data = x_array
self.y_data = y_array
self.y_errors = yerr_array
self.x_errors = xerr_array
self.ndata = len(self.x_data)
self.x_incr = np.zeros(self.ndata)
self.ndf = self.ndata - self.nparams
if len(xerr_array) == 0:
self.xErrorFlag = False
else:
self.xErrorFlag = True
self.CurveFit()
def GetParameterError(self, param_index):
return np.sqrt(abs(self.covMatrix[param_index][param_index]))
def CalculateChiSquare(self, func_vals):
y_eff_errors = np.zeros(self.ndata)
if self.numericFlag:
if self.xErrorFlag:
for i in range(self.ndata):
y_eff_errors[i] = np.sqrt(self.y_errors[i]**2.0 + (self.x_errors[i]*self.function.NumericDeriv(self.x_data[i], func_vals[i], self.x_incr[i], self.parameters))**2.0)
else:
y_eff_errors = self.y_errors
elif self.xErrorFlag:
for i in range(self.ndata):
y_eff_errors[i] = np.sqrt(self.y_errors[i]**2.0 + (self.x_errors[i]*self.function.EvaluateDeriv(self.x_data[i], self.parameters))**2.0)
else:
y_eff_errors = self.y_errors
chisq = np.sum(((self.y_data - func_vals)/y_eff_errors)**2.0)
return chisq, y_eff_errors
def CurveFit(self):
hauptschritt=0
sub_iters=0
residuum = np.zeros(self.ndata)
y = 0.0
dy = 0.0
chi_lastmain = 0.0
chi_lastsub = 0.0
chi_fit = 0.0
schlechter = True
b = np.zeros(self.nparams)
abl = np.zeros(self.nparams)
norm = np.zeros(self.nparams)
func_vals = np.zeros(self.ndata)
y_eff_errors = np.zeros(self.ndata)
a = np.zeros((self.nparams, self.nparams))
afaktor = np.zeros((self.nparams, self.nparams))
inv = np.zeros((self.nparams, self.nparams))
for i in range(self.ndata):
func_vals[i] = self.function.Evaluate(self.x_data[i], self.parameters)
if self.numericFlag:
self.SetParamIncrement(func_vals)
if self.xErrorFlag:
self.SetXIncrement(func_vals)
chi_lastsub, y_eff_errors = self.CalculateChiSquare(func_vals)
chi_lastmain = chi_lastsub
temp_params = np.zeros(self.nparams)
#Main optimization loop
while True:
chi_fit = chi_lastmain
for i in range(self.nparams):
for j in range(self.nparams):
a[i][j] = 0.0
b[i] = 0.0
residuum = (self.y_data - func_vals)/(y_eff_errors**2.0)
for i in range(self.ndata):
for j in range(self.nparams):
if(self.numericFlag):
abl[j] = self.function.NumericParamDeriv(self.x_data[i], func_vals[i], self.parameters, j, param_incr[j])
else:
abl[j] = self.function.EvaluateParamDeriv(self.x_data[i], self.parameters, j)
b[j] += abl[j]*residuum[i]
for j in range(self.nparams):
for k in range(self.nparams):
a[j][k] += abl[j]*abl[k]/(y_eff_errors[i]**2.0)
for i in range(self.nparams):
if a[i][i] < 1e-15:
a[i][i] = 1e-15
norm[i] = np.sqrt(a[i][i])
temp_params = self.parameters
#sub-loop looking for the best next step
sub_iters=0
while True:
chi_lastsub = 0.0
for i in range(self.nparams):
for j in range(self.nparams):
afaktor[i][j] = a[i][j]/(norm[i]*norm[j])
afaktor[i][i] = 1.0 + self.Lambda
inv = np.linalg.inv(afaktor)
for i in range(self.nparams):
for j in range(self.nparams):
self.parameters[i] += b[j]*inv[i][j]/(norm[i]*norm[j])
for i in range(self.ndata):
func_vals[i] = self.function.Evaluate(self.x_data[i], self.parameters)
chi_lastsub, y_eff_errors = self.CalculateChiSquare(func_vals)
schlechter = (chi_lastsub - chi_lastmain > 1e-5) and self.Lambda != 0
sub_iters += 1
if sub_iters > self.MAX_ITERS:
break
elif schlechter:
self.parameters = temp_params
self.Lambda *= 10.0
else:
break
#end sub-loop
chi_lastmain = chi_lastsub
if self.numericFlag:
self.SetParamIncrement(func_vals)
self.Lambda *= 0.1
hauptschritt += 1
if abs(chi_fit - chi_lastmain) < self.PRECISION*chi_fit or hauptschritt > self.MAX_ITERS or self.Lambda == 0.0:
break
#end main loop
for i in range(self.nparams):
for j in range(self.nparams):
afaktor[i][j] = a[i][j]/(norm[i]*norm[j])
self.covMatrix = np.linalg.inv(afaktor)
self.chiSq = chi_lastmain
if self.ndata > self.nparams:
self.redChiSq = self.chiSq/self.ndf
else:
self.redChiSq = 0.0
return hauptschritt
class LinearFunction(FitFunction):
def __init__(self):
super().__init__(2, numericFlag=False)
def Evaluate(self, x, params):
return params[0] + x*params[1]
def EvaluateDeriv(self, x, params):
return params[1]
def EvaluateParamDeriv(self, x, params, this_param):
if this_param == 0:
return 1.0
elif this_param == 1:
return x
else:
return 0.0
class QuadraticFunction(FitFunction):
def __init__(self):
super().__init__(3, numericFlag=False)
def Evaluate(self, x, params):
return params[0] + x*params[1] + (x**2.0)*params[2]
def EvaluateDeriv(self, x, params):
return params[1]+2.0*x*params[2]
def EvaluateParamDeriv(self, x, params, this_param):
if this_param == 0:
return 1.0
elif this_param == 1:
return x
elif this_param == 2:
return x**2.0
else:
return 0.0
class CubicFunction(FitFunction):
def __init__(self):
super().__init__(4, numericFlag=False)
def Evaluate(self, x, params):
return params[0] + x*params[1] + (x**2.0)*params[2] + (x**3.0)*params[3]
def EvaluateDeriv(self, x, params):
return params[1]+2.0*x*params[2]+3.0*(x**2.0)*params[3]
def EvaluateParamDeriv(self, x, params, this_param):
if this_param == 0:
return 1.0
elif this_param == 1:
return x
elif this_param == 2:
return x**2.0
elif this_param == 3:
return x**3.0
else:
return 0.0
def main():
ndata = 100
x = np.zeros(ndata)
y = np.zeros(ndata)
dy = np.zeros(ndata)
dx = np.zeros(ndata)
for i in range(ndata):
x[i] = i
y[i] = i+0.0001*i*i+7.0
dy[i] = 0.1
dx[i] = 0.1
fit_func = LinearFunction()
my_fit = SpancFit(fit_func)
my_fit.parameters[0] = 1.0
my_fit.parameters[1] = 1.0
print("Testing SPANC fitting routine using test data and linear function...")
my_fit.Fit(x, y, dy)
print("Results from fit with y-errors: param[0] =",my_fit.parameters[0],"param[1] =",my_fit.parameters[1],"Reduced chi-sq =",my_fit.redChiSq)
if __name__ == '__main__':
main()

133
spanc/LayeredTarget.py
View File

@ -1,133 +0,0 @@
#!/bin/usr/env python3
import numpy as np
import EnergyLoss as eloss
class Target:
def __init__(self):
self.thickness=0
self.Z = np.empty(0)
self.A = np.empty(0)
self.S = np.empty(0)
self.energy_loss = eloss.EnergyLoss()
def SetElements(self, z, a, s, thick):
self.Z = z
self.A = a
self.S = s
self.thickness = thick
self.energy_loss.SetTargetData(self.Z, self.A, self.S)
def GetComposition(self):
comp_string = "("
for i in range(len(self.Z)):
comp_string += "("+ str(self.Z[i]) +","+ str(self.A[i]) +","+ str(self.S[i]) +")"
if not(i == len(self.Z)-1):
comp_string += ","
comp_string += ")"
return comp_string
def HasElement(self, z, a):
for i in range(len(self.Z)):
if self.Z[i] == z and self.A[i] == a:
return True
return False
def GetEnergyLossTotal(self, zp, ap, start_energy, theta):
if theta == np.pi/2.0:
return start_energy
elif theta > np.pi/2.0:
return self.energy_loss.GetEnergyLoss(zp, ap, start_energy, self.thickness/abs(np.cos(np.pi - theta)))
else:
return self.energy_loss.GetEnergyLoss(zp, ap, start_energy, self.thickness/abs(np.cos(theta)))
def GetEnergyLossHalf(self, zp, ap, start_energy, theta):
if theta == np.pi/2.0:
return start_energy
elif theta > np.pi/2.0:
return self.energy_loss.GetEnergyLoss(zp, ap, start_energy, self.thickness/abs(2.0*np.cos(np.pi - theta)))
else:
return self.energy_loss.GetEnergyLoss(zp, ap, start_energy, self.thickness/abs(2.0*np.cos(theta)))
def GetReverseEnergyLossTotal(self, zp, ap, final_energy, theta):
if theta == np.pi/2.0:
return final_energy
elif theta > np.pi/2.0:
return self.energy_loss.GetReverseEnergyLoss(zp, ap, final_energy, self.thickness/abs(np.cos(np.pi - theta)))
else:
return self.energy_loss.GetReverseEnergyLoss(zp, ap, final_energy, self.thickness/abs(np.cos(theta)))
def GetReverseEnergyLossHalf(self, zp, ap, final_energy, theta):
if theta == np.pi/2.0:
return final_energy
elif theta > np.pi/2.0:
return self.energy_loss.GetReverseEnergyLoss(zp, ap, final_energy, self.thickness/abs(2.0*np.cos(np.pi - theta)))
else:
return self.energy_loss.GetReverseEnergyLoss(zp, ap, final_energy, self.thickness/abs(2.0*np.cos(theta)))
class LayeredTarget:
def __init__(self):
self.targets = [] #Order of layers matters!
self.name = ''
def AddLayer(self, z, a, s, thick):
targ = Target()
targ.SetElements(z, a, s, thick)
self.targets.append(targ)
def AddLayer(self, targ):
if not isinstance(targ, Target) :
print("Cannot add layer if it is not of type Target!")
return
else :
self.targets.append(targ)
def FindLayerContainingElement(self, z, a):
for i in range(len(self.targets)):
if self.targets[i].HasElement(z, a):
return i
return -1
def GetEnergyLoss(self, zp, ap, initial_energy, theta, rxn_layer, kind="projectile"):
if rxn_layer < 0 or rxn_layer > len(self.targets):
print("Bad reaction layer at LayeredTarget::GetEnergyLoss")
return 0.0
e_lost = 0.0
new_energy = initial_energy
if kind == "projectile":
for i in range(rxn_layer+1):
if i == rxn_layer:
e_lost += self.targets[i].GetEnergyLossHalf(zp, ap, new_energy, theta)
new_energy = initial_energy - e_lost
else:
e_lost += self.targets[i].GetEnergyLossTotal(zp, ap, new_energy, theta)
new_energy = initial_energy - e_lost
elif kind == "ejectile":
for i in range(rxn_layer, len(self.targets)):
if i == rxn_layer:
e_lost += self.targets[i].GetEnergyLossHalf(zp, ap, new_energy, theta)
new_energy = initial_energy - e_lost
else:
e_lost = self.targets[i].GetEnergyLossTotal(zp, ap, new_energy, theta)
new_energy = initial_energy - e_lost
else:
print("Invalid kind at LayeredTarget::GetEnergyLoss")
return e_lost
def GetReverseEnergyLoss(self, zp, ap, final_energy, theta, rxn_layer):
if rxn_layer < 0 or rxn_layer > len(self.targets):
print("Bad reaction layer at LayeredTarget::GetReverseEnergyLoss")
return 0.0
e_lost = 0.0
new_energy = final_energy
for i in reversed(range(rxn_layer, len(self.targets))):
if i == rxn_layer:
e_lost += self.targets[i].GetReverseEnergyLossHalf(zp, ap, new_energy, theta)
new_energy = final_energy + e_lost
else:
e_lost += self.targets[i].GetReverseEnergyLossTotal(zp, ap, new_energy, theta)
new_energy = final_energy + e_lost
return e_lost

70
spanc/NucData.py
View File

@ -1,70 +0,0 @@
#!/usr/bin/env python3
import numpy as np
import requests
import lxml.html as xhtml
class MassTable:
def __init__(self):
file = open("./etc/mass.txt","r")
self.mtable = {}
u2mev = 931.4940954
me = 0.000548579909
self.etable = {}
file.readline()
file.readline()
for line in file:
entries = line.split()
n = entries[0]
z = entries[1]
a = entries[2]
element = entries[3]
massBig = float(entries[4])
massSmall = float(entries[5])
key = '('+z+','+a+')'
value = ((massBig+massSmall*1e-6) - float(z)*me)*u2mev
self.mtable[key] = value
self.etable[key] = element
file.close()
def GetMass(self, z, a):
key = '('+str(z)+','+str(a)+')'
if key in self.mtable:
return self.mtable[key]
else:
return 0
def GetSymbol(self, z, a):
key = '('+str(z)+','+str(a)+')'
if key in self.etable:
return str(a)+self.etable[key]
else:
return 'none'
Masses = MassTable()
def GetExcitations(symbol):
levels = np.array(np.empty(0))
text = ''
site = requests.get("https://www.nndc.bnl.gov/nudat2/getdatasetClassic.jsp?nucleus="+symbol+"&unc=nds")
contents = xhtml.fromstring(site.content)
tables = contents.xpath("//table")
rows = tables[2].xpath("./tr")
for row in rows[1:-2]:
entries = row.xpath("./td")
if len(entries) != 0:
entry = entries[0]
data = entry.xpath("./a")
if len(data) == 0:
text = entry.text
else:
text = data[0].text
text = text.replace('?', '')
text = text.replace('\xa0\xa0≈','')
levels = np.append(levels, float(text)/1000.0)
return levels

90
spanc/Reaction.py
View File

@ -1,90 +0,0 @@
#!/usr/bin/env python3
from LayeredTarget import LayeredTarget, Target
from NucData import Masses
import numpy as np
class Nucleus:
def __init__(self, z, a):
self.Z = z
self.A = a
self.Symbol = Masses.GetSymbol(self.Z, self.A)
self.GSMass = Masses.GetMass(self.Z, self.A)
def Minus(self, rhs):
final_Z = self.Z - rhs.Z
final_A = self.A - rhs.A
if final_A < 0 or final_Z < 0:
print("Illegal minus operation on Nuclei!")
return Nucleus(0,0)
else:
return Nucleus(final_Z, final_A)
def Plus(self, rhs):
return Nucleus(self.Z + rhs.Z, self.A + rhs.A)
class Reaction:
DEG2RAD = np.pi/180.0 #degrees to radians
C = 299792458 #speed of light m/s
QBRHO2P = 1.0E-9*C #Converts qbrho to p (kG*cm -> MeV/c)
def __init__(self, zt, at, zp, ap, ze, ae, beamKE, theta, bfield, tdata):
self.Target = Nucleus(zt, at)
self.Projectile = Nucleus(zp, ap)
self.Ejectile = Nucleus(ze, ae)
self.Residual = (self.Target.Plus(self.Projectile)).Minus(self.Ejectile)
self.BKE = beamKE
self.Theta = theta * self.DEG2RAD
self.Bfield = bfield
self.Name = self.Target.Symbol +"("+ self.Projectile.Symbol +","+ self.Ejectile.Symbol +")"+ self.Residual.Symbol
self.target_data = tdata
self.rxn_layer = self.target_data.FindLayerContainingElement(self.Target.Z, self.Target.A)
def GetBKEAtRxn(self):
return self.BKE - self.target_data.GetEnergyLoss(self.Projectile.Z, self.Projectile.A, self.BKE, self.Theta, self.rxn_layer)
def GetEjectileKineticEnergyAtRxn(self, Elevel) :
Q = self.Target.GSMass + self.Projectile.GSMass - (self.Ejectile.GSMass + self.Residual.GSMass + Elevel)
Ethresh = -Q*(self.Ejectile.GSMass+self.Residual.GSMass)/(self.Ejectile.GSMass + self.Residual.GSMass - self.Projectile.GSMass)
BKE_rxn = self.GetBKEAtRxn()
if BKE_rxn < Ethresh:
return 0.0
term1 = np.sqrt(self.Projectile.GSMass*self.Ejectile.GSMass*BKE_rxn)/(self.Ejectile.GSMass + self.Residual.GSMass)*np.cos(self.Theta)
term2 = (BKE_rxn*(self.Residual.GSMass - self.Projectile.GSMass) + self.Residual.GSMass*Q)/(self.Ejectile.GSMass + self.Residual.GSMass)
ke1 = term1 + np.sqrt(term1**2.0 + term2)
ke2 = term1 - np.sqrt(term1**2.0 + term2)
if ke1 > 0:
return ke1**2.0
else :
return ke2**2.0
def GetEjectileKineticEnergyAtDet(self, Elevel):
KE_at_rxn = self.GetEjectileKineticEnergyAtRxn(Elevel)
KE_at_det = KE_at_rxn - self.target_data.GetEnergyLoss(self.Ejectile.Z, self.Ejectile.A, KE_at_rxn, self.Theta, self.rxn_layer, kind="ejectile")
return KE_at_det
def GetEjectileRho(self, Elevel):
KE_at_det = self.GetEjectileKineticEnergyAtDet(Elevel)
p = np.sqrt(KE_at_det*(KE_at_det + 2.0*self.Ejectile.GSMass))
qbrho = p/self.QBRHO2P
return qbrho/(self.Ejectile.Z*self.Bfield)
def GetResidualExcitation(self, rho):
p_eject_at_det = rho*self.Ejectile.Z*self.Bfield*self.QBRHO2P
KE_eject_at_det = np.sqrt(p_eject_at_det**2.0 + self.Ejectile.GSMass**2.0) - self.Ejectile.GSMass
KE_eject_at_rxn = KE_eject_at_det + self.target_data.GetReverseEnergyLoss(self.Ejectile.Z, self.Ejectile.A, KE_eject_at_det, self.Theta, self.rxn_layer)
p_eject_at_rxn = np.sqrt(KE_eject_at_rxn*(KE_eject_at_rxn + 2.0*self.Ejectile.GSMass))
E_eject_at_rxn = KE_eject_at_rxn+self.Ejectile.GSMass
BKE_atRxn = self.GetBKEAtRxn()
E_project = BKE_atRxn + self.Projectile.GSMass
p_project = np.sqrt(BKE_atRxn*(BKE_atRxn + 2.0*self.Projectile.GSMass))
E_resid = E_project + self.Target.GSMass - E_eject_at_rxn
p2_resid = p_project**2.0 + p_eject_at_rxn**2.0 - 2.0*p_project*p_eject_at_rxn*np.cos(self.Theta)
m_resid = np.sqrt(E_resid**2.0 - p2_resid)
return m_resid - self.Residual.GSMass
def ChangeReactionParameters(self, bke, theta, bf) :
self.BKE = bke
self.Theta = theta*self.DEG2RAD
self.Bfield = bf

161
spanc/Spanc.py
View File

@ -1,161 +0,0 @@
#!/usr/bin/env python3
from Reaction import Reaction
from LayeredTarget import LayeredTarget, Target
from IngoFit import SpancFit, LinearFunction, QuadraticFunction, CubicFunction
from Fitter import LinearFit, QuadraticFit, CubicFit
import numpy as np
class Peak :
def __init__(self):
self.Ex = 0.0
self.uEx = 0.0
self.x = 0.0
self.ux_sys = 0.0
self.ux_stat = 0.0
self.rho = 0.0
self.urho = 0.0
self.fwhm_x = 0.0
self.ufwhm_x = 0.0
self.fwhm_Ex = 0.0
self.ufwhm_Ex = 0.0
self.reaction = ""
class Spanc:
def __init__(self):
self.reactions = {}
self.targets = {}
self.calib_peaks = {}
self.output_peaks = {}
self.fitters = {}
self.InitFits()
def WriteConfig(self):
return
def ReadConfig(self):
return
def InitFits(self):
"""
self.fitters["linear"] = SpancFit(LinearFunction())
self.fitters["quadratic"] = SpancFit(QuadraticFunction())
self.fitters["cubic"] = SpancFit(CubicFunction())
"""
self.fitters["linear"] = LinearFit()
self.fitters["quadratic"] = QuadraticFit()
self.fitters["cubic"] = CubicFit()
def PerformFits(self):
xarray = np.empty(0)
yarray = np.empty(0)
uxarray = np.empty(0)
uyarray = np.empty(0)
for peak in self.calib_peaks.values():
xarray = np.append(xarray, peak.x)
uxarray = np.append(uxarray, np.sqrt(peak.ux_sys**2.0 + peak.ux_stat**2.0))
yarray = np.append(yarray, peak.rho)
uyarray = np.append(uyarray, peak.urho)
for key in self.fitters.keys():
#self.fitters[key].Fit(xarray, yarray, uyarray, uxarray)
self.fitters[key].RunFit(xarray, yarray, uyarray, uxarray)
return xarray, yarray, uxarray, uyarray
def CalculateResiduals(self, fit_name):
fit = self.fitters[fit_name]
npeaks = len(self.calib_peaks)
resids = np.empty(npeaks)
student_resids = np.empty(npeaks)
xarray = np.empty(npeaks)
counter=0
for peak in self.calib_peaks.values():
fval = fit.EvaluateFunction(peak.x)
dval = peak.rho
resids[counter] = (dval - fval)
xarray[counter] = peak.x
counter += 1
mean_x = np.average(xarray)/npeaks
rmse = np.sqrt(np.sum(resids**2.0)/fit.GetNDF())
#get leverage
counter=0
sq_diff=0
for peak in self.calib_peaks.values():
sq_diff += (peak.x - mean_x)**2.0
sq_diff = sq_diff/npeaks
for peak in self.calib_peaks.values():
leverage = 1.0/npeaks + (peak.x - mean_x)/sq_diff
student_resids[counter] = resids[counter]/(rmse*np.sqrt(1.0 - leverage))
counter += 1
return xarray, resids, student_resids
def AddCalibrationPeak(self, rxn_name, cal_name, position, ux_stat, ux_sys, ex, uex) :
new_peak = Peak()
new_peak.x = position
new_peak.reaction = rxn_name
new_peak.ux_stat = ux_stat
new_peak.ux_sys = ux_sys
new_peak.Ex = ex
new_peak.uEx = uex
new_peak.rho = self.reactions[rxn_name].GetEjectileRho(ex)
new_peak.urho = abs(self.reactions[rxn_name].GetEjectileRho(ex+uex)-new_peak.rho)
self.calib_peaks[cal_name] = new_peak
def AddOutputPeak(self, rxn_name, out_name, position, ux_stat, ux_sys, fwhm_x, ufwhm_x) :
new_peak = Peak()
new_peak.x = position
new_peak.ux_stat = ux_stat
new_peak.ux_sys = ux_sys
new_peak.reaction = rxn_name
new_peak.fwhm_x = fwhm_x
new_peak.ufwhm_x = ufwhm_x
self.output_peaks[out_name] = new_peak
def CalculateRhoUncertainty(self, peak, fit, deltax=0.0, udeltax=0.0):
urho = 0
for i in range(len(fit.parameters)):
urho += (fit.EvaluateFunctionParamDeriv(peak.x+deltax, i)*fit.GetParameterError(i))**2.0
urho += (fit.EvaluateFunctionDeriv(peak.x+deltax)*np.sqrt(peak.ux_stat**2.0 + peak.ux_sys**2.0 + udeltax**2.0))**2.0
urho = np.sqrt(urho)
return urho
def CalculateOutputs(self, fit_name):
fit = self.fitters[fit_name]
for output in self.output_peaks.values():
output.rho = fit.EvaluateFunction(output.x)
output.urho = self.CalculateRhoUncertainty(output, fit)
output.Ex = self.reactions[output.reaction].GetResidualExcitation(output.rho)
output.uEx = abs(self.reactions[output.reaction].GetResidualExcitation(output.rho + output.urho) - output.Ex)
if output.fwhm_x == 0:
output.fwhm_Ex = 0
output.ufwhm_Ex = 0
else:
rhoLo = fit.EvaluateFunction(output.x - output.fwhm_x/2.0)
urhoLo = self.CalculateRhoUncertainty(output, fit, deltax=-1.0*output.fwhm_x/2.0, udeltax=output.ufwhm_x/2.0)
rhoHi = fit.EvaluateFunction(output.x + output.fwhm_x/2.0)
urhoHi = self.CalculateRhoUncertainty(output, fit, deltax=output.fwhm_x/2.0, udeltax=output.ufwhm_x/2.0)
exLo = self.reactions[output.reaction].GetResidualExcitation(rhoLo)
uexLo = abs(self.reactions[output.reaction].GetResidualExcitation(rhoLo+urhoLo) - exLo)
exHi = self.reactions[output.reaction].GetResidualExcitation(rhoHi)
uexHi = abs(self.reactions[output.reaction].GetResidualExcitation(rhoHi+urhoHi) - exHi)
output.fwhm_Ex = abs(exHi - exLo)
output.ufwhm_Ex = output.ufwhm_x/output.fwhm_x*output.fwhm_Ex
def CalculateCalibrations(self):
for peak in self.calib_peaks.values():
peak.rho = self.reactions[peak.reaction].GetEjectileRho(peak.Ex)
peak.urho = abs(self.reactions[peak.reaction].GetEjectileRho(peak.Ex+peak.uEx) - peak.rho)

900
spanc/SpancGUI.py
View File

@ -1,900 +0,0 @@
#!/usr/bin/env python3
import Spanc as spc
from LayeredTarget import LayeredTarget, Target
from Reaction import Reaction
import sys
from qtpy.QtWidgets import QApplication, QWidget, QMainWindow
from qtpy.QtWidgets import QLabel, QMenuBar, QAction
from qtpy.QtWidgets import QHBoxLayout, QVBoxLayout, QGridLayout, QGroupBox
from qtpy.QtWidgets import QPushButton, QButtonGroup, QRadioButton
from qtpy.QtWidgets import QSpinBox, QDoubleSpinBox, QComboBox
from qtpy.QtWidgets import QDialog, QFileDialog, QDialogButtonBox
from qtpy.QtWidgets import QTableWidget, QTableWidgetItem
from qtpy.QtWidgets import QLineEdit, QTabWidget, QFormLayout
from qtpy.QtCore import Signal
import matplotlib as mpl
import pickle as pickle
import numpy as np
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg
from matplotlib.figure import Figure
class MPLCanvas(FigureCanvasQTAgg):
def __init__(self, parent=None, width=5, height=4, dpi=100):
self.fig = Figure(figsize=(width, height), dpi=dpi, edgecolor="black",linewidth=0.5,constrained_layout=True)
self.axes = self.fig.add_subplot(111)
self.axes.spines['top'].set_visible(False)
super(MPLCanvas, self).__init__(self.fig)
class TargetDialog(QDialog):
new_target = Signal(list, str)
def __init__(self, parent=None, target=None):
super().__init__(parent)
self.setWindowTitle("Add A Target")
nameLabel = QLabel("Target name", self)
self.nameInput = QLineEdit(self)
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
self.buttonBox.accepted.connect(self.SendTarget)
self.buttonBox.rejected.connect(self.reject)
self.layout = QVBoxLayout()
self.setLayout(self.layout)
self.layout.addWidget(nameLabel)
self.layout.addWidget(self.nameInput)
self.CreateTargetInputs()
if target is not None:
self.SetInitialValues(target)
self.layout.addWidget(self.buttonBox)
def CreateTargetInputs(self):
self.layerAInputs = []
self.layerZInputs = []
self.layerSInputs = []
self.layerThickInputs = []
self.layer1GroupBox = QGroupBox("Layer 1", self)
layer1Layout = QVBoxLayout()
thick1Label = QLabel("Thickness(ug/cm^2)", self.layer1GroupBox)
self.layerThickInputs.append(QDoubleSpinBox(self.layer1GroupBox))
self.layerThickInputs[0].setRange(0, 999.0)
self.layerThickInputs[0].setDecimals(4)
self.layer1ComponentsBox = QGroupBox("Layer 1 Components", self.layer1GroupBox)
layer1compLayout = QGridLayout()
layer1compLayout.addWidget(QLabel("Z", self.layer1ComponentsBox), 0, 1)
layer1compLayout.addWidget(QLabel("A", self.layer1ComponentsBox), 0, 2)
layer1compLayout.addWidget(QLabel("Stoich", self.layer1ComponentsBox), 0, 3)
for i in range(3):
layer1compLayout.addWidget(QLabel("Component"+str(i), self.layer1ComponentsBox), i+1, 0)
self.layerZInputs.append(QSpinBox(self.layer1ComponentsBox))
self.layerAInputs.append(QSpinBox(self.layer1ComponentsBox))
self.layerSInputs.append(QSpinBox(self.layer1ComponentsBox))
layer1compLayout.addWidget(self.layerZInputs[i], i+1, 1)
layer1compLayout.addWidget(self.layerAInputs[i], i+1, 2)
layer1compLayout.addWidget(self.layerSInputs[i], i+1, 3)
self.layer1ComponentsBox.setLayout(layer1compLayout)
layer1Layout.addWidget(thick1Label)
layer1Layout.addWidget(self.layerThickInputs[0])
layer1Layout.addWidget(self.layer1ComponentsBox)
self.layer1GroupBox.setLayout(layer1Layout)
self.layer2GroupBox = QGroupBox("Layer 2", self)
layer2Layout = QVBoxLayout()
thick2Label = QLabel("Thickness(ug/cm^2)", self.layer2GroupBox)
self.layerThickInputs.append(QDoubleSpinBox(self.layer2GroupBox))
self.layerThickInputs[1].setRange(0, 999.0)
self.layerThickInputs[1].setDecimals(4)
self.layer2ComponentsBox = QGroupBox("Layer 2 Components", self.layer2GroupBox)
layer2compLayout = QGridLayout()
layer2compLayout.addWidget(QLabel("Z", self.layer2ComponentsBox), 0, 1)
layer2compLayout.addWidget(QLabel("A", self.layer2ComponentsBox), 0, 2)
layer2compLayout.addWidget(QLabel("Stoich", self.layer2ComponentsBox), 0, 3)
for i in range(3):
layer2compLayout.addWidget(QLabel("Component"+str(i), self.layer2ComponentsBox), i+1, 0)
self.layerZInputs.append(QSpinBox(self.layer2ComponentsBox))
self.layerAInputs.append(QSpinBox(self.layer2ComponentsBox))
self.layerSInputs.append(QSpinBox(self.layer2ComponentsBox))
layer2compLayout.addWidget(self.layerZInputs[i+3], i+1, 1)
layer2compLayout.addWidget(self.layerAInputs[i+3], i+1, 2)
layer2compLayout.addWidget(self.layerSInputs[i+3], i+1, 3)
self.layer2ComponentsBox.setLayout(layer2compLayout)
layer2Layout.addWidget(thick2Label)
layer2Layout.addWidget(self.layerThickInputs[1])
layer2Layout.addWidget(self.layer2ComponentsBox)
self.layer2GroupBox.setLayout(layer2Layout)
self.layer3GroupBox = QGroupBox("Layer 3", self)
layer3Layout = QVBoxLayout()
thick3Label = QLabel("Thickness(ug/cm^2)", self.layer3GroupBox)
self.layerThickInputs.append(QDoubleSpinBox(self.layer3GroupBox))
self.layerThickInputs[2].setRange(0, 999.0)
self.layerThickInputs[2].setDecimals(4)
self.layer3ComponentsBox = QGroupBox("Layer 3 Components", self.layer3GroupBox)
layer3compLayout = QGridLayout()
layer3compLayout.addWidget(QLabel("Z", self.layer3ComponentsBox), 0, 1)
layer3compLayout.addWidget(QLabel("A", self.layer3ComponentsBox), 0, 2)
layer3compLayout.addWidget(QLabel("Stoich", self.layer3ComponentsBox), 0, 3)
for i in range(3):
layer3compLayout.addWidget(QLabel("Component"+str(i), self.layer3ComponentsBox), i+1, 0)
self.layerZInputs.append(QSpinBox(self.layer3ComponentsBox))
self.layerAInputs.append(QSpinBox(self.layer3ComponentsBox))
self.layerSInputs.append(QSpinBox(self.layer3ComponentsBox))
layer3compLayout.addWidget(self.layerZInputs[i+6], i+1, 1)
layer3compLayout.addWidget(self.layerAInputs[i+6], i+1, 2)
layer3compLayout.addWidget(self.layerSInputs[i+6], i+1, 3)
self.layer3ComponentsBox.setLayout(layer3compLayout)
layer3Layout.addWidget(thick3Label)
layer3Layout.addWidget(self.layerThickInputs[2])
layer3Layout.addWidget(self.layer3ComponentsBox)
self.layer3GroupBox.setLayout(layer3Layout)
self.layout.addWidget(self.layer1GroupBox)
self.layout.addWidget(self.layer2GroupBox)
self.layout.addWidget(self.layer3GroupBox)
def SetInitialValues(self, target):
self.nameInput.setText(target.name)
self.nameInput.setReadOnly(True)
for i in range(len(target.targets)):
self.layerThickInputs[i].setValue(target.targets[i].thickness)
for j in range(len(target.targets[i].Z)):
self.layerZInputs[j+i*3].setValue(target.targets[i].Z[j])
self.layerAInputs[j+i*3].setValue(target.targets[i].A[j])
self.layerSInputs[j+i*3].setValue(target.targets[i].S[j])
def SendTarget(self):
name = self.nameInput.text()
if name == "":
return
t1 = Target()
t2 = Target()
t3 = Target()
tlist = []
Z1 = []
A1 = []
S1 = []
Z2 = []
A2 = []
S2 = []
Z3 = []
A3 = []
S3 = []
z = 0
a = 0
s = 0
thick1 = self.layerThickInputs[0].value()
thick2 = self.layerThickInputs[1].value()
thick3 = self.layerThickInputs[2].value()
for i in range(3):
z = self.layerZInputs[i].value()
a = self.layerAInputs[i].value()
s = self.layerSInputs[i].value()
if z != 0 and a != 0 and s != 0:
Z1.append(z)
A1.append(a)
S1.append(s)
z = self.layerZInputs[i+3].value()
a = self.layerAInputs[i+3].value()
s = self.layerSInputs[i+3].value()
if z != 0 and a != 0 and s != 0:
Z2.append(z)
A2.append(a)
S2.append(s)
z = self.layerZInputs[i+6].value()
a = self.layerAInputs[i+6].value()
s = self.layerSInputs[i+6].value()
if z != 0 and a != 0 and s != 0:
Z3.append(z)
A3.append(a)
S3.append(s)
if len(Z1) != 0:
t1.SetElements(Z1, A1, S1, thick1)
tlist.append(t1)
if len(Z2) != 0:
t2.SetElements(Z2, A2, S2, thick2)
tlist.append(t2)
if len(Z3) != 0:
t3.SetElements(Z3, A3, S3, thick3)
tlist.append(t3)
if len(tlist) != 0:
self.new_target.emit(tlist, name)
class ReactionDialog(QDialog):
new_reaction = Signal(int, int, int, int, int, int, float, float, float, str)
update_reaction = Signal(float, float, float, str)
def __init__(self, parent=None, rxn=None, rxnKey=None):
super().__init__(parent)
self.setWindowTitle("Add A Reaction")
tnameLabel = QLabel("Target Name", self)
self.targetNameBox = QComboBox(self)
for target in parent.spanc.targets:
self.targetNameBox.addItem(target)
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
if rxn is not None:
self.buttonBox.accepted.connect(self.SendReactionUpdate)
else:
self.buttonBox.accepted.connect(self.SendReaction)
self.buttonBox.rejected.connect(self.reject)
self.layout = QVBoxLayout()
self.setLayout(self.layout)
self.layout.addWidget(tnameLabel)
self.layout.addWidget(self.targetNameBox)
self.CreateReactionInputs()
if rxn is not None:
self.SetInitialValues(rxn)
self.rxnKey = rxnKey
self.layout.addWidget(self.buttonBox)
def SendReaction(self) :
self.new_reaction.emit(self.ztInput.value(),self.atInput.value(),self.zpInput.value(),self.apInput.value(),self.zeInput.value(),self.aeInput.value(), self.bkeInput.value(), self.thetaInput.value(), self.bfieldInput.value(), self.targetNameBox.currentText())
def SendReactionUpdate(self):
self.update_reaction.emit(self.bkeInput.value(), self.thetaInput.value(), self.bfieldInput.value(), self.rxnKey)
def CreateReactionInputs(self) :
self.nucleiGroupBox = QGroupBox("Reaction Nuclei",self)
inputLayout = QFormLayout()
self.ztInput = QSpinBox(self.nucleiGroupBox)
self.ztInput.setRange(1, 110)
self.atInput = QSpinBox(self.nucleiGroupBox)
self.atInput.setRange(1,270)
self.zpInput = QSpinBox(self.nucleiGroupBox)
self.zpInput.setRange(1, 110)
self.apInput = QSpinBox(self.nucleiGroupBox)
self.apInput.setRange(1,270)
self.zeInput = QSpinBox(self.nucleiGroupBox)
self.zeInput.setRange(1, 110)
self.aeInput = QSpinBox(self.nucleiGroupBox)
self.aeInput.setRange(1,270)
inputLayout.addRow("ZT",self.ztInput)
inputLayout.addRow("AT",self.atInput)
inputLayout.addRow("ZP",self.zpInput)
inputLayout.addRow("AP",self.apInput)
inputLayout.addRow("ZE",self.zeInput)
inputLayout.addRow("AE",self.aeInput)
self.parameterGroupBox = QGroupBox("Reaction Parameters", self)
parameterLayout = QFormLayout()
self.bkeInput = QDoubleSpinBox(self.parameterGroupBox)
self.bkeInput.setRange(0.0, 40.0)
self.bkeInput.setDecimals(4)
self.thetaInput = QDoubleSpinBox(self.parameterGroupBox)
self.thetaInput.setRange(0.0, 180.0)
self.thetaInput.setDecimals(4)
self.bfieldInput = QDoubleSpinBox(self.parameterGroupBox)
self.bfieldInput.setRange(0.0, 16.0)
self.bfieldInput.setDecimals(6)
parameterLayout.addRow("Beam KE(Mev)",self.bkeInput)
parameterLayout.addRow("Theta(deg)",self.thetaInput)
parameterLayout.addRow("Bfield(kG)",self.bfieldInput)
self.nucleiGroupBox.setLayout(inputLayout)
self.parameterGroupBox.setLayout(parameterLayout)
self.layout.addWidget(self.nucleiGroupBox)
self.layout.addWidget(self.parameterGroupBox)
def SetInitialValues(self, rxn):
self.targetNameBox.setCurrentIndex(self.targetNameBox.findText(rxn.target_data.name))
self.targetNameBox.setEnabled(False)
self.ztInput.setValue(rxn.Target.Z)
self.ztInput.setEnabled(False)
self.atInput.setValue(rxn.Target.A)
self.atInput.setEnabled(False)
self.zpInput.setValue(rxn.Projectile.Z)
self.zpInput.setEnabled(False)
self.apInput.setValue(rxn.Projectile.A)
self.apInput.setEnabled(False)
self.zeInput.setValue(rxn.Ejectile.Z)
self.zeInput.setEnabled(False)
self.aeInput.setValue(rxn.Ejectile.A)
self.aeInput.setEnabled(False)
self.bkeInput.setValue(rxn.BKE)
self.thetaInput.setValue(rxn.Theta/rxn.DEG2RAD)
self.bfieldInput.setValue(rxn.Bfield)
class PeakDialog(QDialog):
new_calibration = Signal(float, float, float, float, float, str)
update_calibration = Signal(float, float, float, float, float, str, str)
new_output = Signal(float, float, float, float, float, str)
update_output = Signal(float, float, float, float, float, str, str)
def __init__(self, peakType, parent=None, peak=None, peakKey=None):
super().__init__(parent)
self.setWindowTitle("Add A "+peakType+" Peak")
rnameLabel = QLabel("Reaction Name", self)
self.rxnNameBox = QComboBox(self)
for reaction in parent.spanc.reactions:
self.rxnNameBox.addItem(reaction)
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
if peak is not None and peakType == "Calibration":
self.buttonBox.accepted.connect(self.SendUpdateCalibrationPeak)
elif peakType == "Calibration":
self.buttonBox.accepted.connect(self.SendCalibrationPeak)
elif peak is not None and peakType == "Output":
self.buttonBox.accepted.connect(self.SendUpdateOutputPeak)
elif peakType == "Output":
self.buttonBox.accepted.connect(self.SendOutputPeak)
self.buttonBox.rejected.connect(self.reject)
self.layout = QVBoxLayout()
self.setLayout(self.layout)
self.layout.addWidget(rnameLabel)
self.layout.addWidget(self.rxnNameBox)
if peakType == "Calibration":
self.CreateCalibrationInputs()
if peak is not None:
self.peakKey = peakKey
self.SetCalibrationInputs(peak)
elif peakType == "Output":
self.CreateOutputInputs()
if peak is not None:
self.peakKey = peakKey
self.SetOutputInputs(peak)
self.layout.addWidget(self.buttonBox)
def CreateCalibrationInputs(self):
self.inputGroupBox = QGroupBox("Peak Parameters",self)
inputLayout = QFormLayout()
self.xInput = QDoubleSpinBox(self.inputGroupBox)
self.xInput.setRange(-999, 999)
self.xInput.setDecimals(6)
self.uxsysInput = QDoubleSpinBox(self.inputGroupBox)
self.uxsysInput.setRange(-999, 999)
self.uxsysInput.setDecimals(6)
self.uxstatInput = QDoubleSpinBox(self.inputGroupBox)
self.uxstatInput.setRange(-999, 999)
self.uxstatInput.setDecimals(6)
self.exInput = QDoubleSpinBox(self.inputGroupBox)
self.exInput.setDecimals(6)
self.uexInput = QDoubleSpinBox(self.inputGroupBox)
self.uexInput.setDecimals(6)
inputLayout.addRow("Position(mm)", self.xInput)
inputLayout.addRow("Position Stat. Error(mm)", self.uxstatInput)
inputLayout.addRow("Position Sys. Error(mm)", self.uxsysInput)
inputLayout.addRow("Excitation Energy(MeV)", self.exInput)
inputLayout.addRow("Excitation Energy Error(MeV)", self.uexInput)
self.inputGroupBox.setLayout(inputLayout)
self.layout.addWidget(self.inputGroupBox)
def CreateOutputInputs(self):
self.inputGroupBox = QGroupBox("Peak Parameters",self)
inputLayout = QFormLayout()
self.xInput = QDoubleSpinBox(self.inputGroupBox)
self.xInput.setRange(-999, 999)
self.xInput.setDecimals(6)
self.uxsysInput = QDoubleSpinBox(self.inputGroupBox)
self.uxsysInput.setRange(-999, 999)
self.uxsysInput.setDecimals(6)
self.uxstatInput = QDoubleSpinBox(self.inputGroupBox)
self.uxstatInput.setRange(-999, 999)
self.uxstatInput.setDecimals(6)
self.fwhmInput = QDoubleSpinBox(self.inputGroupBox)
self.fwhmInput.setDecimals(6)
self.ufwhmInput = QDoubleSpinBox(self.inputGroupBox)
self.ufwhmInput.setDecimals(6)
inputLayout.addRow("Position(mm)", self.xInput)
inputLayout.addRow("Position Stat. Error(mm)", self.uxstatInput)
inputLayout.addRow("Position Sys. Error(mm)", self.uxsysInput)
inputLayout.addRow("Position FWHM(mm)", self.fwhmInput)
inputLayout.addRow("Position FWHM Error(mm)", self.ufwhmInput)
self.inputGroupBox.setLayout(inputLayout)
self.layout.addWidget(self.inputGroupBox)
def SetCalibrationInputs(self, peak):
self.rxnNameBox.setCurrentIndex(self.rxnNameBox.findText(peak.reaction))
self.xInput.setValue(peak.x)
self.uxsysInput.setValue(peak.ux_sys)
self.uxstatInput.setValue(peak.ux_stat)
self.exInput.setValue(peak.Ex)
self.uexInput.setValue(peak.uEx)
def SetOutputInputs(self, peak):
self.rxnNameBox.setCurrentIndex(self.rxnNameBox.findText(peak.reaction))
self.xInput.setValue(peak.x)
self.uxsysInput.setValue(peak.ux_sys)
self.uxstatInput.setValue(peak.ux_stat)
self.fwhmInput.setValue(peak.fwhm_x)
self.ufwhmInput.setValue(peak.ufwhm_x)
def SendCalibrationPeak(self):
self.new_calibration.emit(self.xInput.value(), self.uxstatInput.value(), self.uxsysInput.value(), self.exInput.value(), self.uexInput.value(), self.rxnNameBox.currentText())
def SendUpdateCalibrationPeak(self):
self.update_calibration.emit(self.xInput.value(), self.uxstatInput.value(), self.uxsysInput.value(), self.exInput.value(), self.uexInput.value(), self.rxnNameBox.currentText(),self.peakKey)
def SendOutputPeak(self):
self.new_output.emit(self.xInput.value(), self.uxstatInput.value(), self.uxsysInput.value(), self.fwhmInput.value(), self.ufwhmInput.value(), self.rxnNameBox.currentText())
def SendUpdateOutputPeak(self):
self.update_output.emit(self.xInput.value(), self.uxstatInput.value(), self.uxsysInput.value(), self.fwhmInput.value(), self.ufwhmInput.value(), self.rxnNameBox.currentText(),self.peakKey)
class SpancGUI(QMainWindow):
def __init__(self, parent=None):
super().__init__(parent)
self.setWindowTitle("SPANC")
self.spanc = spc.Spanc()
self.tablelayout = QVBoxLayout()
self.plotlayout = QGridLayout() #2x2 grid
self.layout = QVBoxLayout()
self.centralWidget = QTabWidget(self)
self.setCentralWidget(self.centralWidget)
self.centralWidget.setLayout(self.layout)
self.tableTab = QWidget(self.centralWidget)
self.tableTab.setLayout(self.tablelayout)
self.plotTab = QWidget(self.centralWidget)
self.plotTab.setLayout(self.plotlayout)
self.centralWidget.addTab(self.tableTab, "Data Tables")
self.centralWidget.addTab(self.plotTab, "Plots and Fits")
self.fitFlag = False
self.CreateMenus()
self.CreateFitCanvas() #(0,0)
self.CreateResidualCanvas() #(1,0)
self.CreateTargetTable()
self.CreateReactionTable()
self.CreateCalibrationTable()
self.CreateOutputTable()
self.CreateFitTable() #(0,1)
self.CreateResidualTable() #(1,1)
self.show()
def CreateMenus(self):
self.fileMenu = self.menuBar().addMenu("&File")
saveAction = QAction("&Save...",self)
openAction = QAction("&Open...",self)
saveFitAction = QAction("Save Fit Plot...", self)
self.fileMenu.addAction(saveAction)
self.fileMenu.addAction(openAction)
self.fileMenu.addAction(saveFitAction)
self.fileMenu.addAction("&Exit", self.close)
saveAction.triggered.connect(self.HandleSave)
openAction.triggered.connect(self.HandleOpen)
saveFitAction.triggered.connect(self.HandleSaveFit)
self.addMenu = self.menuBar().addMenu("&New")
newTargetAction = QAction("New target...", self)
newReactionAction = QAction("New reaction...", self)
newCalibrationAction = QAction("New calibration...", self)
newOutputAction = QAction("New output...", self)
self.addMenu.addAction(newTargetAction)
self.addMenu.addAction(newReactionAction)
self.addMenu.addAction(newCalibrationAction)
self.addMenu.addAction(newOutputAction)
newTargetAction.triggered.connect(self.HandleNewTarget)
newReactionAction.triggered.connect(self.HandleNewReaction)
newCalibrationAction.triggered.connect(self.HandleNewCalibration)
newOutputAction.triggered.connect(self.HandleNewOutput)
def CreateFitCanvas(self):
self.fitGroup = QGroupBox("Calibration Fit", self.plotTab)
fitLayout = QVBoxLayout()
self.fitCanvas = MPLCanvas(self.fitGroup, width=6, height=6, dpi=100)
self.fitOptionGroup = QGroupBox("Fit options", self.fitGroup)
fitOptionLayout = QHBoxLayout()
self.fitButton = QPushButton("Run Fit", self.fitOptionGroup)
self.fitButton.clicked.connect(self.HandleRunFit)
self.fitTypeBox = QComboBox(self.fitOptionGroup)
self.fitTypeBox.addItem("linear")
self.fitTypeBox.addItem("quadratic")
self.fitTypeBox.addItem("cubic")
fitOptionLayout.addWidget(QLabel("Fit type", self.fitOptionGroup))
fitOptionLayout.addWidget(self.fitTypeBox)
fitOptionLayout.addWidget(self.fitButton)
self.fitOptionGroup.setLayout(fitOptionLayout)
fitLayout.addWidget(self.fitCanvas)
fitLayout.addWidget(self.fitOptionGroup)
self.fitGroup.setLayout(fitLayout)
self.plotlayout.addWidget(self.fitGroup,0,0)
def CreateResidualCanvas(self):
self.residGroup = QGroupBox("Fit Residuals", self.plotTab)
residLayout = QVBoxLayout()
self.residCanvas = MPLCanvas(self.residGroup, width=6, height=6, dpi=100)
self.residOptionGroup = QGroupBox("Fit options", self.residGroup)
residOptionLayout = QHBoxLayout()
self.residButton = QPushButton("Run Resids.", self.residOptionGroup)
self.residButton.clicked.connect(self.HandleRunResids)
residOptionLayout.addWidget(self.residButton)
self.residOptionGroup.setLayout(residOptionLayout)
residLayout.addWidget(self.residCanvas)
residLayout.addWidget(self.residOptionGroup)
self.residGroup.setLayout(residLayout)
self.plotlayout.addWidget(self.residGroup,1,0)
def CreateTargetTable(self):
self.targetGroup = QGroupBox("Targets", self.tableTab)
targetLayout = QVBoxLayout()
self.targetTable = QTableWidget(self.targetGroup)
self.targetTable.setColumnCount(6)
self.targetTable.setHorizontalHeaderLabels(["Layer1 Thickness(ug/cm^2", "Layer1 (Z, A, S)","Layer2 Thickness(ug/cm^2", "Layer2 (Z, A, S)","Layer3 Thickness(ug/cm^2", "Layer3 (Z, A, S)"])
targetLayout.addWidget(self.targetTable)
self.targetGroup.setLayout(targetLayout)
self.tablelayout.addWidget(self.targetGroup)
self.targetTable.resizeColumnsToContents()
self.targetTable.cellDoubleClicked.connect(self.HandleUpdateTarget)
def CreateReactionTable(self):
self.rxnGroup = QGroupBox("Reactions", self.tableTab)
rxnLayout = QVBoxLayout()
self.reactionTable = QTableWidget(self.rxnGroup)
self.reactionTable.setColumnCount(12)
self.reactionTable.setHorizontalHeaderLabels(["Target","ZT","AT","ZP","AP","ZE","AE","ZR","AR","Beam KE(MeV)","BField(kG)","Angle(deg)"])
rxnLayout.addWidget(self.reactionTable)
self.rxnGroup.setLayout(rxnLayout)
self.tablelayout.addWidget(self.rxnGroup)
self.reactionTable.resizeColumnsToContents()
self.reactionTable.cellDoubleClicked.connect(self.HandleUpdateReaction)
def CreateCalibrationTable(self):
self.calGroup = QGroupBox("Calibration Peaks", self.tableTab)
calLayout = QVBoxLayout()
self.calibrationTable = QTableWidget(self.calGroup)
self.calibrationTable.setColumnCount(8)
self.calibrationTable.setHorizontalHeaderLabels(["Reaction","x(mm)","ux stat.(mm)","ux sys.(mm)","rho(cm)","urho(cm)","Ex(MeV)","uEx(MeV)"])
calLayout.addWidget(self.calibrationTable)
self.calGroup.setLayout(calLayout)
self.tablelayout.addWidget(self.calGroup)
self.calibrationTable.resizeColumnsToContents()
self.calibrationTable.cellDoubleClicked.connect(self.HandleUpdateCalibration)
def CreateOutputTable(self):
self.outGroup = QGroupBox("Output Peaks", self.tableTab)
outLayout = QVBoxLayout()
self.outputTable = QTableWidget(self.outGroup)
self.outputTable.setColumnCount(12)
self.outputTable.setHorizontalHeaderLabels(["Reaction","x(mm)","ux stat.(mm)","ux sys.(mm)","rho(cm)","urho(cm)","Ex(MeV)","uEx(MeV)","FWHM(mm)","uFWHM(mm)","FWHM(MeV)","uFWHM(MeV)"])
outLayout.addWidget(self.outputTable)
self.outGroup.setLayout(outLayout)
self.tablelayout.addWidget(self.outGroup)
self.outputTable.resizeColumnsToContents()
self.outputTable.cellDoubleClicked.connect(self.HandleUpdateOutput)
def CreateFitTable(self):
self.ftableGroup = QGroupBox("Fit Results", self.plotTab)
ftableLayout = QVBoxLayout()
self.fitTable = QTableWidget(3, 9, self.ftableGroup)
self.fitTable.setHorizontalHeaderLabels(["a0","ua0","a1","ua1","a2","ua2","a3","ua3","Chi Sq./NDF"])
self.fitTable.setVerticalHeaderLabels(["linear","quadratic","cubic"])
ftableLayout.addWidget(self.fitTable)
self.ftableGroup.setLayout(ftableLayout)
self.plotlayout.addWidget(self.ftableGroup,0,1)
self.fitTable.resizeColumnsToContents()
def CreateResidualTable(self):
self.rtableGroup = QGroupBox("Residual Results", self.plotTab)
rtableLayout = QVBoxLayout()
self.residualTable = QTableWidget(self.rtableGroup)
self.residualTable.setColumnCount(5)
self.residualTable.setHorizontalHeaderLabels(["x(mm)","rho calc(cm)","rho fit(cm)","residual(cm)","studentized residual"])
rtableLayout.addWidget(self.residualTable)
self.rtableGroup.setLayout(rtableLayout)
self.plotlayout.addWidget(self.rtableGroup,1,1)
self.residualTable.resizeColumnsToContents()
def HandleSave(self):
fileName = QFileDialog.getSaveFileName(self, "Save Input","./","Text Files (*.pickle)")
if fileName[0]:
#self.spanc.WriteConfig(fileName[0])
with open(fileName[0], "wb") as savefile:
pickle.dump(self.spanc, savefile, pickle.HIGHEST_PROTOCOL)
savefile.close()
def HandleSaveFit(self):
fileName = QFileDialog.getSaveFileName(self, "Save Fit Image","./","Image Files (*.png, *.eps)")
if fileName[0]:
self.fitCanvas.fig.savefig(fileName[0])
def HandleOpen(self):
fileName = QFileDialog.getOpenFileName(self, "Open Input","./","Text Files (*.pickle)")
if fileName[0]:
with open(fileName[0], "rb") as openfile:
self.spanc = pickle.load(openfile)
self.UpdateTargetTable()
self.UpdateReactionTable()
self.UpdateCalibrationTable()
self.UpdateOutputTable()
openfile.close()
def HandleNewTarget(self):
targDia = TargetDialog(self)
targDia.new_target.connect(self.AddTarget)
targDia.exec()
return
def HandleUpdateTarget(self, row, col):
targName = self.targetTable.verticalHeaderItem(row).text()
targDia = TargetDialog(self, target=self.spanc.targets[targName])
targDia.new_target.connect(self.UpdateTarget)
targDia.exec()
return
def HandleNewReaction(self):
rxnDia = ReactionDialog(self)
rxnDia.new_reaction.connect(self.AddReaction)
rxnDia.exec()
return
def HandleUpdateReaction(self, row, col):
rxnName = self.reactionTable.verticalHeaderItem(row).text()
rxnDia = ReactionDialog(self, rxn=self.spanc.reactions[rxnName], rxnKey=rxnName)
rxnDia.update_reaction.connect(self.UpdateReaction)
rxnDia.exec()
return
def HandleNewCalibration(self):
calDia = PeakDialog("Calibration", self)
calDia.new_calibration.connect(self.AddCalibration)
calDia.exec()
return
def HandleUpdateCalibration(self, row, col):
peakName = self.calibrationTable.verticalHeaderItem(row).text()
calDia = PeakDialog("Calibration", self, peak=self.spanc.calib_peaks[peakName], peakKey=peakName)
calDia.update_calibration.connect(self.UpdateCalibration)
calDia.exec()
return
def HandleNewOutput(self):
outDia = PeakDialog("Output", self)
outDia.new_output.connect(self.AddOutput)
outDia.exec()
return
def HandleUpdateOutput(self, row, col):
peakName = self.outputTable.verticalHeaderItem(row).text()
outDia = PeakDialog("Output",self,peak=self.spanc.output_peaks[peakName],peakKey=peakName)
outDia.update_output.connect(self.UpdateOutput)
outDia.exec()
return
def HandleRunFit(self):
fit_type = self.fitTypeBox.currentText()
npoints = len(self.spanc.calib_peaks)
if npoints < 3 and fit_type == "linear":
print("Warning! Too few points to properly fit a linear function. Results are invalid.")
elif npoints < 4 and fit_type == "quadratic":
print("Warning! Too few points to properly fit a quadratic function. Results are invalid")
elif npoints < 5 and fit_type == "cubic":
print("Warning! Too few points to properly fit a cubic function. Results are invalid")
xarray, yarray, uxarray, uyarray = self.spanc.PerformFits()
fitarray = np.linspace(np.amin(xarray), np.amax(xarray), 1000)
self.fitCanvas.axes.cla()
self.fitCanvas.axes.errorbar(xarray, yarray, yerr=uyarray, xerr=uxarray, marker="o", linestyle="None", elinewidth=2.0)
self.fitCanvas.axes.plot(fitarray, self.spanc.fitters[fit_type].EvaluateFunction(fitarray))
self.fitCanvas.axes.set_xlabel(r"$x$ (mm)")
self.fitCanvas.axes.set_ylabel(r"$\rho$ (cm)")
self.fitCanvas.draw()
self.UpdateFitTable()
self.spanc.CalculateOutputs(fit_type)
self.UpdateOutputTable()
self.fitFlag = True
def HandleRunResids(self):
fit_type = self.fitTypeBox.currentText()
npoints = len(self.spanc.calib_peaks)
if npoints < 3 and fit_type == "linear":
print("Warning! Too few points to properly fit a linear function. Results are invalid.")
elif npoints < 4 and fit_type == "quadratic":
print("Warning! Too few points to properly fit a quadratic function. Results are invalid")
elif npoints < 5 and fit_type == "cubic":
print("Warning! Too few points to properly fit a cubic function. Results are invalid")
xarray, resid_array, student_resids = self.spanc.CalculateResiduals(fit_type)
self.residCanvas.axes.cla()
self.residCanvas.axes.plot(xarray, resid_array, marker="o", linestyle="None")
self.residCanvas.axes.set_xlabel(r"$x$ (mm)")
self.residCanvas.axes.set_ylabel(r"Residual (cm)")
self.residCanvas.draw()
self.UpdateResidualTable(resid_array, student_resids)
def AddTarget(self, layers, name):
target = LayeredTarget()
target.name = name
for layer in layers:
target.AddLayer(layer)
self.spanc.targets[target.name] = target
self.UpdateTargetTable()
def UpdateTarget(self, layers, name):
target = LayeredTarget()
target.name = name
for layer in layers:
target.AddLayer(layer)
self.spanc.targets[target.name] = target
for rxn in self.spanc.reactions.values():
if rxn.target_data.name == name:
rxn.target_data = target
self.UpdateTargetTable()
self.spanc.CalculateCalibrations()
self.UpdateReactionTable()
self.UpdateCalibrationTable()
if self.fitFlag is True:
self.spanc.CalculateOutputs(self.fitTypeBox.currentText())
self.UpdateOutputTable()
def UpdateTargetTable(self):
self.targetTable.setRowCount(len(self.spanc.targets))
self.targetTable.setVerticalHeaderLabels(self.spanc.targets.keys())
cur_row = 0
for key in self.spanc.targets:
for i in range(len(self.spanc.targets[key].targets)) :
self.targetTable.setItem(cur_row, 0+i*2, QTableWidgetItem(str(self.spanc.targets[key].targets[i].thickness)))
self.targetTable.setItem(cur_row, 1+i*2, QTableWidgetItem(self.spanc.targets[key].targets[i].GetComposition()))
cur_row += 1
self.targetTable.resizeColumnsToContents()
self.targetTable.resizeRowsToContents()
def AddReaction(self, zt, at, zp, ap, ze, ae, bke, theta, bfield, name):
targ = self.spanc.targets[name]
rxn = Reaction(zt, at, zp, ap, ze, ae, bke, theta, bfield, targ)
count=0
for key in self.spanc.reactions:
if key == rxn.Name:
count += 1
rxn.Name = rxn.Name + "_" + str(count)
self.spanc.reactions[rxn.Name] = rxn
self.UpdateReactionTable()
def UpdateReaction(self, bke, theta, bfield, key):
self.spanc.reactions[key].ChangeReactionParameters(bke, theta, bfield)
self.UpdateReactionTable()
self.spanc.CalculateCalibrations()
self.UpdateCalibrationTable()
if self.fitFlag is True:
self.spanc.CalculateOutputs(self.fitTypeBox.currentText())
self.UpdateOutputTable()
def UpdateReactionTable(self):
self.reactionTable.setRowCount(len(self.spanc.reactions))
self.reactionTable.setVerticalHeaderLabels(self.spanc.reactions.keys())
cur_row = 0
for key in self.spanc.reactions:
self.reactionTable.setItem(cur_row, 0, QTableWidgetItem(str(self.spanc.reactions[key].target_data.name)))
self.reactionTable.setItem(cur_row, 1, QTableWidgetItem(str(self.spanc.reactions[key].Target.Z)))
self.reactionTable.setItem(cur_row, 2, QTableWidgetItem(str(self.spanc.reactions[key].Target.A)))
self.reactionTable.setItem(cur_row, 3, QTableWidgetItem(str(self.spanc.reactions[key].Projectile.Z)))
self.reactionTable.setItem(cur_row, 4, QTableWidgetItem(str(self.spanc.reactions[key].Projectile.A)))
self.reactionTable.setItem(cur_row, 5, QTableWidgetItem(str(self.spanc.reactions[key].Ejectile.Z)))
self.reactionTable.setItem(cur_row, 6, QTableWidgetItem(str(self.spanc.reactions[key].Ejectile.A)))
self.reactionTable.setItem(cur_row, 7, QTableWidgetItem(str(self.spanc.reactions[key].Residual.Z)))
self.reactionTable.setItem(cur_row, 8, QTableWidgetItem(str(self.spanc.reactions[key].Residual.A)))
self.reactionTable.setItem(cur_row, 9, QTableWidgetItem(str(self.spanc.reactions[key].BKE)))
self.reactionTable.setItem(cur_row, 10, QTableWidgetItem(str(self.spanc.reactions[key].Bfield)))
self.reactionTable.setItem(cur_row, 11, QTableWidgetItem(str(self.spanc.reactions[key].Theta/self.spanc.reactions[key].DEG2RAD)))
cur_row += 1
self.reactionTable.resizeColumnsToContents()
self.reactionTable.resizeRowsToContents()
def AddCalibration(self, x, uxstat, uxsys, ex, uex, rxnname):
peak_name = "Cal" + str(len(self.spanc.calib_peaks))
self.spanc.AddCalibrationPeak(rxnname, peak_name, x, uxstat, uxsys, ex, uex)
self.UpdateCalibrationTable()
def UpdateCalibration(self, x, uxstat, uxsys, ex, uex, rxnname, peakname):
self.spanc.AddCalibrationPeak(rxnname, peakname, x, uxstat, uxsys, ex, uex)
self.UpdateCalibrationTable()
def UpdateCalibrationTable(self):
self.calibrationTable.setRowCount(len(self.spanc.calib_peaks))
self.calibrationTable.setVerticalHeaderLabels(self.spanc.calib_peaks.keys())
cur_row = 0
for key in self.spanc.calib_peaks:
self.calibrationTable.setItem(cur_row, 0, QTableWidgetItem(self.spanc.calib_peaks[key].reaction))
self.calibrationTable.setItem(cur_row, 1, QTableWidgetItem(str(self.spanc.calib_peaks[key].x)))
self.calibrationTable.setItem(cur_row, 2, QTableWidgetItem(str(self.spanc.calib_peaks[key].ux_stat)))
self.calibrationTable.setItem(cur_row, 3, QTableWidgetItem(str(self.spanc.calib_peaks[key].ux_sys)))
self.calibrationTable.setItem(cur_row, 4, QTableWidgetItem(str(self.spanc.calib_peaks[key].rho)))
self.calibrationTable.setItem(cur_row, 5, QTableWidgetItem(str(self.spanc.calib_peaks[key].urho)))
self.calibrationTable.setItem(cur_row, 6, QTableWidgetItem(str(self.spanc.calib_peaks[key].Ex)))
self.calibrationTable.setItem(cur_row, 7, QTableWidgetItem(str(self.spanc.calib_peaks[key].uEx)))
cur_row += 1
self.calibrationTable.resizeColumnsToContents()
self.calibrationTable.resizeRowsToContents()
def AddOutput(self, x, uxstat, uxsys, fwhm, ufwhm, rxnname):
peak_name = "Out" + str(len(self.spanc.output_peaks))
self.spanc.AddOutputPeak(rxnname, peak_name, x, uxstat, uxsys, fwhm, ufwhm)
self.UpdateOutputTable()
def UpdateOutput(self, x, uxstat, uxsys, fwhm, ufwhm, rxnname, peakname):
self.spanc.AddOutputPeak(rxnname, peakname, x, uxstat, uxsys, fwhm, ufwhm)
self.UpdateOutputTable()
def UpdateOutputTable(self):
self.outputTable.setRowCount(len(self.spanc.output_peaks))
self.outputTable.setVerticalHeaderLabels(self.spanc.output_peaks.keys())
cur_row = 0
for key in self.spanc.output_peaks:
self.outputTable.setItem(cur_row, 0, QTableWidgetItem(self.spanc.output_peaks[key].reaction))
self.outputTable.setItem(cur_row, 1, QTableWidgetItem(str(self.spanc.output_peaks[key].x)))
self.outputTable.setItem(cur_row, 2, QTableWidgetItem(str(self.spanc.output_peaks[key].ux_stat)))
self.outputTable.setItem(cur_row, 3, QTableWidgetItem(str(self.spanc.output_peaks[key].ux_sys)))
self.outputTable.setItem(cur_row, 4, QTableWidgetItem(str(self.spanc.output_peaks[key].rho)))
self.outputTable.setItem(cur_row, 5, QTableWidgetItem(str(self.spanc.output_peaks[key].urho)))
self.outputTable.setItem(cur_row, 6, QTableWidgetItem(str(self.spanc.output_peaks[key].Ex)))
self.outputTable.setItem(cur_row, 7, QTableWidgetItem(str(self.spanc.output_peaks[key].uEx)))
self.outputTable.setItem(cur_row, 8, QTableWidgetItem(str(self.spanc.output_peaks[key].fwhm_x)))
self.outputTable.setItem(cur_row, 9, QTableWidgetItem(str(self.spanc.output_peaks[key].ufwhm_x)))
self.outputTable.setItem(cur_row, 10, QTableWidgetItem(str(self.spanc.output_peaks[key].fwhm_Ex)))
self.outputTable.setItem(cur_row, 11, QTableWidgetItem(str(self.spanc.output_peaks[key].ufwhm_Ex)))
cur_row += 1
self.outputTable.resizeColumnsToContents()
self.outputTable.resizeRowsToContents()
def UpdateFitTable(self):
cur_row=0
for key in self.spanc.fitters:
for i in range(len(self.spanc.fitters[key].parameters)):
self.fitTable.setItem(cur_row, 0+i*2, QTableWidgetItem(str(self.spanc.fitters[key].parameters[i])))
self.fitTable.setItem(cur_row, 1+i*2, QTableWidgetItem(str(self.spanc.fitters[key].GetParameterError(i))))
#self.fitTable.setItem(cur_row, 8, QTableWidgetItem(str(self.spanc.fitters[key].redChiSq)))
self.fitTable.setItem(cur_row, 8, QTableWidgetItem(str(self.spanc.fitters[key].ReducedChiSquare())))
cur_row += 1
self.fitTable.resizeColumnsToContents()
self.fitTable.resizeRowsToContents()
def UpdateResidualTable(self, resids, student_resids):
self.residualTable.setRowCount(len(self.spanc.calib_peaks))
self.residualTable.setVerticalHeaderLabels(self.spanc.calib_peaks.keys())
cur_row=0
for key in self.spanc.calib_peaks:
self.residualTable.setItem(cur_row, 0, QTableWidgetItem(str(self.spanc.calib_peaks[key].x)))
self.residualTable.setItem(cur_row, 1, QTableWidgetItem(str(self.spanc.calib_peaks[key].rho)))
self.residualTable.setItem(cur_row, 2, QTableWidgetItem(str(self.spanc.calib_peaks[key].rho + resids[cur_row])))
self.residualTable.setItem(cur_row, 3, QTableWidgetItem(str(resids[cur_row])))
self.residualTable.setItem(cur_row, 4, QTableWidgetItem(str(student_resids[cur_row])))
cur_row += 1
self.residualTable.resizeColumnsToContents()
self.residualTable.resizeRowsToContents()
def main() :
mpl.use("Qt5Agg")
myapp = QApplication(sys.argv)
window = SpancGUI()
sys.exit(myapp.exec_())
if __name__ == '__main__':
main()

70
spsplot/NucData.py
View File

@ -1,70 +0,0 @@
#!/usr/bin/env python3
import numpy as np
import requests
import lxml.html as xhtml
class MassTable:
def __init__(self):
file = open("./etc/mass.txt","r")
self.mtable = {}
u2mev = 931.4940954
me = 0.000548579909 #MeV
self.etable = {}
file.readline()
file.readline()
for line in file:
entries = line.split()
n = entries[0]
z = entries[1]
a = entries[2]
element = entries[3]
massBig = float(entries[4])
massSmall = float(entries[5])
key = '('+z+','+a+')'
value = (massBig+massSmall*1e-6)*u2mev - float(z)*me
self.mtable[key] = value
self.etable[key] = element
file.close()
def GetMass(self, z, a):
key = '('+str(z)+','+str(a)+')'
if key in self.mtable:
return self.mtable[key]
else:
return 0
def GetSymbol(self, z, a):
key = '('+str(z)+','+str(a)+')'
if key in self.etable:
return str(a)+self.etable[key]
else:
return 'none'
Masses = MassTable()
def GetExcitations(symbol):
levels = np.array(np.empty(0))
text = ''
site = requests.get("https://www.nndc.bnl.gov/nudat2/getdatasetClassic.jsp?nucleus="+symbol+"&unc=nds")
contents = xhtml.fromstring(site.content)
tables = contents.xpath("//table")
rows = tables[2].xpath("./tr")
for row in rows[1:-2]:
entries = row.xpath("./td")
if len(entries) != 0:
entry = entries[0]
data = entry.xpath("./a")
if len(data) == 0:
text = entry.text
else:
text = data[0].text
text = text.replace('?', '')
text = text.replace('\xa0\xa0≈','')
levels = np.append(levels, float(text)/1000.0)
return levels

83
spsplot/NuclearRxn.py
View File

@ -1,83 +0,0 @@
#!/usr/bin/env python3
import numpy as np
import NucData
class Nucleus:
def __init__(self, z, a):
self.Z = z
self.A = a
self.Symbol = NucData.Masses.GetSymbol(self.Z, self.A)
self.GSMass = NucData.Masses.GetMass(self.Z, self.A)
def Minus(self, rhs):
final_Z = self.Z - rhs.Z
final_A = self.A - rhs.A
if final_A < 0 or final_Z < 0:
print("Illegal minus operation on Nuclei!")
return Nucleus(0,0)
else:
return Nucleus(final_Z, final_A)
def Plus(self, rhs):
return Nucleus(self.Z + rhs.Z, self.A + rhs.A)
class Reaction:
DEG2RAD = np.pi/180.0 #degrees to radians
C = 299792458 #speed of light m/s
QBRHO2P = 1.0E-9*C #Converts qbrho to p (kG*cm -> MeV/c)
def __init__(self, zt, at, zp, ap, ze, ae, beamKE, theta, bfield):
self.Target = Nucleus(zt, at)
self.Projectile = Nucleus(zp, ap)
self.Ejectile = Nucleus(ze, ae)
self.Residual = (self.Target.Plus(self.Projectile)).Minus(self.Ejectile)
self.BKE = beamKE
self.Theta = theta * self.DEG2RAD
self.Bfield = bfield
self.Name = self.Target.Symbol +"("+ self.Projectile.Symbol +","+ self.Ejectile.Symbol +")"+ self.Residual.Symbol
self.residLevels = NucData.GetExcitations(self.Residual.Symbol)
self.ejectKEvals = np.array(np.empty(len(self.residLevels)))
self.ejectRhovals = np.array(np.empty(len(self.residLevels)))
self.SetEjectileData()
def GetEjectileKineticEnergy(self, Elevel) :
Q = self.Target.GSMass + self.Projectile.GSMass - (self.Ejectile.GSMass + self.Residual.GSMass + Elevel)
Ethresh = -Q*(self.Ejectile.GSMass+self.Residual.GSMass)/(self.Ejectile.GSMass + self.Residual.GSMass - self.Projectile.GSMass)
if self.BKE < Ethresh:
return 0.0
term1 = np.sqrt(self.Projectile.GSMass*self.Ejectile.GSMass*self.BKE)/(self.Ejectile.GSMass + self.Residual.GSMass)*np.cos(self.Theta)
term2 = (self.BKE*(self.Residual.GSMass - self.Projectile.GSMass) + self.Residual.GSMass*Q)/(self.Ejectile.GSMass + self.Residual.GSMass)
ke1 = term1 + np.sqrt(term1**2.0 + term2)
ke2 = term1 - np.sqrt(term1**2.0 + term2)
if ke1 > 0:
return ke1**2.0
else :
return ke2**2.0
def GetEjectileRho(self, ke):
p = np.sqrt(ke*(ke + 2.0*self.Ejectile.GSMass))
return p/(self.QBRHO2P*self.Bfield*self.Ejectile.Z)
def SetEjectileData(self):
for index in range(len(self.residLevels)):
self.ejectKEvals[index] = self.GetEjectileKineticEnergy(self.residLevels[index])
self.ejectRhovals[index] = self.GetEjectileRho(self.ejectKEvals[index])
def ChangeReactionParameters(self, bke, theta, bf) :
self.BKE = bke
self.Theta = theta*self.DEG2RAD
self.Bfield = bf
self.SetEjectileData()
def AddLevel(self, Elevel):
ke = self.GetEjectileKineticEnergy(Elevel)
rho = self.GetEjectileRho(ke)
self.residLevels = np.append(self.residLevels, Elevel)
self.ejectKEvals = np.append(self.ejectKEvals, ke)
self.ejectRhovals = np.append(self.ejectRhovals, rho)

76
spsplot/SPSPlot.py
View File

@ -1,76 +0,0 @@
#!/usr/bin/env python3
import NuclearRxn as rxn
class SPSPlot:
def __init__(self) :
self.reactions = {}
self.configfile = ""
self.rhoMin = 0
self.rhoMax = 99
self.beamKE = 0
self.Bfield = 0
self.angle = 0
def ReadConfig(self, filename) :
self.reactions.clear()
self.configfile = filename
file = open(filename, "r")
line = file.readline()
entries = line.split()
self.beamKE = float(entries[1])
line = file.readline()
entries = line.split()
self.Bfield = float(entries[1])
line = file.readline()
entries = line.split()
self.angle = float(entries[1])
line = file.readline()
entries = line.split()
self.rhoMin = float(entries[1])
self.rhoMax = float(entries[3])
line = file.readline()
for line in file:
entries = line.split()
reac = rxn.Reaction(int(entries[0]), int(entries[1]), int(entries[2]), int(entries[3]), int(entries[4]), int(entries[5]), self.beamKE, self.angle, self.Bfield)
self.reactions[reac.Name] = reac
file.close()
def WriteConfig(self, filename) :
file = open(filename, "w")
line = "BeamEnergy(MeV): "+str(self.beamKE)+"\n"
file.write(line)
line = "B-field(kG): "+str(self.Bfield)+"\n"
file.write(line)
line = "Angle(deg): "+str(self.angle)+"\n"
file.write(line)
line = "RhoMin: "+str(self.rhoMin)+" RhoMax: "+str(self.rhoMax)+"\n"
file.write(line)
line = "ZT AT ZP AP ZE AE\n"
file.write(line)
for rxnName in self.reactions :
reaction = self.reactions[rxnName]
line = str(reaction.Target.Z)+" "+str(reaction.Target.A)+" "+str(reaction.Projectile.Z)+" "+str(reaction.Projectile.A)+" "+str(reaction.Ejectile.Z)+" "+str(reaction.Ejectile.A)+"\n"
file.write(line)
file.close()
def ChangeReactionParameters(self, bke, theta, bf) :
self.beamKE = bke
self.Bfield = bf
self.angle = theta
for rxnName in self.reactions :
self.reactions[rxnName].ChangeReactionParameters(bke, theta, bf)
def AddReaction(self, zt, at, zp, ap, ze, ae) :
reac = rxn.Reaction(zt, at, zp, ap, ze, ae, self.beamKE, self.Bfield, self.angle)
self.reactions[reac.Name] = reac
def AddLevel(self, name, level) :
self.reactions[name].AddLevel(level)

295
spsplot/SPSPlotGUI.py
View File

@ -1,295 +0,0 @@
#!/usr/bin/env python3
import SPSPlot as spsplt
import sys
from qtpy.QtWidgets import QApplication, QWidget, QMainWindow
from qtpy.QtWidgets import QLabel, QMenuBar, QAction
from qtpy.QtWidgets import QHBoxLayout, QVBoxLayout, QGroupBox
from qtpy.QtWidgets import QPushButton, QButtonGroup, QRadioButton
from qtpy.QtWidgets import QSpinBox, QDoubleSpinBox, QComboBox
from qtpy.QtWidgets import QDialog, QFileDialog, QDialogButtonBox
from qtpy.QtCore import Signal
import matplotlib as mpl
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg
from matplotlib.figure import Figure
class MPLCanvas(FigureCanvasQTAgg):
def __init__(self, parent=None, width=5, height=4, dpi=100):
self.fig = Figure(figsize=(width, height), dpi=dpi, edgecolor="black",linewidth=0.5)
self.axes = self.fig.add_subplot(111)
self.axes.spines['top'].set_visible(False)
super(MPLCanvas, self).__init__(self.fig)
class ReactionDialog(QDialog):
new_reaction = Signal(int, int, int, int, int, int)
def __init__(self, parent=None):
super().__init__(parent)
self.setWindowTitle("Add A Reaction")
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
self.buttonBox.accepted.connect(self.SendReaction)
self.buttonBox.rejected.connect(self.reject)
self.layout = QVBoxLayout()
self.setLayout(self.layout)
self.CreateReactionInputs()
self.layout.addWidget(self.buttonBox)
def SendReaction(self) :
self.new_reaction.emit(self.ztInput.value(),self.atInput.value(),self.zpInput.value(),self.apInput.value(),self.zeInput.value(),self.aeInput.value())
def CreateReactionInputs(self) :
self.nucleiGroupBox = QGroupBox("Reaction Nuclei",self)
inputLayout = QVBoxLayout()
ztLabel = QLabel("ZT",self.nucleiGroupBox)
self.ztInput = QSpinBox(self.nucleiGroupBox)
self.ztInput.setRange(1, 110)
atLabel = QLabel("AT",self.nucleiGroupBox)
self.atInput = QSpinBox(self.nucleiGroupBox)
self.atInput.setRange(1,270)
zpLabel = QLabel("ZP",self.nucleiGroupBox)
self.zpInput = QSpinBox(self.nucleiGroupBox)
self.zpInput.setRange(1, 110)
apLabel = QLabel("AP",self.nucleiGroupBox)
self.apInput = QSpinBox(self.nucleiGroupBox)
self.apInput.setRange(1,270)
zeLabel = QLabel("ZE",self.nucleiGroupBox)
self.zeInput = QSpinBox(self.nucleiGroupBox)
self.zeInput.setRange(1, 110)
aeLabel = QLabel("AE",self.nucleiGroupBox)
self.aeInput = QSpinBox(self.nucleiGroupBox)
self.aeInput.setRange(1,270)
inputLayout.addWidget(ztLabel)
inputLayout.addWidget(self.ztInput)
inputLayout.addWidget(atLabel)
inputLayout.addWidget(self.atInput)
inputLayout.addWidget(zpLabel)
inputLayout.addWidget(self.zpInput)
inputLayout.addWidget(apLabel)
inputLayout.addWidget(self.apInput)
inputLayout.addWidget(zeLabel)
inputLayout.addWidget(self.zeInput)
inputLayout.addWidget(aeLabel)
inputLayout.addWidget(self.aeInput)
self.nucleiGroupBox.setLayout(inputLayout)
self.layout.addWidget(self.nucleiGroupBox)
class LevelDialog(QDialog):
new_level = Signal(str,float)
def __init__(self, parent) :
super().__init__(parent)
self.setWindowTitle("Add a Level")
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
self.buttonBox.accepted.connect(self.SendLevel)
self.buttonBox.rejected.connect(self.reject)
self.layout = QVBoxLayout()
self.setLayout(self.layout)
rxnLabel = QLabel("Choose a reaction",self)
self.reactionList = QComboBox(self)
for rxnName in parent.sps.reactions:
self.reactionList.addItem(rxnName)
stateLabel = QLabel("New state energy",self)
self.stateInput = QDoubleSpinBox(self)
self.stateInput.setRange(0.0,40.0)
self.stateInput.setSuffix(" MeV")
self.layout.addWidget(rxnLabel)
self.layout.addWidget(self.reactionList)
self.layout.addWidget(stateLabel)
self.layout.addWidget(self.stateInput)
self.layout.addWidget(self.buttonBox)
def SendLevel(self):
self.new_level.emit(self.reactionList.currentText(),self.stateInput.value())
class SPSPlotGUI(QMainWindow):
def __init__(self, parent=None) :
super().__init__(parent)
self.setWindowTitle("SPSPlot")
self.sps = spsplt.SPSPlot()
self.generalLayout = QVBoxLayout()
self.centralWidget = QWidget(self)
self.setCentralWidget(self.centralWidget)
self.centralWidget.setLayout(self.generalLayout)
self.energyFlag = True #True = ex False = ke
self.CreateCanvas()
self.CreateMenus()
self.CreateInputs()
self.show()
def CreateCanvas(self):
self.canvas = MPLCanvas(self, width=14, height=5, dpi=100)
self.generalLayout.addWidget(self.canvas, 5)
def CreateMenus(self):
self.fileMenu = self.menuBar().addMenu("&File")
saveAction = QAction("&Save...",self)
openAction = QAction("&Open...",self)
self.fileMenu.addAction(saveAction)
self.fileMenu.addAction(openAction)
self.fileMenu.addAction("&Exit", self.close)
saveAction.triggered.connect(self.HandleSave)
openAction.triggered.connect(self.HandleOpen)
self.addMenu = self.menuBar().addMenu("&New")
newStateAction = QAction("New state...", self)
newReactionAction = QAction("New reaction...", self)
self.addMenu.addAction(newStateAction)
self.addMenu.addAction(newReactionAction)
newStateAction.triggered.connect(self.HandleNewState)
newReactionAction.triggered.connect(self.HandleNewReaction)
def CreateInputs(self):
inputLayout = QHBoxLayout()
self.inputGroupBox = QGroupBox("Adjustable Inputs", self)
rhoMinLabel = QLabel("Rho Min", self.inputGroupBox)
self.rhoMinInput = QDoubleSpinBox(self.inputGroupBox)
self.rhoMinInput.setRange(0.0, 150.0)
self.rhoMinInput.setSuffix(" cm")
rhoMaxLabel = QLabel("RhoMax", self.inputGroupBox)
self.rhoMaxInput = QDoubleSpinBox(self.inputGroupBox)
self.rhoMaxInput.setRange(0.0,150.0)
self.rhoMaxInput.setSuffix(" cm")
bkeLabel = QLabel("Beam KE", self.inputGroupBox)
self.bkeInput = QDoubleSpinBox(self.inputGroupBox)
self.bkeInput.setRange(0.0, 500.0)
self.bkeInput.setSuffix(" MeV")
bfieldLabel = QLabel("B-field", self.inputGroupBox)
self.bfieldInput = QDoubleSpinBox(self.inputGroupBox)
self.bfieldInput.setRange(0.0, 17.0)
self.bfieldInput.setSuffix(" kG")
angleLabel = QLabel("Angle", self.inputGroupBox)
self.angleInput = QDoubleSpinBox(self.inputGroupBox)
self.angleInput.setRange(0.0, 180.0)
self.angleInput.setSuffix(" deg")
self.runButton = QPushButton("Run", self.inputGroupBox)
self.runButton.clicked.connect(self.HandleRun)
self.energyButtonGroup = QGroupBox("Ex/KE switch",self)
buttonLayout = QHBoxLayout()
self.exButton = QRadioButton("Excitation energy", self.energyButtonGroup)
self.exButton.toggled.connect(self.HandleExSwitch)
self.keButton = QRadioButton("Ejectile Kinetic energy", self.energyButtonGroup)
self.keButton.toggled.connect(self.HandleKESwitch)
buttonLayout.addWidget(self.exButton)
buttonLayout.addWidget(self.keButton)
self.energyButtonGroup.setLayout(buttonLayout)
inputLayout.addWidget(rhoMinLabel)
inputLayout.addWidget(self.rhoMinInput)
inputLayout.addWidget(rhoMaxLabel)
inputLayout.addWidget(self.rhoMaxInput)
inputLayout.addWidget(bkeLabel)
inputLayout.addWidget(self.bkeInput)
inputLayout.addWidget(bfieldLabel)
inputLayout.addWidget(self.bfieldInput)
inputLayout.addWidget(angleLabel)
inputLayout.addWidget(self.angleInput)
inputLayout.addWidget(self.runButton)
self.inputGroupBox.setLayout(inputLayout)
inputLayout.addWidget(self.energyButtonGroup)
self.generalLayout.addWidget(self.inputGroupBox, 1)
def HandleSave(self):
fileName = QFileDialog.getSaveFileName(self, "Save Input","./","Text Files (*.txt *.inp)")
if fileName[0]:
self.sps.WriteConfig(fileName[0])
def HandleOpen(self):
fileName = QFileDialog.getOpenFileName(self, "Open Input","./","Text Files (*.txt *.inp)")
if fileName[0]:
self.sps.ReadConfig(fileName[0])
self.UpdateInputs()
self.UpdatePlot()
def HandleNewState(self):
stDlg = LevelDialog(self)
stDlg.new_level.connect(self.sps.AddLevel)
if stDlg.exec():
self.UpdatePlot()
def HandleNewReaction(self):
rxnDlg = ReactionDialog(self)
rxnDlg.new_reaction.connect(self.sps.AddReaction)
if rxnDlg.exec():
self.UpdatePlot()
def HandleRun(self):
self.sps.ChangeReactionParameters(self.bkeInput.value(), self.angleInput.value(), self.bfieldInput.value())
self.sps.rhoMin = self.rhoMinInput.value()
self.sps.rhoMax = self.rhoMaxInput.value()
self.UpdatePlot()
def HandleExSwitch(self):
if self.exButton.isChecked() and (not self.energyFlag):
self.energyFlag = True
self.UpdatePlot()
def HandleKESwitch(self):
if self.keButton.isChecked() and self.energyFlag:
self.energyFlag = False
self.UpdatePlot()
def UpdatePlot(self):
rxnNumber = 0
rhos = []
exs = []
kes = []
rxns = []
for rxnName in self.sps.reactions:
rxnNumber += 1
rxn = self.sps.reactions[rxnName]
for i in range(len(rxn.residLevels)):
rxns.append(rxnNumber)
rhos.append(rxn.ejectRhovals[i])
exs.append(rxn.residLevels[i])
kes.append(rxn.ejectKEvals[i])
self.canvas.axes.cla()
self.canvas.axes.plot(rhos, rxns, marker="o", linestyle="None")
for i in range(len(rxns)):
y = rxns[i]
x = rhos[i]
label = ''
if self.energyFlag:
label = "{:.2f}".format(exs[i])
else:
label = "{:.2f}".format(kes[i])
self.canvas.axes.annotate(label, (x,y), textcoords="offset points",xytext=(0,10),ha="center",rotation="90")
self.canvas.axes.set_xlim(self.sps.rhoMin, self.sps.rhoMax)
self.canvas.axes.set_yticks(range(1,rxnNumber+1))
self.canvas.axes.set_yticklabels(self.sps.reactions)
self.canvas.draw()
def UpdateInputs(self):
self.rhoMinInput.setValue(self.sps.rhoMin)
self.rhoMaxInput.setValue(self.sps.rhoMax)
self.bfieldInput.setValue(self.sps.Bfield)
self.bkeInput.setValue(self.sps.beamKE)
self.angleInput.setValue(self.sps.angle)
def main() :
mpl.use("Qt5Agg")
myapp = QApplication(sys.argv)
window = SPSPlotGUI()
sys.exit(myapp.exec_())
if __name__ == '__main__':
main()

138
spspy/Fitter.py Normal file
View File

@ -0,0 +1,138 @@
import numpy as np
from numpy.typing import NDArray
from numpy.polynomial import Polynomial
from scipy import odr
from dataclasses import dataclass
from typing import Optional
INVALID_FIT_RESULT = np.inf
INVALID_NDF = -1
@dataclass
class FitPoint:
x: float = 0.0
y: float = 0.0
xError: float = 0.0
yError: float = 0.0
def convert_fit_points_to_arrays(data: list[FitPoint]) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
xArray = np.empty(len(data))
yArray = np.empty(len(data))
xErrorArray = np.empty(len(data))
yErrorArray = np.empty(len(data))
for index, point in enumerate(data):
xArray[index] = point.x
yArray[index] = point.y
xErrorArray[index] = point.xError
yErrorArray[index] = point.yError
return xArray, yArray, xErrorArray, yErrorArray
@dataclass
class FitResidual:
x: float = 0.0
residual: float = 0.0
studentizedResidual: float = 0.0
def convert_resid_points_to_arrays(data: list[FitResidual]) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
xArray = np.empty(len(data))
residArray = np.empty(len(data))
studentResidArray = np.empty(len(data))
for index, point in enumerate(data):
xArray[index] = point.x
residArray[index] = point.residual
studentResidArray[index] = point.studentizedResidual
return xArray, residArray, studentResidArray
class Fitter:
def __init__(self, order: int =1):
self.polynomialOrder: int = order
self.fitResults: Optional[odr.Output] = None
self.fitData: Optional[list[FitPoint]] = None
self.function: Optional[Polynomial] = None
def set_polynomial_order(self, order: int) -> None:
self.polynomialOrder = order
def run(self, data: list[FitPoint] = None) -> None:
if data is not None:
self.fitData = data
if self.fitData is not None:
xArray, yArray, xErrorArray, yErrorArray = convert_fit_points_to_arrays(self.fitData)
modelData = odr.RealData(xArray, y=yArray, sx=xErrorArray, sy=yErrorArray)
model = odr.polynomial(self.polynomialOrder)
self.fitResults = odr.ODR(modelData, model).run()
self.function = Polynomial(self.fitResults.beta)
else:
print("Cannot run fitter without setting data to be fit!")
def get_parameters(self) -> NDArray[np.float64] :
if self.fitResults is not None:
return self.fitResults.beta
return np.array({INVALID_FIT_RESULT})
def get_parameter_errors(self) -> NDArray[np.float64] :
if self.fitResults is not None:
return self.fitResults.sd_beta
return np.array({INVALID_FIT_RESULT})
def get_ndf(self) -> int:
if self.fitResults is not None:
return len(self.fitData) - 1
return INVALID_NDF
def evaluate(self, x: float) -> float:
if self.function is not None:
return self.function(x)
return INVALID_FIT_RESULT
def evaluate_derivative(self, x: float) -> float:
if self.function is not None:
return self.function.deriv()(x)
return INVALID_FIT_RESULT
def evaluate_param_derivative(self, x: float, index: int) -> float:
if self.fitResults is not None and len(self.fitResults.beta) > index:
return x**index
return INVALID_FIT_RESULT
def get_chisquare(self) -> float:
if self.function is not None:
yEffErrorArray = [np.sqrt(point.yError**2.0 + (point.xError * self.evaluate_derivative(point.x))**2.0) for point in self.fitData]
chisq = 0.0
for index, point in enumerate(self.fitData):
chisq = ((point.y - self.evaluate(point.x)) / yEffErrorArray[index])**2.0
return chisq
return INVALID_FIT_RESULT
def get_reduced_chisquare(self) -> float:
ndf = self.get_ndf()
chisq = self.get_chisquare()
if chisq == INVALID_FIT_RESULT or ndf == INVALID_NDF:
return INVALID_FIT_RESULT
else:
return chisq/ndf
def get_residuals(self) -> list[FitResidual]:
if self.fitData is not None:
fitResiduals = [FitResidual(point.x, point.y - self.evaluate(point.x), 0.0) for point in self.fitData]
#compute the leverage and studentize
xMean = 0.0
rmse = 0.0
npoints = len(fitResiduals)
for resid in fitResiduals:
xMean += resid.x
rmse += resid.residual**2.0
xMean /= npoints
rmse /= self.get_ndf()
meanDiffSq = 0.0
for resid in fitResiduals:
meanDiffSq += (resid.x - xMean) ** 2.0
meanDiffSq /= npoints
for resid in fitResiduals:
leverage = 1.0/npoints + (resid.x - xMean)/meanDiffSq
resid.studentizedResidual = resid.residual / (rmse * np.sqrt(1.0 - leverage))
return fitResiduals
return []

48
spspy/Launcher.py Normal file
View File

@ -0,0 +1,48 @@
from PySide6.QtWidgets import QApplication, QWidget, QMainWindow
from PySide6.QtWidgets import QLabel
from PySide6.QtWidgets import QHBoxLayout, QVBoxLayout, QGroupBox
from PySide6.QtWidgets import QPushButton
from PySide6.QtGui import QAction
from .SPSPlotUI import run_spsplot_ui, SPSPlotGUI
from .SpancUI import run_spanc_ui, SpancGUI
import sys
import matplotlib as mpl
from qdarktheme import load_stylesheet
class Launcher(QMainWindow):
def __init__(self, parent=None):
super().__init__(parent)
self.setWindowTitle("SPSPY Launcher")
self.mainLayout = QVBoxLayout()
self.mainWidget = QWidget(self)
self.setCentralWidget(self.mainWidget)
self.mainWidget.setLayout(self.mainLayout)
self.spsplotButton = QPushButton("Launch SPSPlot", self.mainWidget)
self.spsplotButton.clicked.connect(self.handle_spsplot)
self.spancButton = QPushButton("Launch SPANC", self.mainWidget)
self.spancButton.clicked.connect(self.handle_spanc)
self.mainLayout.addWidget(self.spsplotButton)
self.mainLayout.addWidget(self.spancButton)
self.show()
def handle_spsplot(self) -> None:
SPSPlotGUI(self)
def handle_spanc(self) -> None:
#run_spanc_ui()
SpancGUI(self)
def run_launcher() -> None:
mpl.use("Qt5Agg")
app = QApplication.instance()
if not app:
app = QApplication(sys.argv)
app.setStyleSheet(load_stylesheet())
window = Launcher()
sys.exit(app.exec_())

71
spspy/SPSPlot.py Normal file
View File

@ -0,0 +1,71 @@
from .SPSReaction import *
from .SPSTarget import *
from .data.NuclearData import *
from dataclasses import dataclass, field
import csv
@dataclass
class Excitation:
excitation: float = 0.0 #MeV
kineticEnergy: float = 0.0 #MeV
rho: float = 0.0 #cm
fpZ: float = 0.0 #cm
@dataclass
class PlotData:
rxn: Reaction
excitations: list[Excitation] = field(default_factory=list)
class SPSPlot:
def __init__(self):
self.data : dict[str, PlotData] = {}
self.targets : dict[str, SPSTarget] = {}
self.spsAngle : float = 0.0 #deg
self.beamEnergy : float = 0.0 #MeV
self.magneticField : float = 0.0 #kG
self.rhoMin: float = 0.0 #cm
self.rhoMax: float = 0.0 #cm
def add_target(self, targName: str, layers: list[TargetLayer]) -> None:
self.targets[targName] = SPSTarget(layers, name=targName)
def add_reaction(self, params: RxnParameters, targetName: str) -> None:
target = self.targets.get(targetName, None)
if target is None:
target = SPSTarget([TargetLayer([(1,1,1)], 0.0)]) #insert a dummy target if an invalid one is passed (i.e. no energy loss)
plotData = PlotData(Reaction(params, target))
exList = get_excitations(plotData.rxn.residual.Z, plotData.rxn.residual.A)
for ex in exList:
ke = plotData.rxn.calculate_ejectile_KE(ex)
r = plotData.rxn.convert_ejectile_KE_2_rho(ke)
z = plotData.rxn.calculate_focal_plane_offset(ke)
plotData.excitations.append(Excitation(ex, ke, r, z))
self.data[str(plotData.rxn)] = plotData
def update_reactions(self) -> None:
for datum in self.data.values():
datum.rxn.update_parameters(self.beamEnergy, self.spsAngle, self.magneticField)
if datum.rxn.targetMaterial.name in self.targets:
datum.rxn.targetMaterial = self.targets[datum.rxn.targetMaterial.name]
for ex in datum.excitations:
ex.kineticEnergy = datum.rxn.calculate_ejectile_KE(ex.excitation)
ex.rho = datum.rxn.convert_ejectile_KE_2_rho(ex.kineticEnergy)
def add_excitation(self, rxnName: str, excitation: float) -> None:
if rxnName not in self.data:
print("Cannot add excitation to non-existant reaction named ", rxnName)
return
datum = self.data[rxnName]
ke = datum.rxn.calculate_ejectile_KE(excitation)
rho = datum.rxn.convert_ejectile_KE_2_rho(ke)
datum.excitations.append(Excitation(excitation, ke, rho))
def export_reaction_data(self, fileName: str) -> None:
with open(fileName, "w", newline='') as outputFile:
csvWriter = csv.writer(outputFile, dialect="excel")
csvWriter.writerow(["Reaction", "Excitation(MeV)", "Rho(cm)", "EjectileKE(MeV)", "Z-Offset(cm)"])
for datum in self.data.values():
for point in datum.excitations:
csvWriter.writerow([repr(datum.rxn), f"{point.excitation:.3f}", f"{point.rho:.3f}", f"{point.kineticEnergy:.3f}", f"{point.fpZ:.3f}"])

304
spspy/SPSPlotUI.py Normal file
View File

@ -0,0 +1,304 @@
from .SPSPlot import SPSPlot
from .SPSReaction import RxnParameters
from .ui.MPLCanvas import MPLCanvas
from .ui.ReactionDialog import ReactionDialog
from .ui.TargetDialog import TargetDialog
from .ui.ExcitationDialog import ExcitationDialog
from PySide6.QtWidgets import QApplication, QWidget, QMainWindow
from PySide6.QtWidgets import QLabel, QTabWidget, QTableWidget, QTableWidgetItem
from PySide6.QtWidgets import QHBoxLayout, QVBoxLayout, QGroupBox
from PySide6.QtWidgets import QPushButton, QRadioButton
from PySide6.QtWidgets import QDoubleSpinBox
from PySide6.QtWidgets import QFileDialog
from PySide6.QtGui import QAction
from qdarktheme import load_stylesheet
from enum import Enum, auto
import matplotlib as mpl
import sys
import pickle
DEFAULT_RHO_MIN: float = 69.0
DEFAULT_RHO_MAX: float = 87.0
class PlotType(Enum):
PLOT_EX = auto()
PLOT_KE = auto()
PLOT_Z = auto()
class SPSPlotGUI(QMainWindow):
def __init__(self, parent=None) :
super().__init__(parent)
self.setWindowTitle("SPSPlot")
self.sps = SPSPlot()
self.plotLayout = QVBoxLayout()
self.targetLayout = QVBoxLayout()
self.generalLayout = QVBoxLayout()
self.centralWidget = QTabWidget(self)
self.setCentralWidget(self.centralWidget)
self.centralWidget.setLayout(self.generalLayout)
self.targetTab = QWidget(self.centralWidget)
self.targetTab.setLayout(self.targetLayout)
self.plotTab = QWidget(self.centralWidget)
self.plotTab.setLayout(self.plotLayout)
self.centralWidget.addTab(self.plotTab, "Plot")
self.centralWidget.addTab(self.targetTab, "Targets")
self.plotType = PlotType.PLOT_EX
self.create_menus()
self.create_canvas()
self.create_inputs()
self.create_target_table()
self.update_plot()
self.show()
def create_canvas(self) -> None:
self.canvas = MPLCanvas(self.plotTab, width=14, height=5, dpi=100)
self.plotLayout.addWidget(self.canvas, 4)
def create_menus(self) -> None:
self.fileMenu = self.menuBar().addMenu("&File")
saveAction = QAction("&Save...",self)
openAction = QAction("&Open...",self)
self.fileMenu.addAction(saveAction)
self.fileMenu.addAction(openAction)
self.fileMenu.addAction("&Exit", self.close)
saveAction.triggered.connect(self.handle_save)
openAction.triggered.connect(self.handle_open)
self.addMenu = self.menuBar().addMenu("&New")
newTargetAction = QAction("New target...", self)
newReactionAction = QAction("New reaction...", self)
newStateAction = QAction("New state...", self)
self.addMenu.addAction(newTargetAction)
self.addMenu.addAction(newReactionAction)
self.addMenu.addAction(newStateAction)
newStateAction.triggered.connect(self.handle_new_state)
newReactionAction.triggered.connect(self.handle_new_reaction)
newTargetAction.triggered.connect(self.handle_new_target)
self.exportMenu = self.menuBar().addMenu("&Export")
exportLevels = QAction("Export levels to csv...", self)
self.exportMenu.addAction(exportLevels)
exportLevels.triggered.connect(self.handle_export_levels)
def create_inputs(self) -> None:
inputLayout = QVBoxLayout()
self.inputGroupBox = QGroupBox("Adjustable Inputs", self.plotTab)
self.spsGroupBox = QGroupBox("SPS Parameters", self.inputGroupBox)
spsGroupLayout = QHBoxLayout()
rhoMinLabel = QLabel("<p>&rho;<sub>Min<\sub><\p>", self.spsGroupBox)
self.rhoMinInput = QDoubleSpinBox(self.spsGroupBox)
self.rhoMinInput.setRange(0.0, 150.0)
self.rhoMinInput.setValue(DEFAULT_RHO_MIN)
self.sps.rhoMin = DEFAULT_RHO_MIN
self.rhoMinInput.setSuffix(" cm")
rhoMaxLabel = QLabel("<p>&rho;<sub>Max<\sub><\p>", self.spsGroupBox)
self.rhoMaxInput = QDoubleSpinBox(self.spsGroupBox)
self.rhoMaxInput.setRange(0.0,150.0)
self.rhoMaxInput.setValue(DEFAULT_RHO_MAX)
self.sps.rhoMax = DEFAULT_RHO_MAX
self.rhoMaxInput.setSuffix(" cm")
bkeLabel = QLabel("<p>E<sub>beam<\sub><\p>", self.spsGroupBox)
self.bkeInput = QDoubleSpinBox(self.spsGroupBox)
self.bkeInput.setRange(0.0, 500.0)
self.bkeInput.setSuffix(" MeV")
bfieldLabel = QLabel("B", self.spsGroupBox)
self.bfieldInput = QDoubleSpinBox(self.spsGroupBox)
self.bfieldInput.setRange(0.0, 17.0)
self.bfieldInput.setSuffix(" kG")
angleLabel = QLabel("<p>&theta;<sub>SPS<\sub><\p>", self.spsGroupBox)
self.angleInput = QDoubleSpinBox(self.spsGroupBox)
self.angleInput.setRange(0.0, 180.0)
self.angleInput.setSuffix(" deg")
self.runButton = QPushButton("Set", self.spsGroupBox)
self.runButton.clicked.connect(self.handle_run)
spsGroupLayout.addWidget(rhoMinLabel, 1)
spsGroupLayout.addWidget(self.rhoMinInput, 2)
spsGroupLayout.addWidget(rhoMaxLabel,1 )
spsGroupLayout.addWidget(self.rhoMaxInput, 2)
spsGroupLayout.addWidget(bkeLabel, 1)
spsGroupLayout.addWidget(self.bkeInput, 2)
spsGroupLayout.addWidget(bfieldLabel, 1)
spsGroupLayout.addWidget(self.bfieldInput, 2)
spsGroupLayout.addWidget(angleLabel, 1)
spsGroupLayout.addWidget(self.angleInput, 2)
spsGroupLayout.addWidget(self.runButton, 1)
self.spsGroupBox.setLayout(spsGroupLayout)
self.energyButtonGroup = QGroupBox("Labels",self.plotTab)
buttonLayout = QHBoxLayout()
self.exButton = QRadioButton("Excitation Energy(MeV)", self.energyButtonGroup)
self.exButton.toggled.connect(self.handle_ex_switch)
self.exButton.toggle()
self.keButton = QRadioButton("Ejectile KE(MeV)", self.energyButtonGroup)
self.keButton.toggled.connect(self.handle_ke_switch)
self.zButton = QRadioButton("FocalPlane Z-Shift(cm)", self.energyButtonGroup)
self.zButton.toggled.connect(self.handle_z_switch)
buttonLayout.addWidget(self.exButton)
buttonLayout.addWidget(self.keButton)
buttonLayout.addWidget(self.zButton)
self.energyButtonGroup.setLayout(buttonLayout)
inputLayout.addWidget(self.spsGroupBox)
inputLayout.addWidget(self.energyButtonGroup)
self.inputGroupBox.setLayout(inputLayout)
self.plotLayout.addWidget(self.inputGroupBox, 1)
def create_target_table(self) -> None:
self.targetGroup = QGroupBox("Targets", self.targetTab)
tableLayout = QVBoxLayout()
self.targetTable = QTableWidget(self.targetGroup)
self.targetTable.setColumnCount(6)
self.targetTable.setHorizontalHeaderLabels(["L1 Thickness(ug/cm^2)", "L1 Compound",
"L2 Thickness(ug/cm^2)", "L2 Compound",
"L3 Thickness(ug/cm^2)", "Layer3 Compound"])
tableLayout.addWidget(self.targetTable)
self.targetGroup.setLayout(tableLayout)
self.targetLayout.addWidget(self.targetGroup)
self.targetTable.resizeColumnsToContents()
self.targetTable.cellDoubleClicked.connect(self.handle_update_target)
def handle_save(self) -> None:
fileName = QFileDialog.getSaveFileName(self, "Save Input","./","SPSPlot Files (*.sps)")
if fileName[0]:
with open(fileName[0], "wb") as file:
pickle.dump(self.sps, file, pickle.HIGHEST_PROTOCOL)
def handle_open(self) -> None:
fileName = QFileDialog.getOpenFileName(self, "Open Input","./","SPSPlot Files (*.sps)")
if fileName[0]:
with open(fileName[0], "rb") as file:
self.sps = pickle.load(file)
self.update_inputs()
self.update_plot()
self.update_target_table()
def handle_new_state(self) -> None:
stDlg = ExcitationDialog(self, self.sps.data.keys())
stDlg.new_level.connect(self.sps.add_excitation)
if stDlg.exec():
self.update_plot()
def handle_new_reaction(self) -> None:
rxnDlg = ReactionDialog(parent=self, targets=self.sps.targets.keys())
rxnDlg.new_reaction.connect(self.add_reaction)
rxnDlg.exec()
def handle_new_target(self) -> None:
targDlg = TargetDialog(self)
targDlg.new_target.connect(self.sps.add_target)
if targDlg.exec():
self.update_target_table()
def handle_update_target(self, row: int, col: int) -> None:
targName = self.targetTable.verticalHeaderItem(row).text()
targDia = TargetDialog(self, target=self.sps.targets[targName])
targDia.new_target.connect(self.sps.add_target)
if targDia.exec():
self.update_target_table()
self.sps.update_reactions()
self.update_plot() #in case a reaction is using the target
def handle_run(self) -> None:
self.sps.beamEnergy = self.bkeInput.value()
self.sps.spsAngle = self.angleInput.value()
self.sps.magneticField = self.bfieldInput.value()
self.sps.update_reactions()
self.sps.rhoMin = self.rhoMinInput.value()
self.sps.rhoMax = self.rhoMaxInput.value()
self.update_plot()
def handle_ex_switch(self) -> None:
if self.exButton.isChecked() and self.plotType != PlotType.PLOT_EX:
self.plotType = PlotType.PLOT_EX
self.update_plot()
def handle_ke_switch(self) -> None:
if self.keButton.isChecked() and self.plotType != PlotType.PLOT_KE:
self.plotType = PlotType.PLOT_KE
self.update_plot()
def handle_z_switch(self) -> None:
if self.zButton.isChecked() and self.plotType != PlotType.PLOT_Z:
self.plotType = PlotType.PLOT_Z
self.update_plot()
def handle_export_levels(self) -> None:
fileName = QFileDialog.getSaveFileName(self, "Export Levels to CSV","./","Comma-Separated Values File (*.csv)")
if fileName[0]:
self.sps.export_reaction_data(fileName[0])
def add_reaction(self, rxnParams: RxnParameters, targName: str) -> None:
rxnParams.beamEnergy = self.bkeInput.value()
rxnParams.spsAngle = self.angleInput.value()
rxnParams.magneticField = self.bfieldInput.value()
self.sps.add_reaction(rxnParams, targName)
self.update_plot()
def update_plot(self) -> None:
rhos = []
exs = []
kes = []
zs = []
rxns = []
for rxnNumber, rxn in enumerate(self.sps.data.values()):
for point in rxn.excitations:
rxns.append(rxnNumber+1)
rhos.append(point.rho)
exs.append(point.excitation)
kes.append(point.kineticEnergy)
zs.append(point.fpZ)
self.canvas.axes.cla()
self.canvas.axes.plot(rhos, rxns, marker="o", linestyle="None")
for i in range(len(rxns)):
y = rxns[i]
x = rhos[i]
label = ''
if self.plotType == PlotType.PLOT_EX:
label = "{:.2f}".format(exs[i])
elif self.plotType == PlotType.PLOT_KE:
label = "{:.2f}".format(kes[i])
elif self.plotType == PlotType.PLOT_Z:
label = "{:.2f}".format(zs[i])
self.canvas.axes.annotate(label, (x,y), textcoords="offset points",xytext=(0,10),ha="center",rotation="vertical")
ylabels = [r.rxn.get_latex_rep() for r in self.sps.data.values()]
ylabels.append("Reactions")
self.canvas.axes.set_xlim(self.sps.rhoMin, self.sps.rhoMax)
self.canvas.axes.set_yticks(range(1,len(self.sps.data)+2))
self.canvas.axes.set_yticklabels(ylabels)
self.canvas.axes.set_xlabel(r"$\rho$ (cm)")
self.canvas.draw()
def update_inputs(self):
self.rhoMinInput.setValue(self.sps.rhoMin)
self.rhoMaxInput.setValue(self.sps.rhoMax)
self.bfieldInput.setValue(self.sps.magneticField)
self.bkeInput.setValue(self.sps.beamEnergy)
self.angleInput.setValue(self.sps.spsAngle)
def update_target_table(self):
self.targetTable.setRowCount(len(self.sps.targets))
self.targetTable.setVerticalHeaderLabels(self.sps.targets.keys())
for row, key in enumerate(self.sps.targets):
for col, layer in enumerate(self.sps.targets[key].layer_details) :
self.targetTable.setItem(row, col*2, QTableWidgetItem(str(layer.thickness)))
self.targetTable.setCellWidget(row, 1+col*2, QLabel(str(layer)))
self.targetTable.resizeColumnsToContents()
self.targetTable.resizeRowsToContents()
def run_spsplot_ui():
mpl.use("Qt5Agg")
app = QApplication.instance()
if not app:
app = QApplication(sys.argv)
app.setStyleSheet(load_stylesheet())
window = SPSPlotGUI()
sys.exit(app.exec_())

109
spspy/SPSReaction.py Normal file
View File

@ -0,0 +1,109 @@
from .data.NuclearData import global_nuclear_data, NucleusData
from .SPSTarget import SPSTarget
from dataclasses import dataclass
from numpy import sqrt, cos, pi, sin
INVALID_KINETIC_ENERGY: float = -1000.0
@dataclass
class RxnParameters:
target: NucleusData
projectile: NucleusData
ejectile: NucleusData
beamEnergy: float = 0.0 #MeV
magneticField: float = 0.0 #kG
spsAngle: float = 0.0 #rad
def create_reaction_parameters(zt: int, at: int, zp: int, ap: int, ze: int, ae: int) -> RxnParameters:
return RxnParameters(global_nuclear_data.get_data(zt, at), global_nuclear_data.get_data(zp, ap), global_nuclear_data.get_data(ze, ae))
class Reaction:
DEG2RAD: float = pi/180.0 #degrees -> radians
C = 299792458 #speed of light m/s
QBRHO2P = 1.0E-9*C #Converts qbrho to momentum (p) (kG*cm -> MeV/c)
FP_MAGNIFICATION = 0.39
FP_DISPERSION = 1.96
def __init__(self, params: RxnParameters, target: SPSTarget):
self.params = params
self.targetMaterial = target
self.rxnLayer = self.targetMaterial.get_rxn_layer(self.params.target.Z, self.params.target.A)
self.setup_nuclei()
def setup_nuclei(self) -> None:
residZ = self.params.target.Z + self.params.projectile.Z - self.params.ejectile.Z
residA = self.params.target.A + self.params.projectile.A - self.params.ejectile.A
assert residZ > 0 and residA > 0, "Unable to construct residual in Reaction!"
self.residual = global_nuclear_data.get_data(residZ, residA)
self.Qvalue = self.params.target.mass + self.params.projectile.mass - self.params.ejectile.mass - self.residual.mass
def __str__(self) -> str:
return f"{self.params.target.prettyIsotopicSymbol}({self.params.projectile.prettyIsotopicSymbol},{self.params.ejectile.prettyIsotopicSymbol}){self.residual.prettyIsotopicSymbol}"
def __repr__(self) -> str:
return f"{self.params.target.isotopicSymbol}({self.params.projectile.isotopicSymbol},{self.params.ejectile.isotopicSymbol}){self.residual.isotopicSymbol}_{self.params.beamEnergy}MeV_{self.params.spsAngle}deg_{self.params.magneticField}kG"
def get_latex_rep(self) -> str:
return f"{self.params.target.get_latex_rep()}({self.params.projectile.get_latex_rep()},{self.params.ejectile.get_latex_rep()}){self.residual.get_latex_rep()}"
#MeV
def calculate_ejectile_KE(self, excitation: float) -> float:
rxnQ = self.Qvalue - excitation
angleRads = self.params.spsAngle * self.DEG2RAD
beamRxnEnergy = self.params.beamEnergy - self.targetMaterial.get_incoming_energyloss(self.params.projectile.Z, self.params.projectile.mass, self.params.beamEnergy, self.rxnLayer, 0.0)
threshold = -rxnQ*(self.params.ejectile.mass+self.residual.mass)/(self.params.ejectile.mass + self.residual.mass - self.params.projectile.mass)
if beamRxnEnergy < threshold:
return INVALID_KINETIC_ENERGY
term1 = sqrt(self.params.projectile.mass * self.params.ejectile.mass * beamRxnEnergy) / (self.params.ejectile.mass + self.residual.mass) * cos(angleRads * self.DEG2RAD)
term2 = (beamRxnEnergy * (self.residual.mass - self.params.projectile.mass) + self.residual.mass * rxnQ) / (self.params.ejectile.mass + self.residual.mass)
ke1 = term1 + sqrt(term1**2.0 + term2)
ke2 = term1 + sqrt(term1**2.0 + term2)
ejectileEnergy = 0.0
if ke1 > 0.0:
ejectileEnergy = ke1**2.0
else:
ejectileEnergy = ke2**2.0
ejectileEnergy -= self.targetMaterial.get_outgoing_energyloss(self.params.ejectile.Z, self.params.ejectile.mass, ejectileEnergy, self.rxnLayer, angleRads)
return ejectileEnergy
def convert_ejectile_KE_2_rho(self, ejectileEnergy: float) -> float:
if ejectileEnergy == INVALID_KINETIC_ENERGY:
return 0.0
p = sqrt( ejectileEnergy * (ejectileEnergy + 2.0 * self.params.ejectile.mass))
#convert to QBrho
qbrho = p/self.QBRHO2P
return qbrho / (float(self.params.ejectile.Z) * self.params.magneticField)
def calculate_excitation(self, rho: float) -> float:
angleRads = self.params.spsAngle * self.DEG2RAD
ejectileP = rho * float(self.params.ejectile.Z) * self.params.magneticField * self.QBRHO2P
ejectileEnergy = sqrt(ejectileP**2.0 + self.params.ejectile.mass**2.0) - self.params.ejectile.mass
ejectileRxnEnergy = ejectileEnergy + self.targetMaterial.get_outgoing_reverse_energyloss(self.params.ejectile.Z, self.params.ejectile.mass, ejectileEnergy, self.rxnLayer, angleRads)
ejectileRxnP = sqrt(ejectileRxnEnergy * (ejectileRxnEnergy + 2.0 * self.params.ejectile.mass))
beamRxnEnergy = self.params.beamEnergy - self.targetMaterial.get_incoming_energyloss(self.params.projectile.Z, self.params.projectile.mass, self.params.beamEnergy, self.rxnLayer, 0.0)
beamRxnP = sqrt(beamRxnEnergy * (beamRxnEnergy + 2.0 * self.params.projectile.mass))
residRxnEnergy = beamRxnEnergy + self.params.projectile.mass + self.params.target.mass - ejectileRxnEnergy - self.params.ejectile.mass
residRxnP2 = beamRxnP**2.0 + ejectileRxnP**2.0 - 2.0 * ejectileRxnP * beamRxnP * cos(angleRads)
return sqrt(residRxnEnergy**2.0 - residRxnP2) - self.residual.mass
def calculate_focal_plane_offset(self, ejectileEnergy: float) -> float:
if ejectileEnergy == INVALID_KINETIC_ENERGY:
return 0.0
ejectileRho = self.convert_ejectile_KE_2_rho(ejectileEnergy)
k = sqrt(self.params.projectile.mass * self.params.ejectile.mass * self.params.beamEnergy / ejectileEnergy) * sin(self.params.spsAngle * self.DEG2RAD)
k /= self.params.ejectile.mass + self.residual.mass - sqrt(self.params.projectile.mass * self.params.ejectile.mass * self.params.beamEnergy/ejectileEnergy) * cos(self.params.spsAngle * self.DEG2RAD)
return -1.0*k*ejectileRho*self.FP_DISPERSION*self.FP_MAGNIFICATION
def update_parameters(self, beamEnergy: float, spsAngle: float, magenticField: float):
self.params.beamEnergy = beamEnergy
self.params.spsAngle = spsAngle
self.params.magneticField = magenticField

171
spspy/SPSTarget.py Normal file
View File

@ -0,0 +1,171 @@
from .data.NuclearData import global_nuclear_data
import pycatima as catima
from dataclasses import dataclass, field
from numpy import pi, cos
INVALID_RXN_LAYER: int = -1
ADAPTIVE_DEPTH_MAX: int = 100
ENERGY_PERCENT_STEP_MIN: float = 0.001
@dataclass
class TargetLayer:
compound_list: list[tuple[int, int , int]] = field(default_factory=list) #Z,A,S
thickness: float = 0.0 #ug/cm^2
def __str__(self) -> str:
return "".join([f"{global_nuclear_data.get_data(z, a,).prettyIsotopicSymbol}<sub>{s}<\sub>" for z, a, s in self.compound_list])
#integrate energy loss starting from the final energy and running backwards to initial energy
#catima does not natively provide this type of method
#returns the total energy loss (really in this case energy gain) through the material
def get_reverse_energyloss(projectile: catima.Projectile, material: catima.Material) -> float:
depth = 0
e_out = projectile.T() #MeV/u
e_initial = e_out
x_step = 0.25*material.thickness() #g/cm^2
x_traversed = 0.0
e_step = catima.dedx(projectile, material)*x_step
A_recip = 1.0/projectile.A()
if material.thickness() <= 0.0:
return 0.0
#The integration step is adaptive, so use while(true)
while(True):
if e_step/e_initial > ENERGY_PERCENT_STEP_MIN and depth < ADAPTIVE_DEPTH_MAX:
depth += 1
x_step *= 0.5
e_step = catima.dedx(projectile, material)*x_step*A_recip
elif (x_step + x_traversed) >= material.thickness():
x_step = material.thickness() - x_traversed
e_step = catima.dedx(projectile, material)*x_step
e_initial += e_step*A_recip
projectile.T(e_initial)
return (e_initial - e_out)*projectile.A()
elif depth == ADAPTIVE_DEPTH_MAX:
return e_out*projectile.A()
else:
e_step = catima.dedx(projectile, material)*x_step
e_initial += e_step*A_recip
projectile.T(e_initial)
x_traversed += x_step
#integrate energy loss starting from the initial energy to final energy
#catima does not natively provide this type of method, only a calculate function which does a whole bunch of other stuff too
#returns the total energy loss through the material
def get_energyloss(projectile: catima.Projectile, material: catima.Material) -> float:
depth = 0
e_in = projectile.T() # MeV/u
e_final = e_in
x_step = 0.25*material.thickness() #g/cm^2
x_traversed = 0.0
e_step = catima.dedx(projectile, material)*x_step
A_recip = 1.0/projectile.A()
if material.thickness() <= 0.0:
return 0.0
while(True):
if e_step/e_final > ENERGY_PERCENT_STEP_MIN and depth < ADAPTIVE_DEPTH_MAX:
depth += 1
x_step *= 0.5
e_step = catima.dedx(projectile, material)*x_step*A_recip
elif (x_step + x_traversed) >= material.thickness():
x_step = material.thickness() - x_traversed
e_step = catima.dedx(projectile, material)*x_step*A_recip
e_final -= e_step
projectile.T(e_final)
return (e_in - e_final)*projectile.A()
elif depth == ADAPTIVE_DEPTH_MAX:
return e_in*projectile.A()
else:
e_step = catima.dedx(projectile, material)*x_step*A_recip
e_final -= e_step
projectile.T(e_final)
x_traversed += x_step
class SPSTarget:
UG2G: float = 1.0e-6 #convert ug to g
def __init__(self, layers: list[TargetLayer], name: str = "default"):
self.layer_details = layers
self.name = name
def __str__(self):
return self.name
def get_rxn_layer(self, zt: int, at: int) -> int:
for idx, layer in enumerate(self.layer_details):
for (a, z, s) in layer.compound_list:
if at == a and zt == z:
return idx
return INVALID_RXN_LAYER
#Calculate energy loss for a particle coming into the target, up to rxn layer (halfway through rxn layer)
def get_incoming_energyloss(self, zp: int, ap: float, e_initial: float, rxn_layer: int, angle: float) -> float:
if angle == pi*0.5:
return e_initial
projectile = catima.Projectile(zp, ap)
e_current = e_initial/ap
for (idx, layer) in enumerate(self.layer_details):
material = catima.Material([(global_nuclear_data.get_data(z, a).mass, z, float(s)) for (z, a, s) in layer.compound_list])
projectile.T(e_current) #catima wants MeV/u
if idx == rxn_layer:
material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle))))
e_current -= get_energyloss(projectile, material)
return e_initial - e_current*ap
else:
material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle)))
e_current -= get_energyloss(projectile, material)
return e_initial - e_current*ap
#Calculate energy loss for a particle leaving the target, from rxn layer (halfway through rxn layer) to end
def get_outgoing_energyloss(self, zp: int, ap: float, e_initial: float, rxn_layer: int, angle: float) -> float:
if angle == pi*0.5:
return e_initial
projectile = catima.Projectile(zp, ap)
e_current = e_initial/ap
for (idx, layer) in enumerate(self.layer_details[rxn_layer:], start=rxn_layer):
material = catima.Material([(global_nuclear_data.get_data(z, a).mass, z, float(s)) for (z, a, s) in layer.compound_list])
projectile.T(e_current) #catima wants MeV/u
if idx == rxn_layer:
material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle))))
else:
material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle)))
e_current -= get_energyloss(projectile, material)
return e_initial - e_current*ap
#Calculate reverse energy loss (energy gain) for a particle that left the target after a reaction (end -> rxn_layer)
def get_outgoing_reverse_energyloss(self, zp: int, ap: float, e_final: float, rxn_layer: int, angle: float) -> float:
if angle == pi*0.5:
return 0.0
projectile = catima.Projectile(zp, ap)
e_current = e_final/ap
sublist = self.layer_details[rxn_layer:] #only care about rxn_layer -> exit
reveresedRxnLayer = len(sublist) -1 #when reversed rxn_layer is the last layer
for (idx, layer) in reversed(list(enumerate(sublist))):
material = catima.Material([(global_nuclear_data.get_data(z, a).mass, z, float(s)) for (z, a, s) in layer.compound_list])
projectile.T(e_current) #catima wants MeV/u
if idx == reveresedRxnLayer:
material.thickness(self.layer_details[idx].thickness * self.UG2G / (2.0 * abs(cos(angle))))
else:
material.thickness(self.layer_details[idx].thickness * self.UG2G / abs(cos(angle)))
e_current += get_reverse_energyloss(projectile, material)
return e_current*ap - e_final

121
spspy/Spanc.py Normal file
View File

@ -0,0 +1,121 @@
from .SPSReaction import *
from .SPSTarget import *
from .Fitter import *
from .data.NuclearData import *
from dataclasses import dataclass, field
import numpy as np
from enum import Enum
INVALID_PEAK_ID: int = -1
class PeakType(Enum):
CALIBRATION = "Calibration"
OUTPUT = "Output"
@dataclass
class Peak:
excitation: float = 0.0 #MeV
excitationErr: float = 0.0 #MeV
position: float = 0.0 #arb
positionErrStat: float = 0.0 #arb
positionErrSys: float = 0.0 #arb
rho: float = 0.0 #cm
rhoErr: float = 0.0 #cm
positionFWHM: float = 0.0
positionFWHMErr: float = 0.0
excitationFWHM: float = 0.0
excitationFWHMErr: float = 0.0
rxnName: str = ""
peakID: int = INVALID_PEAK_ID
class Spanc:
def __init__(self):
self.targets: dict[str, SPSTarget] = {}
self.reactions: dict[str, Reaction] = {}
self.calibrations: dict[int, Peak] = {}
self.outputs: dict[int, Peak] = {}
self.fitter : Fitter = Fitter()
self.isFit: bool = False
def set_fit_order(self, order: int) -> None:
self.fitter.set_polynomial_order(order)
#return fit data so that the data points can be drawn
def fit(self) -> list[FitPoint]:
fitData = [FitPoint(peak.position, peak.rho, np.sqrt(peak.positionErrStat**2.0 + peak.positionErrSys**2.0), peak.rhoErr) for peak in self.calibrations.values()]
self.fitter.run(data=fitData)
self.isFit = True
return fitData
def get_residuals(self) -> list[FitResidual]:
return self.fitter.get_residuals()
def add_target(self, targName: str, layers: list[TargetLayer]) -> None:
self.targets[targName] = SPSTarget(layers, name=targName)
def add_reaction(self, params: RxnParameters, targetName: str) -> None:
if targetName not in self.targets:
print("Cannot create reaction with non-existant target ", targetName)
return
key = f"Rxn{len(self.reactions)}"
rxn = Reaction(params, target=self.targets[targetName])
self.reactions[key] = rxn
def update_reaction_parameters(self, beamEnergy: float, spsAngle: float, magneticField: float, rxnName: str):
if rxnName in self.reactions:
rxn = self.reactions[rxnName]
rxn.params.beamEnergy = beamEnergy
rxn.params.spsAngle = spsAngle
rxn.params.magneticField = magneticField
def add_calibration(self, data: Peak) -> None:
if data.rxnName in self.reactions:
rxn = self.reactions[data.rxnName]
data.rho = rxn.convert_ejectile_KE_2_rho(rxn.calculate_ejectile_KE(data.excitation))
data.rhoErr = np.abs(rxn.convert_ejectile_KE_2_rho(rxn.calculate_ejectile_KE(data.excitation + data.excitationErr)) - data.rho)
if data.peakID == INVALID_PEAK_ID:
data.peakID = len(self.calibrations)
self.calibrations[data.peakID] = data
return
def add_output(self, data: Peak) -> None:
if data.peakID == INVALID_PEAK_ID:
data.peakID = len(self.outputs)
self.outputs[data.peakID] = data
return
def calculate_output_urho(self, peak: Peak) -> float:
urho = 0.0
paramErrors = self.fitter.get_parameter_errors()
for i, paramErr in enumerate(paramErrors):
urho += (self.fitter.evaluate_param_derivative(peak.position, i) * paramErr)**2.0
urho += (self.fitter.evaluate_derivative(peak.position)*np.sqrt(peak.positionErrStat**2.0 + peak.positionErrSys**2.0))**2.0
return np.sqrt(urho)
def calculate_outputs(self) -> None:
if self.isFit == False:
return
for output in self.outputs.values():
rxn = self.reactions[output.rxnName]
output.rho = self.fitter.evaluate(output.position)
output.rhoErr = self.calculate_output_urho(output)
output.excitation = rxn.calculate_excitation(output.rho)
output.excitationErr = np.abs(rxn.calculate_excitation(output.rho + output.rhoErr) - output.excitation)
if output.positionFWHM == 0:
output.excitationFWHM = 0
output.excitationFWHMErr = 0
else:
rhoLo = self.fitter.evaluate(output.position - output.positionFWHM * 0.5)
rhoHi = self.fitter.evaluate(output.position + output.positionFWHM * 0.5)
exLo = rxn.calculate_excitation(rhoLo)
exHi = rxn.calculate_excitation(rhoHi)
output.excitationFWHM = abs(exHi - exLo)
output.excitationFWHMErr = output.positionFWHMErr/output.positionFWHM*output.excitationFWHM
def calculate_calibrations(self) -> None:
for calibration in self.calibrations.values():
rxn = self.reactions[calibration.rxnName]
calibration.rho = rxn.convert_ejectile_KE_2_rho(rxn.calculate_ejectile_KE(calibration.excitation))
calibration.rhoErr = np.abs(rxn.convert_ejectile_KE_2_rho(rxn.calculate_ejectile_KE(calibration.excitation + calibration.excitationErr)) - calibration.rho)

375
spspy/SpancUI.py Normal file
View File

@ -0,0 +1,375 @@
from .Spanc import Spanc, PeakType
from .Fitter import convert_fit_points_to_arrays, convert_resid_points_to_arrays
from .ui.MPLCanvas import MPLCanvas
from .ui.ReactionDialog import ReactionDialog
from .ui.TargetDialog import TargetDialog
from .ui.PeakDialog import PeakDialog
from PySide6.QtWidgets import QApplication, QWidget, QMainWindow
from PySide6.QtWidgets import QLabel, QTabWidget, QTableWidget, QTableWidgetItem
from PySide6.QtWidgets import QHBoxLayout, QVBoxLayout, QGroupBox
from PySide6.QtWidgets import QPushButton, QTextEdit, QSpinBox
from PySide6.QtWidgets import QFileDialog
from PySide6.QtGui import QAction
from qdarktheme import load_stylesheet, load_palette
import matplotlib as mpl
import numpy as np
from numpy.typing import NDArray
import sys
import pickle
#Get y-value for baseline at 0
def baseline(x: float) -> float:
return 0.0
class SpancGUI(QMainWindow):
def __init__(self, parent=None):
super().__init__(parent)
self.setWindowTitle("SPANC")
self.spanc = Spanc()
self.tablelayout = QVBoxLayout()
self.plotlayout = QHBoxLayout()
self.layout = QVBoxLayout()
self.centralWidget = QTabWidget(self)
self.setCentralWidget(self.centralWidget)
self.centralWidget.setLayout(self.layout)
self.tableTab = QWidget(self.centralWidget)
self.tableTab.setLayout(self.tablelayout)
self.plotTab = QWidget(self.centralWidget)
self.plotTab.setLayout(self.plotlayout)
self.centralWidget.addTab(self.tableTab, "Data Tables")
self.centralWidget.addTab(self.plotTab, "Plots and Fits")
self.fitFlag = False
self.create_menus()
self.create_fit_canvas()
self.create_target_table()
self.create_reaction_table()
self.create_calibration_table()
self.create_output_table()
self.create_fit_result_text()
self.show()
def create_menus(self) -> None:
self.fileMenu = self.menuBar().addMenu("&File")
saveAction = QAction("&Save...",self)
openAction = QAction("&Open...",self)
saveFitAction = QAction("Save Fit Plot...", self)
self.fileMenu.addAction(saveAction)
self.fileMenu.addAction(openAction)
self.fileMenu.addAction(saveFitAction)
self.fileMenu.addAction("&Exit", self.close)
saveAction.triggered.connect(self.handle_save)
openAction.triggered.connect(self.handle_open)
saveFitAction.triggered.connect(self.handle_save_fit)
self.addMenu = self.menuBar().addMenu("&New")
newTargetAction = QAction("New target...", self)
newReactionAction = QAction("New reaction...", self)
newCalibrationAction = QAction("New calibration...", self)
newOutputAction = QAction("New output...", self)
self.addMenu.addAction(newTargetAction)
self.addMenu.addAction(newReactionAction)
self.addMenu.addAction(newCalibrationAction)
self.addMenu.addAction(newOutputAction)
newTargetAction.triggered.connect(self.handle_new_target)
newReactionAction.triggered.connect(self.handle_new_reaction)
newCalibrationAction.triggered.connect(self.handle_new_calibration)
newOutputAction.triggered.connect(self.handle_new_output)
def create_fit_canvas(self) -> None:
self.fitGroup = QGroupBox("Calibration Fit", self.plotTab)
fitLayout = QVBoxLayout()
self.fitCanvas = MPLCanvas(self.fitGroup, width=6, height=6, dpi=100)
self.residCanvas = MPLCanvas(self.fitGroup, width=6, height=6, dpi=100)
self.fitOptionGroup = QGroupBox("Fit options", self.fitGroup)
fitOptionLayout = QHBoxLayout()
self.fitButton = QPushButton("Run Fit", self.fitOptionGroup)
self.fitButton.clicked.connect(self.handle_run_fit)
self.fitOrderBox = QSpinBox(self.fitOptionGroup)
self.fitOrderBox.valueChanged.connect(self.handle_change_fit_order)
self.fitOrderBox.setValue(1)
self.fitOrderBox.setMaximum(10)
self.fitOrderBox.setMinimum(0)
fitOptionLayout.addWidget(QLabel("Polynomial Order", self.fitOptionGroup))
fitOptionLayout.addWidget(self.fitOrderBox)
fitOptionLayout.addWidget(self.fitButton)
self.fitOptionGroup.setLayout(fitOptionLayout)
fitLayout.addWidget(QLabel("Fit", self.fitCanvas))
fitLayout.addWidget(self.fitCanvas)
fitLayout.addWidget(QLabel("Residuals", self.fitCanvas))
fitLayout.addWidget(self.residCanvas)
fitLayout.addWidget(self.fitOptionGroup)
self.fitGroup.setLayout(fitLayout)
self.plotlayout.addWidget(self.fitGroup)
def create_target_table(self) -> None:
self.targetGroup = QGroupBox("Targets", self.tableTab)
targetLayout = QVBoxLayout()
self.targetTable = QTableWidget(self.targetGroup)
self.targetTable.setColumnCount(6)
self.targetTable.setHorizontalHeaderLabels(["L1 Thickness(ug/cm^2)", "L1 Compound","L2 Thickness(ug/cm^2)", "L2 Compound","L3 Thickness(ug/cm^2)", "L3 Compound"])
targetLayout.addWidget(self.targetTable)
self.targetGroup.setLayout(targetLayout)
self.tablelayout.addWidget(self.targetGroup)
self.targetTable.resizeColumnsToContents()
self.targetTable.cellDoubleClicked.connect(self.handle_update_target)
def create_reaction_table(self) -> None:
self.rxnGroup = QGroupBox("Reactions", self.tableTab)
rxnLayout = QVBoxLayout()
self.reactionTable = QTableWidget(self.rxnGroup)
self.reactionTable.setColumnCount(5)
self.reactionTable.setHorizontalHeaderLabels(["Target Material","Reaction Equation","Beam KE(MeV)","BField(kG)","Angle(deg)"])
rxnLayout.addWidget(self.reactionTable)
self.rxnGroup.setLayout(rxnLayout)
self.tablelayout.addWidget(self.rxnGroup)
self.reactionTable.resizeColumnsToContents()
self.reactionTable.cellDoubleClicked.connect(self.handle_update_reaction)
def create_calibration_table(self) -> None:
self.calGroup = QGroupBox("Calibration Peaks", self.tableTab)
calLayout = QVBoxLayout()
self.calibrationTable = QTableWidget(self.calGroup)
self.calibrationTable.setColumnCount(8)
self.calibrationTable.setHorizontalHeaderLabels(["Reaction","x(mm)","ux stat.(mm)","ux sys.(mm)","rho(cm)","urho(cm)","Ex(MeV)","uEx(MeV)"])
calLayout.addWidget(self.calibrationTable)
self.calGroup.setLayout(calLayout)
self.tablelayout.addWidget(self.calGroup)
self.calibrationTable.resizeColumnsToContents()
self.calibrationTable.cellDoubleClicked.connect(self.handle_update_calibration)
def create_output_table(self) -> None:
self.outGroup = QGroupBox("Output Peaks", self.tableTab)
outLayout = QVBoxLayout()
self.outputTable = QTableWidget(self.outGroup)
self.outputTable.setColumnCount(12)
self.outputTable.setHorizontalHeaderLabels(["Reaction","x(mm)","ux stat.(mm)","ux sys.(mm)","rho(cm)","urho(cm)","Ex(MeV)","uEx(MeV)","FWHM(mm)","uFWHM(mm)","FWHM(MeV)","uFWHM(MeV)"])
outLayout.addWidget(self.outputTable)
self.outGroup.setLayout(outLayout)
self.tablelayout.addWidget(self.outGroup)
self.outputTable.resizeColumnsToContents()
self.outputTable.cellDoubleClicked.connect(self.handle_update_output)
def create_fit_result_text(self) -> None:
self.fitTextGroup = QGroupBox("Fit Results", self.plotTab)
fitTextLayout = QVBoxLayout()
self.fitResultText = QTextEdit(self.fitTextGroup)
self.fitResultText.setReadOnly(True)
fitTextLayout.addWidget(self.fitResultText)
self.fitTextGroup.setLayout(fitTextLayout)
self.plotlayout.addWidget(self.fitTextGroup)
def handle_save(self) -> None:
fileName = QFileDialog.getSaveFileName(self, "Save Input","./","SPANC Files (*.spanc)")
if fileName[0]:
#self.spanc.WriteConfig(fileName[0])
with open(fileName[0], "wb") as savefile:
pickle.dump(self.spanc, savefile, pickle.HIGHEST_PROTOCOL)
savefile.close()
def handle_save_fit(self) -> None:
fileName = QFileDialog.getSaveFileName(self, "Save Fit Image","./","Image Files (*.png, *.eps)")
if fileName[0]:
self.fitCanvas.fig.savefig(fileName[0])
def handle_open(self) -> None:
fileName = QFileDialog.getOpenFileName(self, "Open Input","./","SPANC Files (*.spanc)")
if fileName[0]:
with open(fileName[0], "rb") as openfile:
self.spanc = pickle.load(openfile)
self.update_target_table()
self.update_reaction_table()
self.update_calibration_table()
self.update_fit_order()
self.update_output_table()
openfile.close()
def handle_new_target(self) -> None:
targDia = TargetDialog(self)
targDia.new_target.connect(self.spanc.add_target)
if targDia.exec() :
self.update_target_table()
return
def handle_update_target(self, row: int, col: int) -> None:
targName = self.targetTable.verticalHeaderItem(row).text()
targDia = TargetDialog(self, target=self.spanc.targets[targName])
targDia.new_target.connect(self.spanc.add_target)
if targDia.exec():
self.update_target_table()
self.spanc.calculate_calibrations()
self.update_reaction_table()
self.update_calibration_table()
self.spanc.calculate_outputs()
self.update_output_table()
return
def handle_new_reaction(self) -> None:
rxnDia = ReactionDialog(self, targets=self.spanc.targets.keys(), extraParams=True)
rxnDia.new_reaction.connect(self.spanc.add_reaction)
if rxnDia.exec():
self.update_reaction_table()
return
def handle_update_reaction(self, row: int, col: int) -> None:
rxnName = self.reactionTable.verticalHeaderItem(row).text()
rxnDia = ReactionDialog(self, targets=self.spanc.targets.keys(), rxn=self.spanc.reactions[rxnName], rxnKey=rxnName, extraParams=True)
rxnDia.update_reaction.connect(self.spanc.update_reaction_parameters)
if rxnDia.exec():
self.update_reaction_table()
self.spanc.calculate_calibrations()
self.update_calibration_table()
self.spanc.calculate_outputs()
self.update_output_table()
return
def handle_new_calibration(self) -> None:
calDia = PeakDialog(PeakType.CALIBRATION, self.spanc.reactions.keys(), self)
calDia.new_peak.connect(self.spanc.add_calibration)
if calDia.exec():
self.update_calibration_table()
return
def handle_update_calibration(self, row: int, col: int) -> None:
peakData = self.spanc.calibrations[row]
calDia = PeakDialog(PeakType.CALIBRATION, self.spanc.reactions.keys(), self, peak=peakData)
calDia.new_peak.connect(self.spanc.add_calibration)
if calDia.exec():
self.update_calibration_table()
return
def handle_new_output(self) -> None:
outDia = PeakDialog(PeakType.OUTPUT, self.spanc.reactions.keys(), self)
outDia.new_peak.connect(self.spanc.add_output)
if outDia.exec():
self.update_output_table()
return
def handle_update_output(self, row: int, col: int) -> None:
peakData = self.spanc.calibrations[row]
outDia = PeakDialog(PeakType.OUTPUT, self.spanc.reactions.keys(), self, peak=peakData)
outDia.new_peak.connect(self.spanc.add_output)
if outDia.exec():
self.update_output_table()
return
def handle_change_fit_order(self, order: int) -> None:
self.spanc.set_fit_order(order)
def handle_run_fit(self) -> None:
order = self.spanc.fitter.polynomialOrder
npoints = len(self.spanc.calibrations)
if npoints < (order + 2):
print(f"Warning! Attempting to fit {npoints} data points with order {order} polyomial, too few degrees of freedom!")
print(f"Increase number of data points to at minimum {order+2} to use a polynomial of this order.")
return
fitData = self.spanc.fit()
xArray, yArray, xErrArray, yErrArray = convert_fit_points_to_arrays(fitData)
xMin = np.amin(xArray)
xMax = np.amax(xArray)
fitArray = np.linspace(xMin, xMax, 1000)
self.fitCanvas.axes.cla()
self.fitCanvas.axes.errorbar(xArray, yArray, yerr=yErrArray, xerr=xErrArray, marker="o", linestyle="None", elinewidth=2.0)
self.fitCanvas.axes.plot(fitArray, self.spanc.fitter.evaluate(fitArray))
self.fitCanvas.axes.set_xlabel(r"$x$ (mm)")
self.fitCanvas.axes.set_ylabel(r"$\rho$ (cm)")
self.fitCanvas.fig.tight_layout()
self.fitCanvas.draw()
self.spanc.calculate_outputs()
self.update_output_table()
self.fitFlag = True
residData = self.spanc.get_residuals()
xArray, residArray, studentResidArray = convert_resid_points_to_arrays(residData)
self.residCanvas.axes.cla()
self.residCanvas.axes.plot(xArray, residArray, marker="o", linestyle="None")
self.residCanvas.axes.hlines(0.0, xMin, xMax, colors="r", linestyles="dashed")
self.residCanvas.axes.set_xlabel(r"$x$ (mm)")
self.residCanvas.axes.set_ylabel(r"Residual (cm)")
self.residCanvas.fig.tight_layout()
self.residCanvas.draw()
self.update_fit_text(residArray, studentResidArray)
def update_target_table(self) -> None:
self.targetTable.setRowCount(len(self.spanc.targets))
self.targetTable.setVerticalHeaderLabels(self.spanc.targets.keys())
for row, key in enumerate(self.spanc.targets):
for i, layer in enumerate(self.spanc.targets[key].layer_details) :
self.targetTable.setItem(row, 0+i*2, QTableWidgetItem(str(layer.thickness)))
self.targetTable.setCellWidget(row, 1+i*2, QLabel(str(layer)))
self.targetTable.resizeColumnsToContents()
self.targetTable.resizeRowsToContents()
def update_reaction_table(self) -> None:
self.reactionTable.setRowCount(len(self.spanc.reactions))
self.reactionTable.setVerticalHeaderLabels(self.spanc.reactions.keys())
for row, rxn in enumerate(self.spanc.reactions.values()):
self.reactionTable.setItem(row, 0, QTableWidgetItem(str(rxn.targetMaterial)))
self.reactionTable.setCellWidget(row, 1, QLabel(str(rxn)))
self.reactionTable.setItem(row, 2, QTableWidgetItem(str(rxn.params.beamEnergy)))
self.reactionTable.setItem(row, 3, QTableWidgetItem(str(rxn.params.magneticField)))
self.reactionTable.setItem(row, 4, QTableWidgetItem(str(rxn.params.spsAngle)))
self.reactionTable.resizeColumnsToContents()
self.reactionTable.resizeRowsToContents()
def update_calibration_table(self) -> None:
self.calibrationTable.setRowCount(len(self.spanc.calibrations))
self.calibrationTable.setVerticalHeaderLabels(self.spanc.calibrations.keys())
for row, peak in enumerate(self.spanc.calibrations.values()):
self.calibrationTable.setItem(row, 0, QTableWidgetItem(peak.rxnName))
self.calibrationTable.setItem(row, 1, QTableWidgetItem(str(peak.position)))
self.calibrationTable.setItem(row, 2, QTableWidgetItem(str(peak.positionErrStat)))
self.calibrationTable.setItem(row, 3, QTableWidgetItem(str(peak.positionErrSys)))
self.calibrationTable.setItem(row, 4, QTableWidgetItem(str(peak.rho)))
self.calibrationTable.setItem(row, 5, QTableWidgetItem(str(peak.rhoErr)))
self.calibrationTable.setItem(row, 6, QTableWidgetItem(str(peak.excitation)))
self.calibrationTable.setItem(row, 7, QTableWidgetItem(str(peak.excitationErr)))
self.calibrationTable.resizeColumnsToContents()
self.calibrationTable.resizeRowsToContents()
def update_output_table(self) -> None:
self.outputTable.setRowCount(len(self.spanc.outputs))
self.outputTable.setVerticalHeaderLabels(self.spanc.outputs.keys())
for row, peak in enumerate(self.spanc.outputs.values()):
self.outputTable.setItem(row, 0, QTableWidgetItem(peak.rxnName))
self.outputTable.setItem(row, 1, QTableWidgetItem(str(peak.position)))
self.outputTable.setItem(row, 2, QTableWidgetItem(str(peak.positionErrStat)))
self.outputTable.setItem(row, 3, QTableWidgetItem(str(peak.positionErrSys)))
self.outputTable.setItem(row, 4, QTableWidgetItem(str(peak.rho)))
self.outputTable.setItem(row, 5, QTableWidgetItem(str(peak.rhoErr)))
self.outputTable.setItem(row, 6, QTableWidgetItem(str(peak.excitation)))
self.outputTable.setItem(row, 7, QTableWidgetItem(str(peak.excitationErr)))
self.outputTable.setItem(row, 8, QTableWidgetItem(str(peak.positionFWHM)))
self.outputTable.setItem(row, 9, QTableWidgetItem(str(peak.positionFWHMErr)))
self.outputTable.setItem(row, 10, QTableWidgetItem(str(peak.excitationFWHM)))
self.outputTable.setItem(row, 11, QTableWidgetItem(str(peak.excitationFWHMErr)))
self.outputTable.resizeColumnsToContents()
self.outputTable.resizeRowsToContents()
def update_fit_order(self) -> None:
self.fitOrderBox.setValue(self.spanc.fitter.polynomialOrder)
#generate markdown text string and render in the text edit
def update_fit_text(self, residuals: NDArray[np.float64], studentizedResiduals: NDArray[np.float64]) -> None:
markdownString = (f"# Fit Results\n"
f"## Polynomial Order: {self.spanc.fitter.polynomialOrder} \n \n"
f"## Chi-Square: {self.spanc.fitter.get_chisquare():.3f} \n \n"
f"## NDF: {self.spanc.fitter.get_ndf()} \n \n"
f"## Reduced Chi-Square {self.spanc.fitter.get_reduced_chisquare():.3f} \n \n"
f"## Parameter Values (a0 -> aN): {np.array_str(self.spanc.fitter.get_parameters(), precision=3)} \n \n"
f"## Parameter Uncertanties (ua0 -> uaN): {np.array_str(self.spanc.fitter.get_parameter_errors(), precision=3)} \n \n"
f"## Residuals (x0 -> xN): {np.array_str(residuals, precision=3)} \n \n"
f"## Studentized Residuals (x0 -> xN): {np.array_str(studentizedResiduals, precision=3)} \n \n")
self.fitResultText.setMarkdown(markdownString)
def run_spanc_ui() :
mpl.use("Qt5Agg")
app = QApplication.instance()
if not app:
app = QApplication(sys.argv)
app.setStyleSheet(load_stylesheet())
window = SpancGUI()
sys.exit(app.exec_())

0
spspy/__init__.py Normal file
View File

77
spspy/data/NuclearData.py Normal file
View File

@ -0,0 +1,77 @@
import numpy as np
import requests as req
import lxml.html as xhtml
from dataclasses import dataclass
PATH_TO_MASSFILE = "./etc/amdc2016_mass.txt"
@dataclass
class NucleusData:
mass: float = 0.0
elementSymbol: str = ""
isotopicSymbol: str = ""
prettyIsotopicSymbol: str = ""
Z: int = 0
A: int = 0
def __str__(self):
return self.isotopicSymbol
def get_latex_rep(self):
return "$^{" + str(self.A) + "}$" + self.elementSymbol
def generate_nucleus_id(z: np.uint32, a: np.uint32) -> np.uint32 :
return z*z + z + a if z > a else a*a + z
class NuclearDataMap:
U2MEV: float = 931.4940954
ELECTRON_MASS: float = 0.000548579909
def __init__(self):
self.map = {}
with open(PATH_TO_MASSFILE) as massfile:
massfile.readline()
massfile.readline()
for line in massfile:
entries = line.split()
data = NucleusData()
data.Z = int(entries[1])
data.A = int(entries[2])
data.mass = (float(entries[4]) + 1.0e-6 * float(entries[5]) - float(data.Z) * self.ELECTRON_MASS) * self.U2MEV
data.elementSymbol = entries[3]
data.isotopicSymbol = f"{data.A}{entries[3]}"
data.prettyIsotopicSymbol = f"<sup>{data.A}</sup>{entries[3]}"
self.map[generate_nucleus_id(data.Z, data.A)] = data
def get_data(self, z: np.uint32, a: np.uint32) -> NucleusData:
return self.map[generate_nucleus_id(z, a)]
global_nuclear_data = NuclearDataMap()
def get_excitations(Z: int, A: int) -> list[float]:
levels = []
text = ''
symbol = global_nuclear_data.get_data(Z, A).isotopicSymbol
site = req.get(f"https://www.nndc.bnl.gov/nudat2/getdatasetClassic.jsp?nucleus={symbol}&unc=nds")
contents = xhtml.fromstring(site.content)
tables = contents.xpath("//table")
rows = tables[2].xpath("./tr")
for row in rows[1:-2]:
entries = row.xpath("./td")
if len(entries) != 0:
entry = entries[0]
data = entry.xpath("./a")
if len(data) == 0:
text = entry.text
else:
text = data[0].text
text = text.replace('?', '')
text = text.replace('\xa0\xa0≈','')
levels.append(float(text)/1000.0) #convert to MeV
return levels

0
spspy/data/__init__.py Normal file
View File

42
spspy/ui/ExcitationDialog.py Normal file
View File

@ -0,0 +1,42 @@
from PySide6.QtWidgets import QLabel
from PySide6.QtWidgets import QVBoxLayout
from PySide6.QtWidgets import QDoubleSpinBox, QListWidget, QListWidgetItem
from PySide6.QtWidgets import QDialog, QDialogButtonBox
from PySide6.QtWidgets import QAbstractItemView
from PySide6.QtCore import Signal
class ExcitationDialog(QDialog):
new_level = Signal(str,float)
def __init__(self, parent, rxnList: list[str]) :
super().__init__(parent)
self.setWindowTitle("Add an Excitation")
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
self.buttonBox.accepted.connect(self.send_level)
self.buttonBox.rejected.connect(self.reject)
self.layout = QVBoxLayout()
self.setLayout(self.layout)
rxnLabel = QLabel("Choose a reaction",self)
self.reactionList = QListWidget(self)
self.reactionList.setSelectionMode(QAbstractItemView.SelectionMode.SingleSelection)
for rxnName in rxnList:
item = QListWidgetItem(rxnName)
item.setText("")
self.reactionList.addItem(item)
self.reactionList.setItemWidget(item, QLabel(rxnName))
stateLabel = QLabel("New state energy",self)
self.stateInput = QDoubleSpinBox(self)
self.stateInput.setRange(0.0,40.0)
self.stateInput.setSuffix(" MeV")
self.layout.addWidget(rxnLabel)
self.layout.addWidget(self.reactionList)
self.layout.addWidget(stateLabel)
self.layout.addWidget(self.stateInput)
self.layout.addWidget(self.buttonBox)
def send_level(self) -> None:
items = self.reactionList.selectedItems()
if len(items) == 1:
self.new_level.emit(self.reactionList.itemWidget(items[0]).text(),self.stateInput.value())

18
spspy/ui/MPLCanvas.py Normal file
View File

@ -0,0 +1,18 @@
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg
from matplotlib.figure import Figure
import matplotlib
matplotlib.rcParams["axes.facecolor"] = "#0F1016"
matplotlib.rcParams["figure.facecolor"] = "#0F1016"
matplotlib.rcParams["axes.labelcolor"] = "w"
matplotlib.rcParams["axes.edgecolor"] = "w"
matplotlib.rcParams["xtick.color"] = "w"
matplotlib.rcParams["ytick.color"] = "w"
matplotlib.rcParams["text.color"] = "w"
class MPLCanvas(FigureCanvasQTAgg):
def __init__(self, parent=None, width=3, height=4, dpi=100):
self.fig = Figure(figsize=(width, height), dpi=dpi, edgecolor="black",linewidth=0.5)
self.axes = self.fig.add_subplot(111)
self.axes.spines['top'].set_visible(False)
super(MPLCanvas, self).__init__(self.fig)

126
spspy/ui/PeakDialog.py Normal file
View File

@ -0,0 +1,126 @@
from ..Spanc import Peak, PeakType, INVALID_PEAK_ID
from PySide6.QtWidgets import QLabel
from PySide6.QtWidgets import QVBoxLayout, QFormLayout, QGroupBox
from PySide6.QtWidgets import QComboBox, QDoubleSpinBox
from PySide6.QtWidgets import QDialog, QDialogButtonBox
from PySide6.QtCore import Signal
class PeakDialog(QDialog):
new_peak = Signal(Peak)
def __init__(self, peakType: PeakType, rxnList: list[str], parent=None, peak: Peak=None):
super().__init__(parent)
self.setWindowTitle(f"Add A {peakType.value} Peak")
rnameLabel = QLabel("Reaction Name", self)
self.rxnNameBox = QComboBox(self)
for reaction in rxnList:
self.rxnNameBox.addItem(reaction)
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
self.buttonBox.rejected.connect(self.reject)
self.centralLayout = QVBoxLayout()
self.setLayout(self.centralLayout)
self.centralLayout.addWidget(rnameLabel)
self.centralLayout.addWidget(self.rxnNameBox)
if peakType == PeakType.CALIBRATION:
self.create_calibration_inputs()
if peak is not None:
self.set_calibration_inputs(peak)
self.peakID = peak.peakID
self.buttonBox.accepted.connect(self.send_update_calibration_peak)
else:
self.buttonBox.accepted.connect(self.send_calibration_peak)
elif peakType == PeakType.OUTPUT:
self.create_output_inputs()
if peak is not None:
self.set_output_inputs(peak)
self.peakID = peak.peakID
self.buttonBox.accepted.connect(self.send_update_output_peak)
else:
self.buttonBox.accepted.connect(self.send_output_peak)
self.centralLayout.addWidget(self.buttonBox)
def create_calibration_inputs(self) -> None:
self.inputGroupBox = QGroupBox("Peak Parameters",self)
inputLayout = QFormLayout()
self.xInput = QDoubleSpinBox(self.inputGroupBox)
self.xInput.setRange(-999, 999)
self.xInput.setDecimals(6)
self.uxsysInput = QDoubleSpinBox(self.inputGroupBox)
self.uxsysInput.setRange(-999, 999)
self.uxsysInput.setDecimals(6)
self.uxstatInput = QDoubleSpinBox(self.inputGroupBox)
self.uxstatInput.setRange(-999, 999)
self.uxstatInput.setDecimals(6)
self.exInput = QDoubleSpinBox(self.inputGroupBox)
self.exInput.setDecimals(6)
self.uexInput = QDoubleSpinBox(self.inputGroupBox)
self.uexInput.setDecimals(6)
inputLayout.addRow("Position(mm)", self.xInput)
inputLayout.addRow("Position Stat. Error(mm)", self.uxstatInput)
inputLayout.addRow("Position Sys. Error(mm)", self.uxsysInput)
inputLayout.addRow("Excitation Energy(MeV)", self.exInput)
inputLayout.addRow("Excitation Energy Error(MeV)", self.uexInput)
self.inputGroupBox.setLayout(inputLayout)
self.centralLayout.addWidget(self.inputGroupBox)
def create_output_inputs(self) -> None:
self.inputGroupBox = QGroupBox("Peak Parameters",self)
inputLayout = QFormLayout()
self.xInput = QDoubleSpinBox(self.inputGroupBox)
self.xInput.setRange(-999, 999)
self.xInput.setDecimals(6)
self.uxsysInput = QDoubleSpinBox(self.inputGroupBox)
self.uxsysInput.setRange(-999, 999)
self.uxsysInput.setDecimals(6)
self.uxstatInput = QDoubleSpinBox(self.inputGroupBox)
self.uxstatInput.setRange(-999, 999)
self.uxstatInput.setDecimals(6)
self.fwhmInput = QDoubleSpinBox(self.inputGroupBox)
self.fwhmInput.setDecimals(6)
self.ufwhmInput = QDoubleSpinBox(self.inputGroupBox)
self.ufwhmInput.setDecimals(6)
inputLayout.addRow("Position(mm)", self.xInput)
inputLayout.addRow("Position Stat. Error(mm)", self.uxstatInput)
inputLayout.addRow("Position Sys. Error(mm)", self.uxsysInput)
inputLayout.addRow("Position FWHM(mm)", self.fwhmInput)
inputLayout.addRow("Position FWHM Error(mm)", self.ufwhmInput)
self.inputGroupBox.setLayout(inputLayout)
self.centralLayout.addWidget(self.inputGroupBox)
def set_calibration_inputs(self, peak: Peak) -> None:
self.rxnNameBox.setCurrentIndex(self.rxnNameBox.findText(peak.rxnName))
self.xInput.setValue(peak.position)
self.uxsysInput.setValue(peak.positionErrSys)
self.uxstatInput.setValue(peak.positionErrStat)
self.exInput.setValue(peak.excitation)
self.uexInput.setValue(peak.excitationErr)
def set_output_inputs(self, peak: Peak) -> None:
self.rxnNameBox.setCurrentIndex(self.rxnNameBox.findText(peak.rxnName))
self.xInput.setValue(peak.position)
self.uxsysInput.setValue(peak.positionErrSys)
self.uxstatInput.setValue(peak.positionErrStat)
self.fwhmInput.setValue(peak.positionFWHM)
self.ufwhmInput.setValue(peak.positionFWHMErr)
def send_calibration_peak(self) -> None:
peak = Peak(excitation=self.exInput.value(), excitationErr=self.uexInput.value(), position=self.xInput.value(),
positionErrStat=self.uxstatInput.value(), positionErrSys=self.uxsysInput.value(), rxnName=self.rxnNameBox.currentText(), peakID=INVALID_PEAK_ID)
self.new_peak.emit(peak)
def send_update_calibration_peak(self) -> None:
peak = Peak(excitation=self.exInput.value(), excitationErr=self.uexInput.value(), position=self.xInput.value(),
positionErrStat=self.uxstatInput.value(), positionErrSys=self.uxsysInput.value(), rxnName=self.rxnNameBox.currentText(), peakID=self.peakID)
self.new_peak.emit(peak)
def send_output_peak(self) -> None:
peak = Peak(position=self.xInput.value(), positionErrStat=self.uxstatInput.value(), positionErrSys=self.uxsysInput.value(),
positionFWHM=self.fwhmInput.value(), positionFWHMErr=self.ufwhmInput.value(), rxnName=self.rxnNameBox.currentText(), peakID=INVALID_PEAK_ID)
self.new_peak.emit(peak)
def send_update_output_peak(self) -> None:
peak = Peak(position=self.xInput.value(), positionErrStat=self.uxstatInput.value(), positionErrSys=self.uxsysInput.value(),
positionFWHM=self.fwhmInput.value(), positionFWHMErr=self.ufwhmInput.value(), rxnName=self.rxnNameBox.currentText(), peakID=self.peakID)
self.new_peak.emit(peak)

121
spspy/ui/ReactionDialog.py Normal file
View File

@ -0,0 +1,121 @@
from ..SPSReaction import RxnParameters, Reaction, create_reaction_parameters
from ..data.NuclearData import global_nuclear_data
from PySide6.QtWidgets import QLabel
from PySide6.QtWidgets import QVBoxLayout, QFormLayout, QGroupBox
from PySide6.QtWidgets import QSpinBox, QDoubleSpinBox, QComboBox
from PySide6.QtWidgets import QDialog, QDialogButtonBox
from PySide6.QtCore import Signal
MAXIMUM_NUCLEAR_Z: int = 110
MAXIMUM_NUCLEAR_A: int = 270
MINIMUM_BEAM_ENERGY: float = 0.0 #MeV
MAXIMUM_BEAM_ENERGY: float = 100.0 #MeV
MINIMUM_SPS_ANGLE: float = 0.0 #deg
MAXIMUM_SPS_ANGLE: float = 90.0 #deg
MINIMUM_MAG_FIELD: float = 0.0 #kG
MAXIMUM_MAG_FIELD: float = 16.0 #kG
class ReactionDialog(QDialog):
new_reaction = Signal(RxnParameters, str)
update_reaction = Signal(float, float, float, str)
def __init__(self, parent=None, rxn: Reaction =None, rxnKey: str =None, targets: list[str]=None, extraParams: bool = False):
super().__init__(parent)
self.setWindowTitle("Add Reaction")
self.extraParams = extraParams
tnameLabel = QLabel("Target Name", self)
self.targetNameBox = QComboBox(self)
for target in targets:
self.targetNameBox.addItem(target)
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
if rxn is not None:
self.buttonBox.accepted.connect(self.send_reaction_update)
else:
self.buttonBox.accepted.connect(self.send_reaction)
self.buttonBox.rejected.connect(self.reject)
self.layout = QVBoxLayout()
self.setLayout(self.layout)
self.layout.addWidget(tnameLabel)
self.layout.addWidget(self.targetNameBox)
self.create_reaction_inputs()
if rxn is not None:
self.set_initial_values(rxn)
self.rxnKey = rxnKey
self.layout.addWidget(self.buttonBox)
def send_reaction(self) -> None :
params = create_reaction_parameters(self.ztInput.value(), self.atInput.value(), self.zpInput.value(), self.apInput.value(), self.zeInput.value(), self.aeInput.value())
if self.extraParams == True:
params.beamEnergy = self.bkeInput.value()
params.spsAngle = self.thetaInput.value()
params.magneticField = self.bfieldInput.value()
self.new_reaction.emit(params, self.targetNameBox.currentText())
def send_reaction_update(self) -> None:
self.update_reaction.emit(self.bkeInput.value(), self.thetaInput.value(), self.bfieldInput.value(), self.rxnKey)
def create_reaction_inputs(self) -> None:
self.nucleiGroupBox = QGroupBox("Reaction Nuclei",self)
inputLayout = QFormLayout()
self.ztInput = QSpinBox(self.nucleiGroupBox)
self.ztInput.setRange(1, MAXIMUM_NUCLEAR_Z)
self.atInput = QSpinBox(self.nucleiGroupBox)
self.atInput.setRange(1, MAXIMUM_NUCLEAR_A)
self.zpInput = QSpinBox(self.nucleiGroupBox)
self.zpInput.setRange(1, MAXIMUM_NUCLEAR_Z)
self.apInput = QSpinBox(self.nucleiGroupBox)
self.apInput.setRange(1, MAXIMUM_NUCLEAR_A)
self.zeInput = QSpinBox(self.nucleiGroupBox)
self.zeInput.setRange(1, MAXIMUM_NUCLEAR_Z)
self.aeInput = QSpinBox(self.nucleiGroupBox)
self.aeInput.setRange(1, MAXIMUM_NUCLEAR_A)
inputLayout.addRow("ZT",self.ztInput)
inputLayout.addRow("AT",self.atInput)
inputLayout.addRow("ZP",self.zpInput)
inputLayout.addRow("AP",self.apInput)
inputLayout.addRow("ZE",self.zeInput)
inputLayout.addRow("AE",self.aeInput)
self.nucleiGroupBox.setLayout(inputLayout)
self.layout.addWidget(self.nucleiGroupBox)
if self.extraParams == True:
self.parameterGroupBox = QGroupBox("Reaction Parameters", self)
parameterLayout = QFormLayout()
self.bkeInput = QDoubleSpinBox(self.parameterGroupBox)
self.bkeInput.setRange(MINIMUM_BEAM_ENERGY, MAXIMUM_BEAM_ENERGY)
self.bkeInput.setDecimals(4)
self.thetaInput = QDoubleSpinBox(self.parameterGroupBox)
self.thetaInput.setRange(MINIMUM_SPS_ANGLE, MAXIMUM_SPS_ANGLE)
self.thetaInput.setDecimals(4)
self.bfieldInput = QDoubleSpinBox(self.parameterGroupBox)
self.bfieldInput.setRange(MINIMUM_MAG_FIELD, MAXIMUM_MAG_FIELD)
self.bfieldInput.setDecimals(6)
parameterLayout.addRow(QLabel("E<sub>beam</sub>(Mev)"),self.bkeInput)
parameterLayout.addRow("<p>&theta;<sub>SPS</sub>(deg)</p>",self.thetaInput)
parameterLayout.addRow("B(kG)",self.bfieldInput)
self.parameterGroupBox.setLayout(parameterLayout)
self.layout.addWidget(self.parameterGroupBox)
def set_initial_values(self, rxn: Reaction) -> None:
self.targetNameBox.setCurrentIndex(self.targetNameBox.findText(rxn.targetMaterial.name))
self.targetNameBox.setEnabled(False)
self.ztInput.setValue(rxn.params.target.Z)
self.ztInput.setEnabled(False)
self.atInput.setValue(rxn.params.target.A)
self.atInput.setEnabled(False)
self.zpInput.setValue(rxn.params.projectile.Z)
self.zpInput.setEnabled(False)
self.apInput.setValue(rxn.params.projectile.A)
self.apInput.setEnabled(False)
self.zeInput.setValue(rxn.params.ejectile.Z)
self.zeInput.setEnabled(False)
self.aeInput.setValue(rxn.params.ejectile.A)
self.aeInput.setEnabled(False)
self.bkeInput.setValue(rxn.params.beamEnergy)
self.thetaInput.setValue(rxn.params.spsAngle)
self.bfieldInput.setValue(rxn.params.magneticField)

150
spspy/ui/TargetDialog.py Normal file
View File

@ -0,0 +1,150 @@
from ..SPSTarget import SPSTarget, TargetLayer
from PySide6.QtWidgets import QLabel, QLineEdit
from PySide6.QtWidgets import QHBoxLayout, QVBoxLayout, QGridLayout, QGroupBox
from PySide6.QtWidgets import QSpinBox, QDoubleSpinBox
from PySide6.QtWidgets import QDialog, QDialogButtonBox
from PySide6.QtCore import Signal
class TargetDialog(QDialog):
new_target = Signal(str, list)
def __init__(self, parent=None, target: SPSTarget =None):
super().__init__(parent)
self.setWindowTitle("Add A Target")
nameLabel = QLabel("Target name", self)
self.nameInput = QLineEdit(self)
QBtn = QDialogButtonBox.Ok | QDialogButtonBox.Cancel
self.buttonBox = QDialogButtonBox(QBtn)
self.buttonBox.accepted.connect(self.accept)
self.buttonBox.accepted.connect(self.send_target)
self.buttonBox.rejected.connect(self.reject)
self.layout = QVBoxLayout()
self.setLayout(self.layout)
self.layout.addWidget(nameLabel)
self.layout.addWidget(self.nameInput)
self.create_target_inputs()
if target is not None:
self.set_initial_values(target)
self.layout.addWidget(self.buttonBox)
def create_target_inputs(self) -> None:
self.layerGroup = QGroupBox("Target Layers", self)
layerGroupLayout = QHBoxLayout()
self.layerAInputs = []
self.layerZInputs = []
self.layerSInputs = []
self.layerThickInputs = []
self.layer1GroupBox = QGroupBox("Layer 1", self.layerGroup)
layer1Layout = QVBoxLayout()
thick1Label = QLabel("Thickness(ug/cm^2)", self.layer1GroupBox)
self.layerThickInputs.append(QDoubleSpinBox(self.layer1GroupBox))
self.layerThickInputs[0].setRange(0, 999.0)
self.layerThickInputs[0].setDecimals(4)
self.layer1ComponentsBox = QGroupBox("Layer 1 Components", self.layer1GroupBox)
layer1compLayout = QGridLayout()
layer1compLayout.addWidget(QLabel("Z", self.layer1ComponentsBox), 0, 1)
layer1compLayout.addWidget(QLabel("A", self.layer1ComponentsBox), 0, 2)
layer1compLayout.addWidget(QLabel("Stoich", self.layer1ComponentsBox), 0, 3)
for i in range(3):
layer1compLayout.addWidget(QLabel("Component"+str(i), self.layer1ComponentsBox), i+1, 0)
self.layerZInputs.append(QSpinBox(self.layer1ComponentsBox))
self.layerAInputs.append(QSpinBox(self.layer1ComponentsBox))
self.layerSInputs.append(QSpinBox(self.layer1ComponentsBox))
layer1compLayout.addWidget(self.layerZInputs[i], i+1, 1)
layer1compLayout.addWidget(self.layerAInputs[i], i+1, 2)
layer1compLayout.addWidget(self.layerSInputs[i], i+1, 3)
self.layer1ComponentsBox.setLayout(layer1compLayout)
layer1Layout.addWidget(thick1Label)
layer1Layout.addWidget(self.layerThickInputs[0])
layer1Layout.addWidget(self.layer1ComponentsBox)
self.layer1GroupBox.setLayout(layer1Layout)
self.layer2GroupBox = QGroupBox("Layer 2", self.layerGroup)
layer2Layout = QVBoxLayout()
thick2Label = QLabel("Thickness(ug/cm^2)", self.layer2GroupBox)
self.layerThickInputs.append(QDoubleSpinBox(self.layer2GroupBox))
self.layerThickInputs[1].setRange(0, 999.0)
self.layerThickInputs[1].setDecimals(4)
self.layer2ComponentsBox = QGroupBox("Layer 2 Components", self.layer2GroupBox)
layer2compLayout = QGridLayout()
layer2compLayout.addWidget(QLabel("Z", self.layer2ComponentsBox), 0, 1)
layer2compLayout.addWidget(QLabel("A", self.layer2ComponentsBox), 0, 2)
layer2compLayout.addWidget(QLabel("Stoich", self.layer2ComponentsBox), 0, 3)
for i in range(3):
layer2compLayout.addWidget(QLabel("Component"+str(i), self.layer2ComponentsBox), i+1, 0)
self.layerZInputs.append(QSpinBox(self.layer2ComponentsBox))
self.layerAInputs.append(QSpinBox(self.layer2ComponentsBox))
self.layerSInputs.append(QSpinBox(self.layer2ComponentsBox))
layer2compLayout.addWidget(self.layerZInputs[i+3], i+1, 1)
layer2compLayout.addWidget(self.layerAInputs[i+3], i+1, 2)
layer2compLayout.addWidget(self.layerSInputs[i+3], i+1, 3)
self.layer2ComponentsBox.setLayout(layer2compLayout)
layer2Layout.addWidget(thick2Label)
layer2Layout.addWidget(self.layerThickInputs[1])
layer2Layout.addWidget(self.layer2ComponentsBox)
self.layer2GroupBox.setLayout(layer2Layout)
self.layer3GroupBox = QGroupBox("Layer 3", self.layerGroup)
layer3Layout = QVBoxLayout()
thick3Label = QLabel("Thickness(ug/cm^2)", self.layer3GroupBox)
self.layerThickInputs.append(QDoubleSpinBox(self.layer3GroupBox))
self.layerThickInputs[2].setRange(0, 999.0)
self.layerThickInputs[2].setDecimals(4)
self.layer3ComponentsBox = QGroupBox("Layer 3 Components", self.layer3GroupBox)
layer3compLayout = QGridLayout()
layer3compLayout.addWidget(QLabel("Z", self.layer3ComponentsBox), 0, 1)
layer3compLayout.addWidget(QLabel("A", self.layer3ComponentsBox), 0, 2)
layer3compLayout.addWidget(QLabel("Stoich", self.layer3ComponentsBox), 0, 3)
for i in range(3):
layer3compLayout.addWidget(QLabel("Component"+str(i), self.layer3ComponentsBox), i+1, 0)
self.layerZInputs.append(QSpinBox(self.layer3ComponentsBox))
self.layerAInputs.append(QSpinBox(self.layer3ComponentsBox))
self.layerSInputs.append(QSpinBox(self.layer3ComponentsBox))
layer3compLayout.addWidget(self.layerZInputs[i+6], i+1, 1)
layer3compLayout.addWidget(self.layerAInputs[i+6], i+1, 2)
layer3compLayout.addWidget(self.layerSInputs[i+6], i+1, 3)
self.layer3ComponentsBox.setLayout(layer3compLayout)
layer3Layout.addWidget(thick3Label)
layer3Layout.addWidget(self.layerThickInputs[2])
layer3Layout.addWidget(self.layer3ComponentsBox)
self.layer3GroupBox.setLayout(layer3Layout)
layerGroupLayout.addWidget(self.layer1GroupBox)
layerGroupLayout.addWidget(self.layer2GroupBox)
layerGroupLayout.addWidget(self.layer3GroupBox)
self.layerGroup.setLayout(layerGroupLayout)
self.layout.addWidget(self.layerGroup)
def set_initial_values(self, target: SPSTarget) -> None:
self.nameInput.setText(target.name)
self.nameInput.setReadOnly(True)
for i, layer in enumerate(target.layer_details):
self.layerThickInputs[i].setValue(layer.thickness)
for j, (z, a, s) in enumerate(layer.compound_list):
self.layerZInputs[j+i*3].setValue(z)
self.layerAInputs[j+i*3].setValue(a)
self.layerSInputs[j+i*3].setValue(s)
def send_target(self) -> None:
name = self.nameInput.text()
if name == "":
return
tlist = []
z = 0
a = 0
s = 0
for i in range(3):
tempLayer = TargetLayer()
tempLayer.thickness = self.layerThickInputs[i].value()
for j in range(3):
z = self.layerZInputs[i*3 + j].value()
a = self.layerAInputs[i*3 + j].value()
s = self.layerSInputs[i*3 + j].value()
if z != 0 and a != 0 and s != 0:
tempLayer.compound_list.append((z, a, s))
if len(tempLayer.compound_list) != 0:
tlist.append(tempLayer)
if len(tlist) != 0:
self.new_target.emit(name, tlist)

8
test_config.txt
View File

@ -1,8 +0,0 @@
BeamEnergy(Mev): 16.0
B-field(kG): 7.8
Angle(deg): 27.5
RhoMin: 69.5 RhoMax: 83.5
ZT AT ZP AP ZE AE
6 12 1 2 1 1
6 12 1 2 1 2
8 16 1 2 1 1