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

Overhaul of basically everything. Added custom math and completed SABRE efficiency

Browse Source
This commit is contained in:
Gordon McCann 2020-12-09 17:07:58 -05:00
parent b009f56ba0
commit 3aac64acb4
43 changed files with 581 additions and 1670 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
.DS_Store vendored
View File

Binary file not shown.

8
.gitignore vendored
View File

@ -1,8 +1,16 @@
###Files to ignore### ###Files to ignore###
*.root *.root
*.png *.png
*.jpg
*.sublime-workspace *.sublime-workspace
*.sublime-project *.sublime-project
*.C
*.cxx
*.pcm
*.o
*.swp
.DS_Store
###Keep this file### ###Keep this file###
!.gitignore !.gitignore

3
bin/.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
###keep only the directory, not the contents###
*
!.gitignore

BIN
bin/kinematics
View File

Binary file not shown.

BIN
bin/kinematics_dict_rdict.pcm
View File

Binary file not shown.

484
checkfile.txt
View File

@ -1,484 +0,0 @@
RINGS
det ch corner x_fl y_fl z_fl r_fl t_fl p_fl x_ti y_ti z_ti r_ti t_ti p_ti
0 0 0 1.78297 3.46928 0 3.90062 90 62.8 3.00456 1.10206 -10.22 10.7093 162.613 20.1428
0 0 1 -1.78297 3.46928 0 3.90062 90 117.2 0.11966 3.19806 -10.22 10.7093 162.613 87.8572
0 0 2 -1.49014 2.8995 0 3.26 90 117.2 0.100008 2.67282 -10.5862 10.9189 165.821 87.8572
0 0 3 1.49014 2.8995 0 3.26 90 62.8 2.5111 0.921061 -10.5862 10.9189 165.821 20.1428
0 1 0 2.0758 4.03906 0 4.54125 90 62.8 3.49802 1.28306 -9.85374 10.5346 159.287 20.1428
0 1 1 -2.0758 4.03906 0 4.54125 90 117.2 0.139313 3.7233 -9.85374 10.5346 159.287 87.8572
0 1 2 -1.78297 3.46928 0 3.90062 90 117.2 0.11966 3.19806 -10.22 10.7093 162.613 87.8572
0 1 3 1.78297 3.46928 0 3.90062 90 62.8 3.00456 1.10206 -10.22 10.7093 162.613 20.1428
0 2 0 2.36862 4.60884 0 5.18187 90 62.8 3.99148 1.46406 -9.48749 10.3965 155.862 20.1428
0 2 1 -2.36862 4.60884 0 5.18187 90 117.2 0.158965 4.24854 -9.48749 10.3965 155.862 87.8572
0 2 2 -2.0758 4.03906 0 4.54125 90 117.2 0.139313 3.7233 -9.85374 10.5346 159.287 87.8572
0 2 3 2.0758 4.03906 0 4.54125 90 62.8 3.49802 1.28306 -9.85374 10.5346 159.287 20.1428
0 3 0 2.66145 5.17863 0 5.8225 90 62.8 4.48494 1.64505 -9.12124 10.2965 152.357 20.1428
0 3 1 -2.66145 5.17863 0 5.8225 90 117.2 0.178618 4.77378 -9.12124 10.2965 152.357 87.8572
0 3 2 -2.36862 4.60884 0 5.18187 90 117.2 0.158965 4.24854 -9.48749 10.3965 155.862 87.8572
0 3 3 2.36862 4.60884 0 5.18187 90 62.8 3.99148 1.46406 -9.48749 10.3965 155.862 20.1428
0 4 0 2.95428 5.74841 0 6.46312 90 62.8 4.9784 1.82605 -8.75499 10.2357 148.798 20.1428
0 4 1 -2.95428 5.74841 0 6.46312 90 117.2 0.19827 5.29902 -8.75499 10.2357 148.798 87.8572
0 4 2 -2.66145 5.17863 0 5.8225 90 117.2 0.178618 4.77378 -9.12124 10.2965 152.357 87.8572
0 4 3 2.66145 5.17863 0 5.8225 90 62.8 4.48494 1.64505 -9.12124 10.2965 152.357 20.1428
0 5 0 3.24711 6.31819 0 7.10375 90 62.8 5.47186 2.00705 -8.38874 10.2147 145.209 20.1428
0 5 1 -3.24711 6.31819 0 7.10375 90 117.2 0.217923 5.82426 -8.38874 10.2147 145.209 87.8572
0 5 2 -2.95428 5.74841 0 6.46312 90 117.2 0.19827 5.29902 -8.75499 10.2357 148.798 87.8572
0 5 3 2.95428 5.74841 0 6.46312 90 62.8 4.9784 1.82605 -8.75499 10.2357 148.798 20.1428
0 6 0 3.53994 6.88797 0 7.74438 90 62.8 5.96532 2.18805 -8.0225 10.2339 141.62 20.1428
0 6 1 -3.53994 6.88797 0 7.74438 90 117.2 0.237576 6.3495 -8.0225 10.2339 141.62 87.8572
0 6 2 -3.24711 6.31819 0 7.10375 90 117.2 0.217923 5.82426 -8.38874 10.2147 145.209 87.8572
0 6 3 3.24711 6.31819 0 7.10375 90 62.8 5.47186 2.00705 -8.38874 10.2147 145.209 20.1428
0 7 0 3.83277 7.45776 0 8.385 90 62.8 6.45877 2.36905 -7.65625 10.293 138.059 20.1428
0 7 1 -3.83277 7.45776 0 8.385 90 117.2 0.257228 6.87474 -7.65625 10.293 138.059 87.8572
0 7 2 -3.53994 6.88797 0 7.74438 90 117.2 0.237576 6.3495 -8.0225 10.2339 141.62 87.8572
0 7 3 3.53994 6.88797 0 7.74438 90 62.8 5.96532 2.18805 -8.0225 10.2339 141.62 20.1428
0 8 0 4.12559 8.02754 0 9.02562 90 62.8 6.95223 2.55005 -7.29 10.3914 134.551 20.1428
0 8 1 -4.12559 8.02754 0 9.02562 90 117.2 0.276881 7.39997 -7.29 10.3914 134.551 87.8572
0 8 2 -3.83277 7.45776 0 8.385 90 117.2 0.257228 6.87474 -7.65625 10.293 138.059 87.8572
0 8 3 3.83277 7.45776 0 8.385 90 62.8 6.45877 2.36905 -7.65625 10.293 138.059 20.1428
0 9 0 4.41842 8.59732 0 9.66625 90 62.8 7.44569 2.73105 -6.92375 10.5278 131.122 20.1428
0 9 1 -4.41842 8.59732 0 9.66625 90 117.2 0.296533 7.92521 -6.92375 10.5278 131.122 87.8572
0 9 2 -4.12559 8.02754 0 9.02562 90 117.2 0.276881 7.39997 -7.29 10.3914 134.551 87.8572
0 9 3 4.12559 8.02754 0 9.02562 90 62.8 6.95223 2.55005 -7.29 10.3914 134.551 20.1428
0 10 0 4.71125 9.1671 0 10.3069 90 62.8 7.93915 2.91204 -6.5575 10.701 127.792 20.1428
0 10 1 -4.71125 9.1671 0 10.3069 90 117.2 0.316186 8.45045 -6.5575 10.701 127.792 87.8572
0 10 2 -4.41842 8.59732 0 9.66625 90 117.2 0.296533 7.92521 -6.92375 10.5278 131.122 87.8572
0 10 3 4.41842 8.59732 0 9.66625 90 62.8 7.44569 2.73105 -6.92375 10.5278 131.122 20.1428
0 11 0 5.00408 9.73689 0 10.9475 90 62.8 8.43261 3.09304 -6.19125 10.9091 124.578 20.1428
0 11 1 -5.00408 9.73689 0 10.9475 90 117.2 0.335839 8.97569 -6.19125 10.9091 124.578 87.8572
0 11 2 -4.71125 9.1671 0 10.3069 90 117.2 0.316186 8.45045 -6.5575 10.701 127.792 87.8572
0 11 3 4.71125 9.1671 0 10.3069 90 62.8 7.93915 2.91204 -6.5575 10.701 127.792 20.1428
0 12 0 5.29691 10.3067 0 11.5881 90 62.8 8.92607 3.27404 -5.825 11.1501 121.495 20.1428
0 12 1 -5.29691 10.3067 0 11.5881 90 117.2 0.355491 9.50093 -5.825 11.1501 121.495 87.8572
0 12 2 -5.00408 9.73689 0 10.9475 90 117.2 0.335839 8.97569 -6.19125 10.9091 124.578 87.8572
0 12 3 5.00408 9.73689 0 10.9475 90 62.8 8.43261 3.09304 -6.19125 10.9091 124.578 20.1428
0 13 0 5.58974 10.8765 0 12.2287 90 62.8 9.41953 3.45504 -5.45875 11.422 118.549 20.1428
0 13 1 -5.58974 10.8765 0 12.2287 90 117.2 0.375144 10.0262 -5.45875 11.422 118.549 87.8572
0 13 2 -5.29691 10.3067 0 11.5881 90 117.2 0.355491 9.50093 -5.825 11.1501 121.495 87.8572
0 13 3 5.29691 10.3067 0 11.5881 90 62.8 8.92607 3.27404 -5.825 11.1501 121.495 20.1428
0 14 0 5.88256 11.4462 0 12.8694 90 62.8 9.91299 3.63604 -5.0925 11.7227 115.748 20.1428
0 14 1 -5.88256 11.4462 0 12.8694 90 117.2 0.394796 10.5514 -5.0925 11.7227 115.748 87.8572
0 14 2 -5.58974 10.8765 0 12.2287 90 117.2 0.375144 10.0262 -5.45875 11.422 118.549 87.8572
0 14 3 5.58974 10.8765 0 12.2287 90 62.8 9.41953 3.45504 -5.45875 11.422 118.549 20.1428
0 15 0 6.17539 12.016 0 13.51 90 62.8 10.4064 3.81704 -4.72625 12.05 113.093 20.1428
0 15 1 -6.17539 12.016 0 13.51 90 117.2 0.414449 11.0766 -4.72625 12.05 113.093 87.8572
0 15 2 -5.88256 11.4462 0 12.8694 90 117.2 0.394796 10.5514 -5.0925 11.7227 115.748 87.8572
0 15 3 5.88256 11.4462 0 12.8694 90 62.8 9.91299 3.63604 -5.0925 11.7227 115.748 20.1428
1 0 0 1.78297 3.46928 0 3.90062 90 62.8 1.97658 -2.51695 -10.22 10.7093 162.613 308.143
1 0 1 -1.78297 3.46928 0 3.90062 90 117.2 3.07852 0.874452 -10.22 10.7093 162.613 15.8572
1 0 2 -1.49014 2.8995 0 3.26 90 117.2 2.57291 0.730835 -10.5862 10.9189 165.821 15.8572
1 0 3 1.49014 2.8995 0 3.26 90 62.8 1.65195 -2.10358 -10.5862 10.9189 165.821 308.143
1 1 0 2.0758 4.03906 0 4.54125 90 62.8 2.30121 -2.93033 -9.85374 10.5346 159.287 308.143
1 1 1 -2.0758 4.03906 0 4.54125 90 117.2 3.58412 1.01807 -9.85374 10.5346 159.287 15.8572
1 1 2 -1.78297 3.46928 0 3.90062 90 117.2 3.07852 0.874452 -10.22 10.7093 162.613 15.8572
1 1 3 1.78297 3.46928 0 3.90062 90 62.8 1.97658 -2.51695 -10.22 10.7093 162.613 308.143
1 2 0 2.36862 4.60884 0 5.18187 90 62.8 2.62584 -3.3437 -9.48749 10.3965 155.862 308.143
1 2 1 -2.36862 4.60884 0 5.18187 90 117.2 4.08973 1.16169 -9.48749 10.3965 155.862 15.8572
1 2 2 -2.0758 4.03906 0 4.54125 90 117.2 3.58412 1.01807 -9.85374 10.5346 159.287 15.8572
1 2 3 2.0758 4.03906 0 4.54125 90 62.8 2.30121 -2.93033 -9.85374 10.5346 159.287 308.143
1 3 0 2.66145 5.17863 0 5.8225 90 62.8 2.95046 -3.75708 -9.12124 10.2965 152.357 308.143
1 3 1 -2.66145 5.17863 0 5.8225 90 117.2 4.59533 1.3053 -9.12124 10.2965 152.357 15.8572
1 3 2 -2.36862 4.60884 0 5.18187 90 117.2 4.08973 1.16169 -9.48749 10.3965 155.862 15.8572
1 3 3 2.36862 4.60884 0 5.18187 90 62.8 2.62584 -3.3437 -9.48749 10.3965 155.862 308.143
1 4 0 2.95428 5.74841 0 6.46312 90 62.8 3.27509 -4.17046 -8.75499 10.2357 148.798 308.143
1 4 1 -2.95428 5.74841 0 6.46312 90 117.2 5.10094 1.44892 -8.75499 10.2357 148.798 15.8572
1 4 2 -2.66145 5.17863 0 5.8225 90 117.2 4.59533 1.3053 -9.12124 10.2965 152.357 15.8572
1 4 3 2.66145 5.17863 0 5.8225 90 62.8 2.95046 -3.75708 -9.12124 10.2965 152.357 308.143
1 5 0 3.24711 6.31819 0 7.10375 90 62.8 3.59972 -4.58383 -8.38874 10.2147 145.209 308.143
1 5 1 -3.24711 6.31819 0 7.10375 90 117.2 5.60654 1.59254 -8.38874 10.2147 145.209 15.8572
1 5 2 -2.95428 5.74841 0 6.46312 90 117.2 5.10094 1.44892 -8.75499 10.2357 148.798 15.8572
1 5 3 2.95428 5.74841 0 6.46312 90 62.8 3.27509 -4.17046 -8.75499 10.2357 148.798 308.143
1 6 0 3.53994 6.88797 0 7.74438 90 62.8 3.92434 -4.99721 -8.0225 10.2339 141.62 308.143
1 6 1 -3.53994 6.88797 0 7.74438 90 117.2 6.11215 1.73615 -8.0225 10.2339 141.62 15.8572
1 6 2 -3.24711 6.31819 0 7.10375 90 117.2 5.60654 1.59254 -8.38874 10.2147 145.209 15.8572
1 6 3 3.24711 6.31819 0 7.10375 90 62.8 3.59972 -4.58383 -8.38874 10.2147 145.209 308.143
1 7 0 3.83277 7.45776 0 8.385 90 62.8 4.24897 -5.41058 -7.65625 10.293 138.059 308.143
1 7 1 -3.83277 7.45776 0 8.385 90 117.2 6.61775 1.87977 -7.65625 10.293 138.059 15.8572
1 7 2 -3.53994 6.88797 0 7.74438 90 117.2 6.11215 1.73615 -8.0225 10.2339 141.62 15.8572
1 7 3 3.53994 6.88797 0 7.74438 90 62.8 3.92434 -4.99721 -8.0225 10.2339 141.62 308.143
1 8 0 4.12559 8.02754 0 9.02562 90 62.8 4.5736 -5.82396 -7.29 10.3914 134.551 308.143
1 8 1 -4.12559 8.02754 0 9.02562 90 117.2 7.12335 2.02339 -7.29 10.3914 134.551 15.8572
1 8 2 -3.83277 7.45776 0 8.385 90 117.2 6.61775 1.87977 -7.65625 10.293 138.059 15.8572
1 8 3 3.83277 7.45776 0 8.385 90 62.8 4.24897 -5.41058 -7.65625 10.293 138.059 308.143
1 9 0 4.41842 8.59732 0 9.66625 90 62.8 4.89822 -6.23733 -6.92375 10.5278 131.122 308.143
1 9 1 -4.41842 8.59732 0 9.66625 90 117.2 7.62896 2.16701 -6.92375 10.5278 131.122 15.8572
1 9 2 -4.12559 8.02754 0 9.02562 90 117.2 7.12335 2.02339 -7.29 10.3914 134.551 15.8572
1 9 3 4.12559 8.02754 0 9.02562 90 62.8 4.5736 -5.82396 -7.29 10.3914 134.551 308.143
1 10 0 4.71125 9.1671 0 10.3069 90 62.8 5.22285 -6.65071 -6.5575 10.701 127.792 308.143
1 10 1 -4.71125 9.1671 0 10.3069 90 117.2 8.13456 2.31062 -6.5575 10.701 127.792 15.8572
1 10 2 -4.41842 8.59732 0 9.66625 90 117.2 7.62896 2.16701 -6.92375 10.5278 131.122 15.8572
1 10 3 4.41842 8.59732 0 9.66625 90 62.8 4.89822 -6.23733 -6.92375 10.5278 131.122 308.143
1 11 0 5.00408 9.73689 0 10.9475 90 62.8 5.54748 -7.06409 -6.19125 10.9091 124.578 308.143
1 11 1 -5.00408 9.73689 0 10.9475 90 117.2 8.64017 2.45424 -6.19125 10.9091 124.578 15.8572
1 11 2 -4.71125 9.1671 0 10.3069 90 117.2 8.13456 2.31062 -6.5575 10.701 127.792 15.8572
1 11 3 4.71125 9.1671 0 10.3069 90 62.8 5.22285 -6.65071 -6.5575 10.701 127.792 308.143
1 12 0 5.29691 10.3067 0 11.5881 90 62.8 5.8721 -7.47746 -5.825 11.1501 121.495 308.143
1 12 1 -5.29691 10.3067 0 11.5881 90 117.2 9.14577 2.59786 -5.825 11.1501 121.495 15.8572
1 12 2 -5.00408 9.73689 0 10.9475 90 117.2 8.64017 2.45424 -6.19125 10.9091 124.578 15.8572
1 12 3 5.00408 9.73689 0 10.9475 90 62.8 5.54748 -7.06409 -6.19125 10.9091 124.578 308.143
1 13 0 5.58974 10.8765 0 12.2287 90 62.8 6.19673 -7.89084 -5.45875 11.422 118.549 308.143
1 13 1 -5.58974 10.8765 0 12.2287 90 117.2 9.65138 2.74147 -5.45875 11.422 118.549 15.8572
1 13 2 -5.29691 10.3067 0 11.5881 90 117.2 9.14577 2.59786 -5.825 11.1501 121.495 15.8572
1 13 3 5.29691 10.3067 0 11.5881 90 62.8 5.8721 -7.47746 -5.825 11.1501 121.495 308.143
1 14 0 5.88256 11.4462 0 12.8694 90 62.8 6.52136 -8.30421 -5.0925 11.7227 115.748 308.143
1 14 1 -5.88256 11.4462 0 12.8694 90 117.2 10.157 2.88509 -5.0925 11.7227 115.748 15.8572
1 14 2 -5.58974 10.8765 0 12.2287 90 117.2 9.65138 2.74147 -5.45875 11.422 118.549 15.8572
1 14 3 5.58974 10.8765 0 12.2287 90 62.8 6.19673 -7.89084 -5.45875 11.422 118.549 308.143
1 15 0 6.17539 12.016 0 13.51 90 62.8 6.84599 -8.71759 -4.72625 12.05 113.093 308.143
1 15 1 -6.17539 12.016 0 13.51 90 117.2 10.6626 3.02871 -4.72625 12.05 113.093 15.8572
1 15 2 -5.88256 11.4462 0 12.8694 90 117.2 10.157 2.88509 -5.0925 11.7227 115.748 15.8572
1 15 3 5.88256 11.4462 0 12.8694 90 62.8 6.52136 -8.30421 -5.0925 11.7227 115.748 308.143
2 0 0 1.78297 3.46928 0 3.90062 90 62.8 -1.78297 -2.65762 -10.22 10.7093 162.613 236.143
2 0 1 -1.78297 3.46928 0 3.90062 90 117.2 1.78297 -2.65762 -10.22 10.7093 162.613 303.857
2 0 2 -1.49014 2.8995 0 3.26 90 117.2 1.49014 -2.22114 -10.5862 10.9189 165.821 303.857
2 0 3 1.49014 2.8995 0 3.26 90 62.8 -1.49014 -2.22114 -10.5862 10.9189 165.821 236.143
2 1 0 2.0758 4.03906 0 4.54125 90 62.8 -2.0758 -3.0941 -9.85374 10.5346 159.287 236.143
2 1 1 -2.0758 4.03906 0 4.54125 90 117.2 2.0758 -3.0941 -9.85374 10.5346 159.287 303.857
2 1 2 -1.78297 3.46928 0 3.90062 90 117.2 1.78297 -2.65762 -10.22 10.7093 162.613 303.857
2 1 3 1.78297 3.46928 0 3.90062 90 62.8 -1.78297 -2.65762 -10.22 10.7093 162.613 236.143
2 2 0 2.36862 4.60884 0 5.18187 90 62.8 -2.36862 -3.53058 -9.48749 10.3965 155.862 236.143
2 2 1 -2.36862 4.60884 0 5.18187 90 117.2 2.36862 -3.53058 -9.48749 10.3965 155.862 303.857
2 2 2 -2.0758 4.03906 0 4.54125 90 117.2 2.0758 -3.0941 -9.85374 10.5346 159.287 303.857
2 2 3 2.0758 4.03906 0 4.54125 90 62.8 -2.0758 -3.0941 -9.85374 10.5346 159.287 236.143
2 3 0 2.66145 5.17863 0 5.8225 90 62.8 -2.66145 -3.96706 -9.12124 10.2965 152.357 236.143
2 3 1 -2.66145 5.17863 0 5.8225 90 117.2 2.66145 -3.96706 -9.12124 10.2965 152.357 303.857
2 3 2 -2.36862 4.60884 0 5.18187 90 117.2 2.36862 -3.53058 -9.48749 10.3965 155.862 303.857
2 3 3 2.36862 4.60884 0 5.18187 90 62.8 -2.36862 -3.53058 -9.48749 10.3965 155.862 236.143
2 4 0 2.95428 5.74841 0 6.46312 90 62.8 -2.95428 -4.40354 -8.75499 10.2357 148.798 236.143
2 4 1 -2.95428 5.74841 0 6.46312 90 117.2 2.95428 -4.40354 -8.75499 10.2357 148.798 303.857
2 4 2 -2.66145 5.17863 0 5.8225 90 117.2 2.66145 -3.96706 -9.12124 10.2965 152.357 303.857
2 4 3 2.66145 5.17863 0 5.8225 90 62.8 -2.66145 -3.96706 -9.12124 10.2965 152.357 236.143
2 5 0 3.24711 6.31819 0 7.10375 90 62.8 -3.24711 -4.84002 -8.38874 10.2147 145.209 236.143
2 5 1 -3.24711 6.31819 0 7.10375 90 117.2 3.24711 -4.84002 -8.38874 10.2147 145.209 303.857
2 5 2 -2.95428 5.74841 0 6.46312 90 117.2 2.95428 -4.40354 -8.75499 10.2357 148.798 303.857
2 5 3 2.95428 5.74841 0 6.46312 90 62.8 -2.95428 -4.40354 -8.75499 10.2357 148.798 236.143
2 6 0 3.53994 6.88797 0 7.74438 90 62.8 -3.53994 -5.27649 -8.0225 10.2339 141.62 236.143
2 6 1 -3.53994 6.88797 0 7.74438 90 117.2 3.53994 -5.27649 -8.0225 10.2339 141.62 303.857
2 6 2 -3.24711 6.31819 0 7.10375 90 117.2 3.24711 -4.84002 -8.38874 10.2147 145.209 303.857
2 6 3 3.24711 6.31819 0 7.10375 90 62.8 -3.24711 -4.84002 -8.38874 10.2147 145.209 236.143
2 7 0 3.83277 7.45776 0 8.385 90 62.8 -3.83277 -5.71297 -7.65625 10.293 138.059 236.143
2 7 1 -3.83277 7.45776 0 8.385 90 117.2 3.83277 -5.71297 -7.65625 10.293 138.059 303.857
2 7 2 -3.53994 6.88797 0 7.74438 90 117.2 3.53994 -5.27649 -8.0225 10.2339 141.62 303.857
2 7 3 3.53994 6.88797 0 7.74438 90 62.8 -3.53994 -5.27649 -8.0225 10.2339 141.62 236.143
2 8 0 4.12559 8.02754 0 9.02562 90 62.8 -4.12559 -6.14945 -7.29 10.3914 134.551 236.143
2 8 1 -4.12559 8.02754 0 9.02562 90 117.2 4.12559 -6.14945 -7.29 10.3914 134.551 303.857
2 8 2 -3.83277 7.45776 0 8.385 90 117.2 3.83277 -5.71297 -7.65625 10.293 138.059 303.857
2 8 3 3.83277 7.45776 0 8.385 90 62.8 -3.83277 -5.71297 -7.65625 10.293 138.059 236.143
2 9 0 4.41842 8.59732 0 9.66625 90 62.8 -4.41842 -6.58593 -6.92375 10.5278 131.122 236.143
2 9 1 -4.41842 8.59732 0 9.66625 90 117.2 4.41842 -6.58593 -6.92375 10.5278 131.122 303.857
2 9 2 -4.12559 8.02754 0 9.02562 90 117.2 4.12559 -6.14945 -7.29 10.3914 134.551 303.857
2 9 3 4.12559 8.02754 0 9.02562 90 62.8 -4.12559 -6.14945 -7.29 10.3914 134.551 236.143
2 10 0 4.71125 9.1671 0 10.3069 90 62.8 -4.71125 -7.02241 -6.5575 10.701 127.792 236.143
2 10 1 -4.71125 9.1671 0 10.3069 90 117.2 4.71125 -7.02241 -6.5575 10.701 127.792 303.857
2 10 2 -4.41842 8.59732 0 9.66625 90 117.2 4.41842 -6.58593 -6.92375 10.5278 131.122 303.857
2 10 3 4.41842 8.59732 0 9.66625 90 62.8 -4.41842 -6.58593 -6.92375 10.5278 131.122 236.143
2 11 0 5.00408 9.73689 0 10.9475 90 62.8 -5.00408 -7.45889 -6.19125 10.9091 124.578 236.143
2 11 1 -5.00408 9.73689 0 10.9475 90 117.2 5.00408 -7.45889 -6.19125 10.9091 124.578 303.857
2 11 2 -4.71125 9.1671 0 10.3069 90 117.2 4.71125 -7.02241 -6.5575 10.701 127.792 303.857
2 11 3 4.71125 9.1671 0 10.3069 90 62.8 -4.71125 -7.02241 -6.5575 10.701 127.792 236.143
2 12 0 5.29691 10.3067 0 11.5881 90 62.8 -5.29691 -7.89537 -5.825 11.1501 121.495 236.143
2 12 1 -5.29691 10.3067 0 11.5881 90 117.2 5.29691 -7.89537 -5.825 11.1501 121.495 303.857
2 12 2 -5.00408 9.73689 0 10.9475 90 117.2 5.00408 -7.45889 -6.19125 10.9091 124.578 303.857
2 12 3 5.00408 9.73689 0 10.9475 90 62.8 -5.00408 -7.45889 -6.19125 10.9091 124.578 236.143
2 13 0 5.58974 10.8765 0 12.2287 90 62.8 -5.58974 -8.33184 -5.45875 11.422 118.549 236.143
2 13 1 -5.58974 10.8765 0 12.2287 90 117.2 5.58974 -8.33184 -5.45875 11.422 118.549 303.857
2 13 2 -5.29691 10.3067 0 11.5881 90 117.2 5.29691 -7.89537 -5.825 11.1501 121.495 303.857
2 13 3 5.29691 10.3067 0 11.5881 90 62.8 -5.29691 -7.89537 -5.825 11.1501 121.495 236.143
2 14 0 5.88256 11.4462 0 12.8694 90 62.8 -5.88256 -8.76832 -5.0925 11.7227 115.748 236.143
2 14 1 -5.88256 11.4462 0 12.8694 90 117.2 5.88256 -8.76832 -5.0925 11.7227 115.748 303.857
2 14 2 -5.58974 10.8765 0 12.2287 90 117.2 5.58974 -8.33184 -5.45875 11.422 118.549 303.857
2 14 3 5.58974 10.8765 0 12.2287 90 62.8 -5.58974 -8.33184 -5.45875 11.422 118.549 236.143
2 15 0 6.17539 12.016 0 13.51 90 62.8 -6.17539 -9.2048 -4.72625 12.05 113.093 236.143
2 15 1 -6.17539 12.016 0 13.51 90 117.2 6.17539 -9.2048 -4.72625 12.05 113.093 303.857
2 15 2 -5.88256 11.4462 0 12.8694 90 117.2 5.88256 -8.76832 -5.0925 11.7227 115.748 303.857
2 15 3 5.88256 11.4462 0 12.8694 90 62.8 -5.88256 -8.76832 -5.0925 11.7227 115.748 236.143
3 0 0 1.78297 3.46928 0 3.90062 90 62.8 -3.07852 0.874452 -10.22 10.7093 162.613 164.143
3 0 1 -1.78297 3.46928 0 3.90062 90 117.2 -1.97658 -2.51695 -10.22 10.7093 162.613 231.857
3 0 2 -1.49014 2.8995 0 3.26 90 117.2 -1.65195 -2.10358 -10.5862 10.9189 165.821 231.857
3 0 3 1.49014 2.8995 0 3.26 90 62.8 -2.57291 0.730835 -10.5862 10.9189 165.821 164.143
3 1 0 2.0758 4.03906 0 4.54125 90 62.8 -3.58412 1.01807 -9.85374 10.5346 159.287 164.143
3 1 1 -2.0758 4.03906 0 4.54125 90 117.2 -2.30121 -2.93033 -9.85374 10.5346 159.287 231.857
3 1 2 -1.78297 3.46928 0 3.90062 90 117.2 -1.97658 -2.51695 -10.22 10.7093 162.613 231.857
3 1 3 1.78297 3.46928 0 3.90062 90 62.8 -3.07852 0.874452 -10.22 10.7093 162.613 164.143
3 2 0 2.36862 4.60884 0 5.18187 90 62.8 -4.08973 1.16169 -9.48749 10.3965 155.862 164.143
3 2 1 -2.36862 4.60884 0 5.18187 90 117.2 -2.62584 -3.3437 -9.48749 10.3965 155.862 231.857
3 2 2 -2.0758 4.03906 0 4.54125 90 117.2 -2.30121 -2.93033 -9.85374 10.5346 159.287 231.857
3 2 3 2.0758 4.03906 0 4.54125 90 62.8 -3.58412 1.01807 -9.85374 10.5346 159.287 164.143
3 3 0 2.66145 5.17863 0 5.8225 90 62.8 -4.59533 1.3053 -9.12124 10.2965 152.357 164.143
3 3 1 -2.66145 5.17863 0 5.8225 90 117.2 -2.95046 -3.75708 -9.12124 10.2965 152.357 231.857
3 3 2 -2.36862 4.60884 0 5.18187 90 117.2 -2.62584 -3.3437 -9.48749 10.3965 155.862 231.857
3 3 3 2.36862 4.60884 0 5.18187 90 62.8 -4.08973 1.16169 -9.48749 10.3965 155.862 164.143
3 4 0 2.95428 5.74841 0 6.46312 90 62.8 -5.10094 1.44892 -8.75499 10.2357 148.798 164.143
3 4 1 -2.95428 5.74841 0 6.46312 90 117.2 -3.27509 -4.17046 -8.75499 10.2357 148.798 231.857
3 4 2 -2.66145 5.17863 0 5.8225 90 117.2 -2.95046 -3.75708 -9.12124 10.2965 152.357 231.857
3 4 3 2.66145 5.17863 0 5.8225 90 62.8 -4.59533 1.3053 -9.12124 10.2965 152.357 164.143
3 5 0 3.24711 6.31819 0 7.10375 90 62.8 -5.60654 1.59254 -8.38874 10.2147 145.209 164.143
3 5 1 -3.24711 6.31819 0 7.10375 90 117.2 -3.59972 -4.58383 -8.38874 10.2147 145.209 231.857
3 5 2 -2.95428 5.74841 0 6.46312 90 117.2 -3.27509 -4.17046 -8.75499 10.2357 148.798 231.857
3 5 3 2.95428 5.74841 0 6.46312 90 62.8 -5.10094 1.44892 -8.75499 10.2357 148.798 164.143
3 6 0 3.53994 6.88797 0 7.74438 90 62.8 -6.11215 1.73615 -8.0225 10.2339 141.62 164.143
3 6 1 -3.53994 6.88797 0 7.74438 90 117.2 -3.92434 -4.99721 -8.0225 10.2339 141.62 231.857
3 6 2 -3.24711 6.31819 0 7.10375 90 117.2 -3.59972 -4.58383 -8.38874 10.2147 145.209 231.857
3 6 3 3.24711 6.31819 0 7.10375 90 62.8 -5.60654 1.59254 -8.38874 10.2147 145.209 164.143
3 7 0 3.83277 7.45776 0 8.385 90 62.8 -6.61775 1.87977 -7.65625 10.293 138.059 164.143
3 7 1 -3.83277 7.45776 0 8.385 90 117.2 -4.24897 -5.41058 -7.65625 10.293 138.059 231.857
3 7 2 -3.53994 6.88797 0 7.74438 90 117.2 -3.92434 -4.99721 -8.0225 10.2339 141.62 231.857
3 7 3 3.53994 6.88797 0 7.74438 90 62.8 -6.11215 1.73615 -8.0225 10.2339 141.62 164.143
3 8 0 4.12559 8.02754 0 9.02562 90 62.8 -7.12335 2.02339 -7.29 10.3914 134.551 164.143
3 8 1 -4.12559 8.02754 0 9.02562 90 117.2 -4.5736 -5.82396 -7.29 10.3914 134.551 231.857
3 8 2 -3.83277 7.45776 0 8.385 90 117.2 -4.24897 -5.41058 -7.65625 10.293 138.059 231.857
3 8 3 3.83277 7.45776 0 8.385 90 62.8 -6.61775 1.87977 -7.65625 10.293 138.059 164.143
3 9 0 4.41842 8.59732 0 9.66625 90 62.8 -7.62896 2.16701 -6.92375 10.5278 131.122 164.143
3 9 1 -4.41842 8.59732 0 9.66625 90 117.2 -4.89822 -6.23733 -6.92375 10.5278 131.122 231.857
3 9 2 -4.12559 8.02754 0 9.02562 90 117.2 -4.5736 -5.82396 -7.29 10.3914 134.551 231.857
3 9 3 4.12559 8.02754 0 9.02562 90 62.8 -7.12335 2.02339 -7.29 10.3914 134.551 164.143
3 10 0 4.71125 9.1671 0 10.3069 90 62.8 -8.13456 2.31062 -6.5575 10.701 127.792 164.143
3 10 1 -4.71125 9.1671 0 10.3069 90 117.2 -5.22285 -6.65071 -6.5575 10.701 127.792 231.857
3 10 2 -4.41842 8.59732 0 9.66625 90 117.2 -4.89822 -6.23733 -6.92375 10.5278 131.122 231.857
3 10 3 4.41842 8.59732 0 9.66625 90 62.8 -7.62896 2.16701 -6.92375 10.5278 131.122 164.143
3 11 0 5.00408 9.73689 0 10.9475 90 62.8 -8.64017 2.45424 -6.19125 10.9091 124.578 164.143
3 11 1 -5.00408 9.73689 0 10.9475 90 117.2 -5.54748 -7.06409 -6.19125 10.9091 124.578 231.857
3 11 2 -4.71125 9.1671 0 10.3069 90 117.2 -5.22285 -6.65071 -6.5575 10.701 127.792 231.857
3 11 3 4.71125 9.1671 0 10.3069 90 62.8 -8.13456 2.31062 -6.5575 10.701 127.792 164.143
3 12 0 5.29691 10.3067 0 11.5881 90 62.8 -9.14577 2.59786 -5.825 11.1501 121.495 164.143
3 12 1 -5.29691 10.3067 0 11.5881 90 117.2 -5.8721 -7.47746 -5.825 11.1501 121.495 231.857
3 12 2 -5.00408 9.73689 0 10.9475 90 117.2 -5.54748 -7.06409 -6.19125 10.9091 124.578 231.857
3 12 3 5.00408 9.73689 0 10.9475 90 62.8 -8.64017 2.45424 -6.19125 10.9091 124.578 164.143
3 13 0 5.58974 10.8765 0 12.2287 90 62.8 -9.65138 2.74147 -5.45875 11.422 118.549 164.143
3 13 1 -5.58974 10.8765 0 12.2287 90 117.2 -6.19673 -7.89084 -5.45875 11.422 118.549 231.857
3 13 2 -5.29691 10.3067 0 11.5881 90 117.2 -5.8721 -7.47746 -5.825 11.1501 121.495 231.857
3 13 3 5.29691 10.3067 0 11.5881 90 62.8 -9.14577 2.59786 -5.825 11.1501 121.495 164.143
3 14 0 5.88256 11.4462 0 12.8694 90 62.8 -10.157 2.88509 -5.0925 11.7227 115.748 164.143
3 14 1 -5.88256 11.4462 0 12.8694 90 117.2 -6.52136 -8.30421 -5.0925 11.7227 115.748 231.857
3 14 2 -5.58974 10.8765 0 12.2287 90 117.2 -6.19673 -7.89084 -5.45875 11.422 118.549 231.857
3 14 3 5.58974 10.8765 0 12.2287 90 62.8 -9.65138 2.74147 -5.45875 11.422 118.549 164.143
3 15 0 6.17539 12.016 0 13.51 90 62.8 -10.6626 3.02871 -4.72625 12.05 113.093 164.143
3 15 1 -6.17539 12.016 0 13.51 90 117.2 -6.84599 -8.71759 -4.72625 12.05 113.093 231.857
3 15 2 -5.88256 11.4462 0 12.8694 90 117.2 -6.52136 -8.30421 -5.0925 11.7227 115.748 231.857
3 15 3 5.88256 11.4462 0 12.8694 90 62.8 -10.157 2.88509 -5.0925 11.7227 115.748 164.143
4 0 0 1.78297 3.46928 0 3.90062 90 62.8 -0.11966 3.19806 -10.22 10.7093 162.613 92.1428
4 0 1 -1.78297 3.46928 0 3.90062 90 117.2 -3.00456 1.10206 -10.22 10.7093 162.613 159.857
4 0 2 -1.49014 2.8995 0 3.26 90 117.2 -2.5111 0.921061 -10.5862 10.9189 165.821 159.857
4 0 3 1.49014 2.8995 0 3.26 90 62.8 -0.100008 2.67282 -10.5862 10.9189 165.821 92.1428
4 1 0 2.0758 4.03906 0 4.54125 90 62.8 -0.139313 3.7233 -9.85374 10.5346 159.287 92.1428
4 1 1 -2.0758 4.03906 0 4.54125 90 117.2 -3.49802 1.28306 -9.85374 10.5346 159.287 159.857
4 1 2 -1.78297 3.46928 0 3.90062 90 117.2 -3.00456 1.10206 -10.22 10.7093 162.613 159.857
4 1 3 1.78297 3.46928 0 3.90062 90 62.8 -0.11966 3.19806 -10.22 10.7093 162.613 92.1428
4 2 0 2.36862 4.60884 0 5.18187 90 62.8 -0.158965 4.24854 -9.48749 10.3965 155.862 92.1428
4 2 1 -2.36862 4.60884 0 5.18187 90 117.2 -3.99148 1.46406 -9.48749 10.3965 155.862 159.857
4 2 2 -2.0758 4.03906 0 4.54125 90 117.2 -3.49802 1.28306 -9.85374 10.5346 159.287 159.857
4 2 3 2.0758 4.03906 0 4.54125 90 62.8 -0.139313 3.7233 -9.85374 10.5346 159.287 92.1428
4 3 0 2.66145 5.17863 0 5.8225 90 62.8 -0.178618 4.77378 -9.12124 10.2965 152.357 92.1428
4 3 1 -2.66145 5.17863 0 5.8225 90 117.2 -4.48494 1.64505 -9.12124 10.2965 152.357 159.857
4 3 2 -2.36862 4.60884 0 5.18187 90 117.2 -3.99148 1.46406 -9.48749 10.3965 155.862 159.857
4 3 3 2.36862 4.60884 0 5.18187 90 62.8 -0.158965 4.24854 -9.48749 10.3965 155.862 92.1428
4 4 0 2.95428 5.74841 0 6.46312 90 62.8 -0.19827 5.29902 -8.75499 10.2357 148.798 92.1428
4 4 1 -2.95428 5.74841 0 6.46312 90 117.2 -4.9784 1.82605 -8.75499 10.2357 148.798 159.857
4 4 2 -2.66145 5.17863 0 5.8225 90 117.2 -4.48494 1.64505 -9.12124 10.2965 152.357 159.857
4 4 3 2.66145 5.17863 0 5.8225 90 62.8 -0.178618 4.77378 -9.12124 10.2965 152.357 92.1428
4 5 0 3.24711 6.31819 0 7.10375 90 62.8 -0.217923 5.82426 -8.38874 10.2147 145.209 92.1428
4 5 1 -3.24711 6.31819 0 7.10375 90 117.2 -5.47186 2.00705 -8.38874 10.2147 145.209 159.857
4 5 2 -2.95428 5.74841 0 6.46312 90 117.2 -4.9784 1.82605 -8.75499 10.2357 148.798 159.857
4 5 3 2.95428 5.74841 0 6.46312 90 62.8 -0.19827 5.29902 -8.75499 10.2357 148.798 92.1428
4 6 0 3.53994 6.88797 0 7.74438 90 62.8 -0.237576 6.3495 -8.0225 10.2339 141.62 92.1428
4 6 1 -3.53994 6.88797 0 7.74438 90 117.2 -5.96532 2.18805 -8.0225 10.2339 141.62 159.857
4 6 2 -3.24711 6.31819 0 7.10375 90 117.2 -5.47186 2.00705 -8.38874 10.2147 145.209 159.857
4 6 3 3.24711 6.31819 0 7.10375 90 62.8 -0.217923 5.82426 -8.38874 10.2147 145.209 92.1428
4 7 0 3.83277 7.45776 0 8.385 90 62.8 -0.257228 6.87474 -7.65625 10.293 138.059 92.1428
4 7 1 -3.83277 7.45776 0 8.385 90 117.2 -6.45877 2.36905 -7.65625 10.293 138.059 159.857
4 7 2 -3.53994 6.88797 0 7.74438 90 117.2 -5.96532 2.18805 -8.0225 10.2339 141.62 159.857
4 7 3 3.53994 6.88797 0 7.74438 90 62.8 -0.237576 6.3495 -8.0225 10.2339 141.62 92.1428
4 8 0 4.12559 8.02754 0 9.02562 90 62.8 -0.276881 7.39997 -7.29 10.3914 134.551 92.1428
4 8 1 -4.12559 8.02754 0 9.02562 90 117.2 -6.95223 2.55005 -7.29 10.3914 134.551 159.857
4 8 2 -3.83277 7.45776 0 8.385 90 117.2 -6.45877 2.36905 -7.65625 10.293 138.059 159.857
4 8 3 3.83277 7.45776 0 8.385 90 62.8 -0.257228 6.87474 -7.65625 10.293 138.059 92.1428
4 9 0 4.41842 8.59732 0 9.66625 90 62.8 -0.296533 7.92521 -6.92375 10.5278 131.122 92.1428
4 9 1 -4.41842 8.59732 0 9.66625 90 117.2 -7.44569 2.73105 -6.92375 10.5278 131.122 159.857
4 9 2 -4.12559 8.02754 0 9.02562 90 117.2 -6.95223 2.55005 -7.29 10.3914 134.551 159.857
4 9 3 4.12559 8.02754 0 9.02562 90 62.8 -0.276881 7.39997 -7.29 10.3914 134.551 92.1428
4 10 0 4.71125 9.1671 0 10.3069 90 62.8 -0.316186 8.45045 -6.5575 10.701 127.792 92.1428
4 10 1 -4.71125 9.1671 0 10.3069 90 117.2 -7.93915 2.91204 -6.5575 10.701 127.792 159.857
4 10 2 -4.41842 8.59732 0 9.66625 90 117.2 -7.44569 2.73105 -6.92375 10.5278 131.122 159.857
4 10 3 4.41842 8.59732 0 9.66625 90 62.8 -0.296533 7.92521 -6.92375 10.5278 131.122 92.1428
4 11 0 5.00408 9.73689 0 10.9475 90 62.8 -0.335839 8.97569 -6.19125 10.9091 124.578 92.1428
4 11 1 -5.00408 9.73689 0 10.9475 90 117.2 -8.43261 3.09304 -6.19125 10.9091 124.578 159.857
4 11 2 -4.71125 9.1671 0 10.3069 90 117.2 -7.93915 2.91204 -6.5575 10.701 127.792 159.857
4 11 3 4.71125 9.1671 0 10.3069 90 62.8 -0.316186 8.45045 -6.5575 10.701 127.792 92.1428
4 12 0 5.29691 10.3067 0 11.5881 90 62.8 -0.355491 9.50093 -5.825 11.1501 121.495 92.1428
4 12 1 -5.29691 10.3067 0 11.5881 90 117.2 -8.92607 3.27404 -5.825 11.1501 121.495 159.857
4 12 2 -5.00408 9.73689 0 10.9475 90 117.2 -8.43261 3.09304 -6.19125 10.9091 124.578 159.857
4 12 3 5.00408 9.73689 0 10.9475 90 62.8 -0.335839 8.97569 -6.19125 10.9091 124.578 92.1428
4 13 0 5.58974 10.8765 0 12.2287 90 62.8 -0.375144 10.0262 -5.45875 11.422 118.549 92.1428
4 13 1 -5.58974 10.8765 0 12.2287 90 117.2 -9.41953 3.45504 -5.45875 11.422 118.549 159.857
4 13 2 -5.29691 10.3067 0 11.5881 90 117.2 -8.92607 3.27404 -5.825 11.1501 121.495 159.857
4 13 3 5.29691 10.3067 0 11.5881 90 62.8 -0.355491 9.50093 -5.825 11.1501 121.495 92.1428
4 14 0 5.88256 11.4462 0 12.8694 90 62.8 -0.394796 10.5514 -5.0925 11.7227 115.748 92.1428
4 14 1 -5.88256 11.4462 0 12.8694 90 117.2 -9.91299 3.63604 -5.0925 11.7227 115.748 159.857
4 14 2 -5.58974 10.8765 0 12.2287 90 117.2 -9.41953 3.45504 -5.45875 11.422 118.549 159.857
4 14 3 5.58974 10.8765 0 12.2287 90 62.8 -0.375144 10.0262 -5.45875 11.422 118.549 92.1428
4 15 0 6.17539 12.016 0 13.51 90 62.8 -0.414449 11.0766 -4.72625 12.05 113.093 92.1428
4 15 1 -6.17539 12.016 0 13.51 90 117.2 -10.4064 3.81704 -4.72625 12.05 113.093 159.857
4 15 2 -5.88256 11.4462 0 12.8694 90 117.2 -9.91299 3.63604 -5.0925 11.7227 115.748 159.857
4 15 3 5.88256 11.4462 0 12.8694 90 62.8 -0.394796 10.5514 -5.0925 11.7227 115.748 92.1428
WEDGES
det ch corner x_fl y_fl z_fl r_fl t_fl p_fl x_ti y_ti z_ti r_ti t_ti p_ti
0 0 0 -4.70921 12.6627 0 13.51 90 110.4 1.89179 10.6156 -4.31059 11.6125 111.79 79.8955
0 0 1 -6.17539 12.016 0 13.51 90 117.2 0.414449 11.0766 -4.72625 12.05 113.093 87.8572
0 0 2 -1.49014 2.8995 0 3.26 90 117.2 0.100008 2.67282 -10.5862 10.9189 165.821 87.8572
0 0 3 -1.13634 3.05554 0 3.26 90 110.4 0.456494 2.56158 -10.4859 10.8039 166.064 79.8955
0 1 0 -3.17677 13.1312 0 13.51 90 103.6 3.34252 10.0052 -4.00943 11.2851 110.811 71.5267
0 1 1 -4.70921 12.6627 0 13.51 90 110.4 1.89179 10.6156 -4.31059 11.6125 111.79 79.8955
0 1 2 -1.13634 3.05554 0 3.26 90 110.4 0.456494 2.56158 -10.4859 10.8039 166.064 79.8955
0 1 3 -0.766563 3.16859 0 3.26 90 103.6 0.806558 2.41429 -10.4133 10.7199 166.264 71.5267
0 2 0 -1.59964 13.415 0 13.51 90 96.8 4.74622 9.25407 -3.82703 11.082 110.202 62.8477
0 2 1 -3.17677 13.1312 0 13.51 90 103.6 3.34252 10.0052 -4.00943 11.2851 110.811 71.5267
0 2 2 -0.766563 3.16859 0 3.26 90 103.6 0.806558 2.41429 -10.4133 10.7199 166.264 71.5267
0 2 3 -0.385997 3.23707 0 3.26 90 96.8 1.14528 2.23303 -10.3693 10.6686 166.395 62.8477
0 3 0 0 13.51 0 13.51 90 90 6.08314 8.37273 -3.76594 11.0132 109.996 54
0 3 1 -1.59964 13.415 0 13.51 90 96.8 4.74622 9.25407 -3.82703 11.082 110.202 62.8477
0 3 2 -0.385997 3.23707 0 3.26 90 96.8 1.14528 2.23303 -10.3693 10.6686 166.395 62.8477
0 3 3 0 3.26 0 3.26 90 90 1.46788 2.02036 -10.3545 10.6514 166.44 54
0 4 0 1.59964 13.415 0 13.51 90 83.2 7.33448 7.37359 -3.82703 11.082 110.202 45.1523
0 4 1 0 13.51 0 13.51 90 90 6.08314 8.37273 -3.76594 11.0132 109.996 54
0 4 2 0 3.26 0 3.26 90 90 1.46788 2.02036 -10.3545 10.6514 166.44 54
0 4 3 0.385997 3.23707 0 3.26 90 83.2 1.76983 1.77927 -10.3693 10.6686 166.395 45.1523
0 5 0 3.17677 13.1312 0 13.51 90 76.4 8.48264 6.27071 -4.00943 11.2851 110.811 36.4733
0 5 1 1.59964 13.415 0 13.51 90 83.2 7.33448 7.37359 -3.82703 11.082 110.202 45.1523
0 5 2 0.385997 3.23707 0 3.26 90 83.2 1.76983 1.77927 -10.3693 10.6686 166.395 45.1523
0 5 3 0.766563 3.16859 0 3.26 90 76.4 2.04688 1.51314 -10.4133 10.7199 166.264 36.4733
0 6 0 4.70921 12.6627 0 13.51 90 69.6 9.51145 5.0796 -4.31059 11.6125 111.79 28.1045
0 6 1 3.17677 13.1312 0 13.51 90 76.4 8.48264 6.27071 -4.00943 11.2851 110.811 36.4733
0 6 2 0.766563 3.16859 0 3.26 90 76.4 2.04688 1.51314 -10.4133 10.7199 166.264 36.4733
0 6 3 1.13634 3.05554 0 3.26 90 69.6 2.29514 1.22572 -10.4859 10.8039 166.064 28.1045
0 7 0 6.17539 12.016 0 13.51 90 62.8 10.4064 3.81704 -4.72625 12.05 113.093 20.1428
0 7 1 4.70921 12.6627 0 13.51 90 69.6 9.51145 5.0796 -4.31059 11.6125 111.79 28.1045
0 7 2 1.13634 3.05554 0 3.26 90 69.6 2.29514 1.22572 -10.4859 10.8039 166.064 28.1045
0 7 3 1.49014 2.8995 0 3.26 90 62.8 2.5111 0.921061 -10.5862 10.9189 165.821 20.1428
1 0 0 -4.70921 12.6627 0 13.51 90 110.4 10.6806 1.4812 -4.31059 11.6125 111.79 7.89549
1 0 1 -6.17539 12.016 0 13.51 90 117.2 10.6626 3.02871 -4.72625 12.05 113.093 15.8572
1 0 2 -1.49014 2.8995 0 3.26 90 117.2 2.57291 0.730835 -10.5862 10.9189 165.821 15.8572
1 0 3 -1.13634 3.05554 0 3.26 90 110.4 2.57727 0.357419 -10.4859 10.8039 166.064 7.89549
1 1 0 -3.17677 13.1312 0 13.51 90 103.6 10.5484 -0.0871381 -4.00943 11.2851 110.811 359.527
1 1 1 -4.70921 12.6627 0 13.51 90 110.4 10.6806 1.4812 -4.31059 11.6125 111.79 7.89549
1 1 2 -1.13634 3.05554 0 3.26 90 110.4 2.57727 0.357419 -10.4859 10.8039 166.064 7.89549
1 1 3 -0.766563 3.16859 0 3.26 90 103.6 2.54536 -0.0210267 -10.4133 10.7199 166.264 359.527
1 2 0 -1.59964 13.415 0 13.51 90 96.8 10.2678 -1.65425 -3.82703 11.082 110.202 350.848
1 2 1 -3.17677 13.1312 0 13.51 90 103.6 10.5484 -0.0871381 -4.00943 11.2851 110.811 359.527
1 2 2 -0.766563 3.16859 0 3.26 90 103.6 2.54536 -0.0210267 -10.4133 10.7199 166.264 359.527
1 2 3 -0.385997 3.23707 0 3.26 90 96.8 2.47765 -0.399176 -10.3693 10.6686 166.395 350.848
1 3 0 0 13.51 0 13.51 90 90 9.84273 -3.1981 -3.76594 11.0132 109.996 342
1 3 1 -1.59964 13.415 0 13.51 90 96.8 10.2678 -1.65425 -3.82703 11.082 110.202 350.848
1 3 2 -0.385997 3.23707 0 3.26 90 96.8 2.47765 -0.399176 -10.3693 10.6686 166.395 350.848
1 3 3 0 3.26 0 3.26 90 90 2.37508 -0.77171 -10.3545 10.6514 166.44 342
1 4 0 1.59964 13.415 0 13.51 90 83.2 9.27918 -4.69695 -3.82703 11.082 110.202 333.152
1 4 1 0 13.51 0 13.51 90 90 9.84273 -3.1981 -3.76594 11.0132 109.996 342
1 4 2 0 3.26 0 3.26 90 90 2.37508 -0.77171 -10.3545 10.6514 166.44 342
1 4 3 0.385997 3.23707 0 3.26 90 83.2 2.23909 -1.13339 -10.3693 10.6686 166.395 333.152
1 5 0 3.17677 13.1312 0 13.51 90 76.4 8.58508 -6.12971 -4.00943 11.2851 110.811 324.473
1 5 1 1.59964 13.415 0 13.51 90 83.2 9.27918 -4.69695 -3.82703 11.082 110.202 333.152
1 5 2 0.385997 3.23707 0 3.26 90 83.2 2.23909 -1.13339 -10.3693 10.6686 166.395 333.152
1 5 3 0.766563 3.16859 0 3.26 90 76.4 2.0716 -1.47912 -10.4133 10.7199 166.264 324.473
1 6 0 4.70921 12.6627 0 13.51 90 69.6 7.77019 -7.47624 -4.31059 11.6125 111.79 316.105
1 6 1 3.17677 13.1312 0 13.51 90 76.4 8.58508 -6.12971 -4.00943 11.2851 110.811 324.473
1 6 2 0.766563 3.16859 0 3.26 90 76.4 2.0716 -1.47912 -10.4133 10.7199 166.264 324.473
1 6 3 1.13634 3.05554 0 3.26 90 69.6 1.87497 -1.80404 -10.4859 10.8039 166.064 316.105
1 7 0 6.17539 12.016 0 13.51 90 62.8 6.84599 -8.71759 -4.72625 12.05 113.093 308.143
1 7 1 4.70921 12.6627 0 13.51 90 69.6 7.77019 -7.47624 -4.31059 11.6125 111.79 316.105
1 7 2 1.13634 3.05554 0 3.26 90 69.6 1.87497 -1.80404 -10.4859 10.8039 166.064 316.105
1 7 3 1.49014 2.8995 0 3.26 90 62.8 1.65195 -2.10358 -10.5862 10.9189 165.821 308.143
2 0 0 -4.70921 12.6627 0 13.51 90 110.4 4.70921 -9.70018 -4.31059 11.6125 111.79 295.895
2 0 1 -6.17539 12.016 0 13.51 90 117.2 6.17539 -9.2048 -4.72625 12.05 113.093 303.857
2 0 2 -1.49014 2.8995 0 3.26 90 117.2 1.49014 -2.22114 -10.5862 10.9189 165.821 303.857
2 0 3 -1.13634 3.05554 0 3.26 90 110.4 1.13634 -2.34068 -10.4859 10.8039 166.064 295.895
2 1 0 -3.17677 13.1312 0 13.51 90 103.6 3.17677 -10.0591 -4.00943 11.2851 110.811 287.527
2 1 1 -4.70921 12.6627 0 13.51 90 110.4 4.70921 -9.70018 -4.31059 11.6125 111.79 295.895
2 1 2 -1.13634 3.05554 0 3.26 90 110.4 1.13634 -2.34068 -10.4859 10.8039 166.064 295.895
2 1 3 -0.766563 3.16859 0 3.26 90 103.6 0.766563 -2.42728 -10.4133 10.7199 166.264 287.527
2 2 0 -1.59964 13.415 0 13.51 90 96.8 1.59964 -10.2765 -3.82703 11.082 110.202 278.848
2 2 1 -3.17677 13.1312 0 13.51 90 103.6 3.17677 -10.0591 -4.00943 11.2851 110.811 287.527
2 2 2 -0.766563 3.16859 0 3.26 90 103.6 0.766563 -2.42728 -10.4133 10.7199 166.264 287.527
2 2 3 -0.385997 3.23707 0 3.26 90 96.8 0.385997 -2.47974 -10.3693 10.6686 166.395 278.848
2 3 0 0 13.51 0 13.51 90 90 1.26742e-15 -10.3493 -3.76594 11.0132 109.996 270
2 3 1 -1.59964 13.415 0 13.51 90 96.8 1.59964 -10.2765 -3.82703 11.082 110.202 278.848
2 3 2 -0.385997 3.23707 0 3.26 90 96.8 0.385997 -2.47974 -10.3693 10.6686 166.395 278.848
2 3 3 0 3.26 0 3.26 90 90 3.05832e-16 -2.4973 -10.3545 10.6514 166.44 270
2 4 0 1.59964 13.415 0 13.51 90 83.2 -1.59964 -10.2765 -3.82703 11.082 110.202 261.152
2 4 1 0 13.51 0 13.51 90 90 1.26742e-15 -10.3493 -3.76594 11.0132 109.996 270
2 4 2 0 3.26 0 3.26 90 90 3.05832e-16 -2.4973 -10.3545 10.6514 166.44 270
2 4 3 0.385997 3.23707 0 3.26 90 83.2 -0.385997 -2.47974 -10.3693 10.6686 166.395 261.152
2 5 0 3.17677 13.1312 0 13.51 90 76.4 -3.17677 -10.0591 -4.00943 11.2851 110.811 252.473
2 5 1 1.59964 13.415 0 13.51 90 83.2 -1.59964 -10.2765 -3.82703 11.082 110.202 261.152
2 5 2 0.385997 3.23707 0 3.26 90 83.2 -0.385997 -2.47974 -10.3693 10.6686 166.395 261.152
2 5 3 0.766563 3.16859 0 3.26 90 76.4 -0.766563 -2.42728 -10.4133 10.7199 166.264 252.473
2 6 0 4.70921 12.6627 0 13.51 90 69.6 -4.70921 -9.70018 -4.31059 11.6125 111.79 244.105
2 6 1 3.17677 13.1312 0 13.51 90 76.4 -3.17677 -10.0591 -4.00943 11.2851 110.811 252.473
2 6 2 0.766563 3.16859 0 3.26 90 76.4 -0.766563 -2.42728 -10.4133 10.7199 166.264 252.473
2 6 3 1.13634 3.05554 0 3.26 90 69.6 -1.13634 -2.34068 -10.4859 10.8039 166.064 244.105
2 7 0 6.17539 12.016 0 13.51 90 62.8 -6.17539 -9.2048 -4.72625 12.05 113.093 236.143
2 7 1 4.70921 12.6627 0 13.51 90 69.6 -4.70921 -9.70018 -4.31059 11.6125 111.79 244.105
2 7 2 1.13634 3.05554 0 3.26 90 69.6 -1.13634 -2.34068 -10.4859 10.8039 166.064 244.105
2 7 3 1.49014 2.8995 0 3.26 90 62.8 -1.49014 -2.22114 -10.5862 10.9189 165.821 236.143
3 0 0 -4.70921 12.6627 0 13.51 90 110.4 -7.77019 -7.47624 -4.31059 11.6125 111.79 223.895
3 0 1 -6.17539 12.016 0 13.51 90 117.2 -6.84599 -8.71759 -4.72625 12.05 113.093 231.857
3 0 2 -1.49014 2.8995 0 3.26 90 117.2 -1.65195 -2.10358 -10.5862 10.9189 165.821 231.857
3 0 3 -1.13634 3.05554 0 3.26 90 110.4 -1.87497 -1.80404 -10.4859 10.8039 166.064 223.895
3 1 0 -3.17677 13.1312 0 13.51 90 103.6 -8.58508 -6.12971 -4.00943 11.2851 110.811 215.527
3 1 1 -4.70921 12.6627 0 13.51 90 110.4 -7.77019 -7.47624 -4.31059 11.6125 111.79 223.895
3 1 2 -1.13634 3.05554 0 3.26 90 110.4 -1.87497 -1.80404 -10.4859 10.8039 166.064 223.895
3 1 3 -0.766563 3.16859 0 3.26 90 103.6 -2.0716 -1.47912 -10.4133 10.7199 166.264 215.527
3 2 0 -1.59964 13.415 0 13.51 90 96.8 -9.27918 -4.69695 -3.82703 11.082 110.202 206.848
3 2 1 -3.17677 13.1312 0 13.51 90 103.6 -8.58508 -6.12971 -4.00943 11.2851 110.811 215.527
3 2 2 -0.766563 3.16859 0 3.26 90 103.6 -2.0716 -1.47912 -10.4133 10.7199 166.264 215.527
3 2 3 -0.385997 3.23707 0 3.26 90 96.8 -2.23909 -1.13339 -10.3693 10.6686 166.395 206.848
3 3 0 0 13.51 0 13.51 90 90 -9.84273 -3.1981 -3.76594 11.0132 109.996 198
3 3 1 -1.59964 13.415 0 13.51 90 96.8 -9.27918 -4.69695 -3.82703 11.082 110.202 206.848
3 3 2 -0.385997 3.23707 0 3.26 90 96.8 -2.23909 -1.13339 -10.3693 10.6686 166.395 206.848
3 3 3 0 3.26 0 3.26 90 90 -2.37508 -0.77171 -10.3545 10.6514 166.44 198
3 4 0 1.59964 13.415 0 13.51 90 83.2 -10.2678 -1.65425 -3.82703 11.082 110.202 189.152
3 4 1 0 13.51 0 13.51 90 90 -9.84273 -3.1981 -3.76594 11.0132 109.996 198
3 4 2 0 3.26 0 3.26 90 90 -2.37508 -0.77171 -10.3545 10.6514 166.44 198
3 4 3 0.385997 3.23707 0 3.26 90 83.2 -2.47765 -0.399176 -10.3693 10.6686 166.395 189.152
3 5 0 3.17677 13.1312 0 13.51 90 76.4 -10.5484 -0.0871381 -4.00943 11.2851 110.811 180.473
3 5 1 1.59964 13.415 0 13.51 90 83.2 -10.2678 -1.65425 -3.82703 11.082 110.202 189.152
3 5 2 0.385997 3.23707 0 3.26 90 83.2 -2.47765 -0.399176 -10.3693 10.6686 166.395 189.152
3 5 3 0.766563 3.16859 0 3.26 90 76.4 -2.54536 -0.0210267 -10.4133 10.7199 166.264 180.473
3 6 0 4.70921 12.6627 0 13.51 90 69.6 -10.6806 1.4812 -4.31059 11.6125 111.79 172.105
3 6 1 3.17677 13.1312 0 13.51 90 76.4 -10.5484 -0.0871381 -4.00943 11.2851 110.811 180.473
3 6 2 0.766563 3.16859 0 3.26 90 76.4 -2.54536 -0.0210267 -10.4133 10.7199 166.264 180.473
3 6 3 1.13634 3.05554 0 3.26 90 69.6 -2.57727 0.357419 -10.4859 10.8039 166.064 172.105
3 7 0 6.17539 12.016 0 13.51 90 62.8 -10.6626 3.02871 -4.72625 12.05 113.093 164.143
3 7 1 4.70921 12.6627 0 13.51 90 69.6 -10.6806 1.4812 -4.31059 11.6125 111.79 172.105
3 7 2 1.13634 3.05554 0 3.26 90 69.6 -2.57727 0.357419 -10.4859 10.8039 166.064 172.105
3 7 3 1.49014 2.8995 0 3.26 90 62.8 -2.57291 0.730835 -10.5862 10.9189 165.821 164.143
4 0 0 -4.70921 12.6627 0 13.51 90 110.4 -9.51145 5.0796 -4.31059 11.6125 111.79 151.895
4 0 1 -6.17539 12.016 0 13.51 90 117.2 -10.4064 3.81704 -4.72625 12.05 113.093 159.857
4 0 2 -1.49014 2.8995 0 3.26 90 117.2 -2.5111 0.921061 -10.5862 10.9189 165.821 159.857
4 0 3 -1.13634 3.05554 0 3.26 90 110.4 -2.29514 1.22572 -10.4859 10.8039 166.064 151.895
4 1 0 -3.17677 13.1312 0 13.51 90 103.6 -8.48264 6.27071 -4.00943 11.2851 110.811 143.527
4 1 1 -4.70921 12.6627 0 13.51 90 110.4 -9.51145 5.0796 -4.31059 11.6125 111.79 151.895
4 1 2 -1.13634 3.05554 0 3.26 90 110.4 -2.29514 1.22572 -10.4859 10.8039 166.064 151.895
4 1 3 -0.766563 3.16859 0 3.26 90 103.6 -2.04688 1.51314 -10.4133 10.7199 166.264 143.527
4 2 0 -1.59964 13.415 0 13.51 90 96.8 -7.33448 7.37359 -3.82703 11.082 110.202 134.848
4 2 1 -3.17677 13.1312 0 13.51 90 103.6 -8.48264 6.27071 -4.00943 11.2851 110.811 143.527
4 2 2 -0.766563 3.16859 0 3.26 90 103.6 -2.04688 1.51314 -10.4133 10.7199 166.264 143.527
4 2 3 -0.385997 3.23707 0 3.26 90 96.8 -1.76983 1.77927 -10.3693 10.6686 166.395 134.848
4 3 0 0 13.51 0 13.51 90 90 -6.08314 8.37273 -3.76594 11.0132 109.996 126
4 3 1 -1.59964 13.415 0 13.51 90 96.8 -7.33448 7.37359 -3.82703 11.082 110.202 134.848
4 3 2 -0.385997 3.23707 0 3.26 90 96.8 -1.76983 1.77927 -10.3693 10.6686 166.395 134.848
4 3 3 0 3.26 0 3.26 90 90 -1.46788 2.02036 -10.3545 10.6514 166.44 126
4 4 0 1.59964 13.415 0 13.51 90 83.2 -4.74622 9.25407 -3.82703 11.082 110.202 117.152
4 4 1 0 13.51 0 13.51 90 90 -6.08314 8.37273 -3.76594 11.0132 109.996 126
4 4 2 0 3.26 0 3.26 90 90 -1.46788 2.02036 -10.3545 10.6514 166.44 126
4 4 3 0.385997 3.23707 0 3.26 90 83.2 -1.14528 2.23303 -10.3693 10.6686 166.395 117.152
4 5 0 3.17677 13.1312 0 13.51 90 76.4 -3.34252 10.0052 -4.00943 11.2851 110.811 108.473
4 5 1 1.59964 13.415 0 13.51 90 83.2 -4.74622 9.25407 -3.82703 11.082 110.202 117.152
4 5 2 0.385997 3.23707 0 3.26 90 83.2 -1.14528 2.23303 -10.3693 10.6686 166.395 117.152
4 5 3 0.766563 3.16859 0 3.26 90 76.4 -0.806558 2.41429 -10.4133 10.7199 166.264 108.473
4 6 0 4.70921 12.6627 0 13.51 90 69.6 -1.89179 10.6156 -4.31059 11.6125 111.79 100.105
4 6 1 3.17677 13.1312 0 13.51 90 76.4 -3.34252 10.0052 -4.00943 11.2851 110.811 108.473
4 6 2 0.766563 3.16859 0 3.26 90 76.4 -0.806558 2.41429 -10.4133 10.7199 166.264 108.473
4 6 3 1.13634 3.05554 0 3.26 90 69.6 -0.456494 2.56158 -10.4859 10.8039 166.064 100.105
4 7 0 6.17539 12.016 0 13.51 90 62.8 -0.414449 11.0766 -4.72625 12.05 113.093 92.1428
4 7 1 4.70921 12.6627 0 13.51 90 69.6 -1.89179 10.6156 -4.31059 11.6125 111.79 100.105
4 7 2 1.13634 3.05554 0 3.26 90 69.6 -0.456494 2.56158 -10.4859 10.8039 166.064 100.105
4 7 3 1.49014 2.8995 0 3.26 90 62.8 -0.100008 2.67282 -10.5862 10.9189 165.821 92.1428

1
include/EnergyLoss.h
View File

@ -17,7 +17,6 @@ Written by G.W. McCann Aug. 2020
#include <vector> #include <vector>
#include <fstream> #include <fstream>
#include <cmath> #include <cmath>
#include "Eloss_Tables.h"
#include "MassLookup.h" #include "MassLookup.h"
class EnergyLoss { class EnergyLoss {

50
include/G3Vec.h Normal file
View File

@ -0,0 +1,50 @@
#ifndef G3VEC_H
#define G3VEC_H
#include <cmath>
class G3Vec {
public:
G3Vec();
G3Vec(double x, double y, double z);
~G3Vec();
void SetVectorCartesian(double x, double y, double z);
void SetVectorSpherical(double r, double theta, double phi);
inline double GetX() const { return m_data[0]; };
inline double GetY() const { return m_data[1]; };
inline double GetZ() const { return m_data[2]; };
inline double GetRho() const { return std::sqrt(std::pow(m_data[0], 2.0) + std::pow(m_data[1], 2.0)); };
inline double GetR() const { return std::sqrt(std::pow(m_data[0], 2.0) + std::pow(m_data[1], 2.0) + std::pow(m_data[2], 2.0)); }
inline double GetTheta() const { return Atan2(GetRho(), GetZ()); };
inline double GetPhi() const {
double phi = Atan2(GetY(), GetX());
if(phi < 0) phi += M_PI*2.0;
return phi;
};
inline const double operator[](int index) const { return index>2 || index<0 ? 0.0 : m_data[index]; };
inline G3Vec& operator=(const G3Vec& rhs) { SetVectorCartesian(rhs.GetX(), rhs.GetY(), rhs.GetZ()); return *this; };
inline G3Vec operator+(const G3Vec& rhs) const { return G3Vec(this->GetX()+rhs.GetX(), this->GetY()+rhs.GetY(), this->GetZ()+rhs.GetZ()); };
inline G3Vec operator-(const G3Vec& rhs) const { return G3Vec(this->GetX()-rhs.GetX(), this->GetY()-rhs.GetY(), this->GetZ()-rhs.GetZ()); };
double Dot(const G3Vec& rhs) const;
G3Vec Cross(const G3Vec& rhs) const;
private:
inline double Atan2(double y, double x) const {
if(x != 0.0) return std::atan2(y, x);
else if(y > 0.0) return M_PI/2.0;
else if(y < 0.0) return -M_PI/2.0;
else return 0.0;
}
double m_data[3];
};
#endif

17
include/G4Vec.h
View File

@ -15,13 +15,12 @@ public:
inline double GetPy() const {return m_data[1];}; inline double GetPy() const {return m_data[1];};
inline double GetPz() const {return m_data[2];}; inline double GetPz() const {return m_data[2];};
inline double GetP() const {return std::sqrt(m_data[0]*m_data[0] + m_data[1]*m_data[1] + m_data[2]*m_data[2]);}; inline double GetP() const {return std::sqrt(m_data[0]*m_data[0] + m_data[1]*m_data[1] + m_data[2]*m_data[2]);};
inline double GetTheta() const {return GetP() == 0.0 ? 0.0 : acos(GetPz()/GetP());}; inline double GetPxy() const {return std::sqrt(m_data[0]*m_data[0] + m_data[1]*m_data[1]); };
inline double GetTheta() const {return GetPxy() == 0.0 && GetPz() == 0.0 ? 0.0 : Atan2(GetPxy(), GetPz());};
inline double GetPhi() const { inline double GetPhi() const {
if(GetPx() == 0) return M_PI/2.0; double phi = Atan2(GetPy(), GetPx());
double phi = std::atan(GetPy()/GetPx()); if(phi<0) phi += 2.0*M_PI;
if(GetPx()<0) phi += M_PI; return GetPx() == 0.0 && GetPy() == 0.0 ? 0.0 : phi;
else if(GetPy()<0) phi += 2.0*M_PI;
return phi;
}; };
inline double GetInvMass() const {return std::sqrt(GetE()*GetE() - GetP()*GetP());}; inline double GetInvMass() const {return std::sqrt(GetE()*GetE() - GetP()*GetP());};
inline double GetKE() const {return GetE() - GetInvMass();}; inline double GetKE() const {return GetE() - GetInvMass();};
@ -39,6 +38,12 @@ public:
private: private:
void CalcBoostToCM(); void CalcBoostToCM();
inline double Atan2(double y, double x) const {
if(x != 0) return std::atan2(y, x);
else if( y > 0 ) return M_PI/2.0;
else if( y < 0 ) return -M_PI/2.0;
else return 0.0;
};
double m_data[4]; double m_data[4];
double m_boost[3]; double m_boost[3];

69
include/GRotation.h Normal file
View File

@ -0,0 +1,69 @@
#ifndef GROTATION_H
#define GROTATION_H
#include "G3Vec.h"
class GXRotation {
public:
GXRotation();
GXRotation(double ang);
~GXRotation();
G3Vec Rotate(const G3Vec& vector);
inline void SetAngle(double ang) { m_angle = ang; GenerateMatrix(); };
inline GXRotation GetInverse() { return GXRotation(-m_angle); };
inline G3Vec operator*(const G3Vec& vector) {
double x = m_matrix[0][0]*vector[0] + m_matrix[0][1]*vector[1] + m_matrix[0][2]*vector[2];
double y = m_matrix[1][0]*vector[0] + m_matrix[1][1]*vector[1] + m_matrix[1][2]*vector[2];
double z = m_matrix[2][0]*vector[0] + m_matrix[2][1]*vector[1] + m_matrix[2][2]*vector[2];
return G3Vec(x, y, z);
};
private:
void GenerateMatrix();
double m_angle;
double m_matrix[3][3];
};
class GYRotation {
public:
GYRotation();
GYRotation(double ang);
~GYRotation();
G3Vec Rotate(const G3Vec& vector);
inline void SetAngle(double ang) { m_angle = ang; GenerateMatrix(); };
inline GYRotation GetInverse() { return GYRotation(-m_angle); };
inline G3Vec operator*(const G3Vec& vector) {
double x = m_matrix[0][0]*vector[0] + m_matrix[0][1]*vector[1] + m_matrix[0][2]*vector[2];
double y = m_matrix[1][0]*vector[0] + m_matrix[1][1]*vector[1] + m_matrix[1][2]*vector[2];
double z = m_matrix[2][0]*vector[0] + m_matrix[2][1]*vector[1] + m_matrix[2][2]*vector[2];
return G3Vec(x, y, z);
};
private:
void GenerateMatrix();
double m_angle;
double m_matrix[3][3];
};
class GZRotation {
public:
GZRotation();
GZRotation(double ang);
~GZRotation();
G3Vec Rotate(const G3Vec& vector);
inline void SetAngle(double ang) { m_angle = ang; GenerateMatrix();};
inline GZRotation GetInverse() { return GZRotation(-m_angle); };
inline G3Vec operator*(const G3Vec& vector) {
double x = m_matrix[0][0]*vector[0] + m_matrix[0][1]*vector[1] + m_matrix[0][2]*vector[2];
double y = m_matrix[1][0]*vector[0] + m_matrix[1][1]*vector[1] + m_matrix[1][2]*vector[2];
double z = m_matrix[2][0]*vector[0] + m_matrix[2][1]*vector[1] + m_matrix[2][2]*vector[2];
return G3Vec(x, y, z);
};
private:
void GenerateMatrix();
double m_angle;
double m_matrix[3][3];
};
#endif

170
include/SabreDetector.h
View File

@ -1,136 +1,68 @@
#ifndef __SABREDETECTOR_H #ifndef SABREDETECTOR_H
#define __SABREDETECTOR_H #define SABREDETECTOR_H
#include <TRandom3.h> #include <vector>
#include <cmath>
struct CartCoords { #include "G3Vec.h"
#include "GRotation.h"
double x; class SabreDetector {
double y; public:
double z;
double operator[](int i) { //Overloaded for compatibility with Get_Cart. Only for access SabreDetector();
switch(i) { SabreDetector(double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist=0, double ydist=0);
case(0): return this->x; ~SabreDetector();
case(1): return this->y; inline G3Vec GetRingFlatCoords(int ch, int corner) { return CheckRingLocation(ch, corner) ? m_ringCoords_flat[ch][corner] : G3Vec(); };
case(2): return this->z; inline G3Vec GetWedgeFlatCoords(int ch, int corner) { return CheckWedgeLocation(ch, corner) ? m_wedgeCoords_flat[ch][corner] : G3Vec(); };
default: return 0; inline G3Vec GetRingTiltCoords(int ch, int corner) { return CheckRingLocation(ch, corner) ? m_ringCoords_tilt[ch][corner] : G3Vec(); };
} inline G3Vec GetWedgeTiltCoords(int ch, int corner) { return CheckWedgeLocation(ch, corner) ? m_wedgeCoords_tilt[ch][corner] : G3Vec(); };
} G3Vec GetTrajectoryCoordinates(double theta, double phi);
G3Vec GetHitCoordinates(int ringch, int wedgech);
double GetTheta(); inline int GetNumberOfWedges() { return m_nWedges; };
double GetR(); inline int GetNumberOfRings() { return m_nRings; };
double GetPhi(); inline double GetInnerRadius() { return m_Rinner; };
double GetDetectorPhi(); inline double GetOuterRadius() { return m_Router; };
}; inline double GetPhiCentral() { return m_phiCentral; };
inline double GetTiltAngle() { return m_tilt; };
class SabreDetGeometry { inline G3Vec GetTranslation() { return m_translation; };
public:
SabreDetGeometry(double iRinner_flat,
double iRouter_flat,
double ideltaPhi_flat,
double ibeamPhi_central,
double itiltFromVertical,
double idistFromTarget,
double xoffset=0,
double yoffset=0);
~SabreDetGeometry();
double Get_Ring_Flat_X(int ch, int corner);
double Get_Ring_Flat_Y(int ch, int corner);
double Get_Ring_Flat_Z(int ch, int corner);
double Get_Ring_Flat_R(int ch, int corner);
double Get_Ring_Flat_Theta(int ch, int corner);
double Get_Ring_Flat_Phi(int ch, int corner);
double Get_Wedge_Flat_X(int ch, int corner);
double Get_Wedge_Flat_Y(int ch, int corner);
double Get_Wedge_Flat_Z(int ch, int corner);
double Get_Wedge_Flat_R(int ch, int corner);
double Get_Wedge_Flat_Theta(int ch, int corner);
double Get_Wedge_Flat_Phi(int ch, int corner);
double Get_Ring_Tilted_X(int ch, int corner);
double Get_Ring_Tilted_Y(int ch, int corner);
double Get_Ring_Tilted_Z(int ch, int corner);
double Get_Ring_Tilted_R(int ch, int corner);
double Get_Ring_Tilted_Theta(int ch, int corner);
double Get_Ring_Tilted_Phi(int ch, int corner);
double Get_Wedge_Tilted_X(int ch, int corner);
double Get_Wedge_Tilted_Y(int ch, int corner);
double Get_Wedge_Tilted_Z(int ch, int corner);
double Get_Wedge_Tilted_R(int ch, int corner);
double Get_Wedge_Tilted_Theta(int ch, int corner);
double Get_Wedge_Tilted_Phi(int ch, int corner);
bool IsInside(double theta, double phi);
CartCoords GetTrajectoryCoordinates(double theta, double phi);
void Recenter(double x, double y);
int NumRings(); private:
int NumWedges();
/*** Determine coordinates of the hit (ringch, wedgech) ***/ static constexpr int m_nRings = 16;
CartCoords GetCoordinates(int ringch, int wedgech); static constexpr int m_nWedges = 8;
double GetScatteringAngle(int ringch, int wedgech); //shortcut to Scattering angle static constexpr double deg2rad = M_PI/180.0;
static constexpr double POSITION_TOL = 0.0001;
static constexpr double ANGULAR_TOL = M_PI/180.0;
private: void CalculateCorners();
inline G3Vec TransformToTiltedFrame(G3Vec& vector) { return (m_ZRot*(m_YRot*vector)) + m_translation; };
const int NUMRINGS; inline bool CheckRingChannel(int ch) { return (ch<m_nRings && ch>=0) ? true : false; };
const int NUMWEDGES; inline bool CheckWedgeChannel(int ch) { return (ch<m_nWedges && ch >=0) ? true : false; };
inline bool CheckCorner(int corner) { return (corner < 4 && corner >=0) ? true : false; };
inline bool CheckRingLocation(int ch, int corner) { return CheckRingChannel(ch) && CheckCorner(corner); };
inline bool CheckWedgeLocation(int ch, int corner) { return CheckWedgeChannel(ch) && CheckCorner(corner); };
//detector corners inline bool CheckPositionEqual(double val1,double val2) { return fabs(val1-val2) > POSITION_TOL ? false : true; };
int rbc; //ring bottom channel inline bool CheckAngleEqual(double val1,double val2) { return fabs(val1-val2) > ANGULAR_TOL ? false : true; };
int rtc; //ring top channel
int wrc; //wedge right channel
int wlc; //wedge left channel
double Rinner_flat; inline bool IsInside(double r, double phi) {
double Router_flat; double phi_1 = m_deltaPhi_flat/2.0;
double deltaR_flat; double phi_2 = M_PI*2.0 - m_deltaPhi_flat/2.0;
double deltaR_flat_ring; return (((r > m_Rinner && r < m_Router) || CheckPositionEqual(r, m_Rinner) || CheckPositionEqual(r, m_Router)) && (phi > phi_2 || phi < phi_1 || CheckAngleEqual(phi, phi_1) || CheckAngleEqual(phi, phi_2)));
};
double deltaPhi_flat; double m_Router, m_Rinner, m_deltaPhi_flat, m_phiCentral, m_tilt;
double deltaPhi_flat_wedge; G3Vec m_translation;
GYRotation m_YRot;
GZRotation m_ZRot;
double m_deltaR_flat, m_deltaR_flat_ring;
double beamPhi_central; std::vector<std::vector<G3Vec>> m_ringCoords_flat, m_wedgeCoords_flat;
double tiltFromVertical; std::vector<std::vector<G3Vec>> m_ringCoords_tilt, m_wedgeCoords_tilt;
double ZdistFromTarget;
double XdistFromTarget;
double YdistFromTarget;
TRandom3* random;
//default storage is cartesian
//0=x, 1=y, 2=z
CartCoords **ringch_flat_cart;
CartCoords **wedgech_flat_cart;
CartCoords **ringch_tilted_cart;
CartCoords **wedgech_tilted_cart;
double Get_Cart(int fot, int row, int ch, int corner, int cart);
//fot = flat (0) or tilted (1)
//row = ring (0) or wedge (1)
bool CheckRingChannel(int);
bool CheckWedgeChannel(int);
bool CheckCorner(int);
bool CheckBothRing(int,int);
bool CheckBothWedge(int,int);
/*** Perform transformation for arbitrary point on plane ***/
double **XRotMatrix; //rotation about x-axis
double **ZRotMatrix; //rotation about z-axis
double **XRotMatrixInv; //inverse of x-rotation
double **ZRotMatrixInv; //inverse of z-rotation
CartCoords TransformVector(CartCoords vector);
CartCoords InverseTransformVector(CartCoords vector);
}; };

16
include/SabreEfficiency.h
View File

@ -16,7 +16,9 @@ private:
void RunDecay(const char*); void RunDecay(const char*);
int m_rxn_type; int m_rxn_type;
std::vector<SabreDetGeometry> detectors; std::vector<SabreDetector> detectors;
std::vector<double> ringxs, ringys, ringzs;
std::vector<double> wedgexs, wedgeys, wedgezs;
//Sabre constants //Sabre constants
@ -24,13 +26,13 @@ private:
const double OUTER_R = 0.1351; const double OUTER_R = 0.1351;
const double TILT = 40.0; const double TILT = 40.0;
//const double DIST_2_TARG = 0.14549; //const double DIST_2_TARG = 0.14549;
const double DIST_2_TARG = 0.1245; const double DIST_2_TARG = -0.1245;
const double PHI_COVERAGE = 54.4; //delta phi for each det const double PHI_COVERAGE = 54.4; //delta phi for each det
const double PHI0 = 36.0; //center phi values for each det in array const double PHI0 = 234.0; //center phi values for each det in array
const double PHI1 = 108.0; //# is equal to detID in channel map const double PHI1 = 162.0; //# is equal to detID in channel map
const double PHI2 = 324.0; const double PHI2 = 306.0;
const double PHI3 = 252.0; const double PHI3 = 18.0;
const double PHI4 = 180.0; const double PHI4 = 90.0;
const double DEG2RAD = M_PI/180.0; const double DEG2RAD = M_PI/180.0;
const double ENERGY_THRESHOLD = 0.1; //in MeV const double ENERGY_THRESHOLD = 0.1; //in MeV

21
input.txt
View File

@ -1,31 +1,32 @@
----------Data Information---------- ----------Data Information----------
OutputFile: /Volumes/LaCie/test_newkine.root OutputFile: /media/gordon/b6414c35-ec1f-4fc1-83bc-a6b68ca4325a/gwm17/test_newkine.root
SaveTree: yes SaveTree: yes
SavePlots: yes SavePlots: yes
----------Reaction Information---------- ----------Reaction Information----------
ReactionType: 0 ReactionType: 2
Z A (order is target, projectile, ejectile, break1, break3) Z A (order is target, projectile, ejectile, break1, break3)
5 9 5 10
0 0 2 3
1 1 2 4
1 2
----------Target Information---------- ----------Target Information----------
Name: test_targ Name: test_targ
Layers: 2 Layers: 2
~Layer1 ~Layer1
Thickness(ug/cm^2): 0 Thickness(ug/cm^2): 10
Z A Stoich Z A Stoich
6 12 1 6 12 1
0 0
~ ~
~Layer2 ~Layer2
Thickness(ug/cm^2): 0 Thickness(ug/cm^2): 80
Z A Stoich Z A Stoich
5 9 1 5 10 1
0 0
~ ~
----------Sampling Information---------- ----------Sampling Information----------
NumberOfSamples: 1000000 NumberOfSamples: 10000
BeamMeanEnergy(MeV): 24 BeamEnergySigma(MeV): 0.001 BeamMeanEnergy(MeV): 24 BeamEnergySigma(MeV): 0.001
EjectileThetaMin(deg): 20.0 EjectileThetaMax(deg): 20.0 EjectileThetaMin(deg): 20.0 EjectileThetaMax(deg): 20.0
ResidualExMean(MeV): 16.708 ResidualExSigma(MeV): 0.038 ResidualExMean(MeV): 16.8 ResidualExSigma(MeV): 0.038
-------------------------------------- --------------------------------------

8
makefile
View File

@ -13,16 +13,13 @@ BINDIR=./bin
SRC=$(wildcard $(SRCDIR)/*.cpp) SRC=$(wildcard $(SRCDIR)/*.cpp)
OBJS=$(SRC:$(SRCDIR)/%.cpp=$(OBJDIR)/%.o) OBJS=$(SRC:$(SRCDIR)/%.cpp=$(OBJDIR)/%.o)
TPOBJ=objs/testplots.o
TPEXE=tp
DICTOBJ=$(OBJDIR)/kinematics_dict.o DICTOBJ=$(OBJDIR)/kinematics_dict.o
DICTSRC=$(SRCDIR)/kinematics_dict.cxx DICTSRC=$(SRCDIR)/kinematics_dict.cxx
DICT_PAGES=$(INCLDIR)/Kinematics.h $(INCLDIR)/LinkDef_Kinematics.h DICT_PAGES=$(INCLDIR)/Kinematics.h $(INCLDIR)/LinkDef_Kinematics.h
EXE=$(BINDIR)/kinematics EXE=$(BINDIR)/kinematics
CLEANUP=$(EXE) $(OBJS) $(DICTOBJ) $(DICTSRC) $(TPOBJ) CLEANUP=$(EXE) $(OBJS) $(DICTOBJ) $(DICTSRC)
.PHONY: all clean .PHONY: all clean
@ -31,9 +28,6 @@ all: $(EXE) $(TPEXE)
$(EXE): $(DICTOBJ) $(OBJS) $(EXE): $(DICTOBJ) $(OBJS)
$(CC) $(CFLAGS) $^ -o $@ $(LDFLAGS) $(CC) $(CFLAGS) $^ -o $@ $(LDFLAGS)
$(TPEXE): ./objs/SabreDetector.o $(TPOBJ)
$(CC) $(CFLAGS) $^ -o $@ $(LDFLAGS)
$(DICTOBJ): $(DICTSRC) $(DICTOBJ): $(DICTSRC)
$(CC) $(CFLAGS) $(CPPFLAGS) -I $(ROOTINCLDIR) -c $^ -o $@ $(CC) $(CFLAGS) $(CPPFLAGS) -I $(ROOTINCLDIR) -c $^ -o $@
mv $(SRCDIR)/*.pcm $(BINDIR) mv $(SRCDIR)/*.pcm $(BINDIR)

3
objs/.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
###keep only the directory, not the contents###
*
!.gitignore

BIN
objs/EnergyLoss.o
View File

Binary file not shown.

BIN
objs/G4Vec.o
View File

Binary file not shown.

BIN
objs/Kinematics.o
View File

Binary file not shown.

BIN
objs/LayeredTarget.o
View File

Binary file not shown.

BIN
objs/MassLookup.o
View File

Binary file not shown.

BIN
objs/Nucleus.o
View File

Binary file not shown.

BIN
objs/Plotter.o
View File

Binary file not shown.

BIN
objs/Reaction.o
View File

Binary file not shown.

BIN
objs/ReactionSystem.o
View File

Binary file not shown.

BIN
objs/SabreDetector.o
View File

Binary file not shown.

BIN
objs/SabreEfficiency.o
View File

Binary file not shown.

BIN
objs/Target.o
View File

Binary file not shown.

BIN
objs/ThreeStepSystem.o
View File

Binary file not shown.

BIN
objs/TwoStepSystem.o
View File

Binary file not shown.

BIN
objs/kinematics_dict.o
View File

Binary file not shown.

BIN
objs/main.o
View File

Binary file not shown.

BIN
objs/testplots.o
View File

Binary file not shown.

1
src/EnergyLoss.cpp
View File

@ -9,6 +9,7 @@ Written by G.W. McCann Aug. 2020
*/ */
#include "Eloss_Tables.h"
#include "EnergyLoss.h" #include "EnergyLoss.h"
#include "KinematicsExceptions.h" #include "KinematicsExceptions.h"
#include <iostream> #include <iostream>

36
src/G3Vec.cpp Normal file
View File

@ -0,0 +1,36 @@
#include "G3Vec.h"
G3Vec::G3Vec() {
m_data[0] = 0.;
m_data[1] = 0.;
m_data[2] = 0.;
}
G3Vec::G3Vec(double x, double y, double z) {
m_data[0] = x;
m_data[1] = y;
m_data[2] = z;
}
G3Vec::~G3Vec() {}
void G3Vec::SetVectorCartesian(double x, double y, double z) {
m_data[0] = x;
m_data[1] = y;
m_data[2] = z;
}
void G3Vec::SetVectorSpherical(double r, double theta, double phi) {
m_data[0] = r*std::cos(phi)*std::sin(theta);
m_data[1] = r*std::sin(phi)*std::sin(theta);
m_data[2] = r*std::cos(theta);
}
double G3Vec::Dot(const G3Vec& rhs) const {
return GetX()*rhs.GetX() + GetY()*rhs.GetY() + GetZ()*rhs.GetZ();
}
//Unimplemented
G3Vec G3Vec::Cross(const G3Vec& rhs) const {
return G3Vec(0.,0.,0.);
}

62
src/GRotation.cpp Normal file
View File

@ -0,0 +1,62 @@
#include "GRotation.h"
GXRotation::GXRotation() :
m_angle(0)
{
GenerateMatrix();
}
GXRotation::GXRotation(double angle) :
m_angle(angle)
{
GenerateMatrix();
}
GXRotation::~GXRotation() {}
void GXRotation::GenerateMatrix() {
m_matrix[0][0] = 1.0; m_matrix[0][1] = 0.0; m_matrix[0][2] = 0.0;
m_matrix[1][0] = 0.0; m_matrix[1][1] = std::cos(m_angle); m_matrix[1][2] = -std::sin(m_angle);
m_matrix[2][0] = 0.0; m_matrix[2][1] = std::sin(m_angle); m_matrix[2][2] = std::cos(m_angle);
}
GYRotation::GYRotation() :
m_angle(0)
{
GenerateMatrix();
}
GYRotation::GYRotation(double angle) :
m_angle(angle)
{
GenerateMatrix();
}
GYRotation::~GYRotation() {}
void GYRotation::GenerateMatrix() {
m_matrix[0][0] = std::cos(m_angle); m_matrix[0][1] = 0.0; m_matrix[0][2] = -std::sin(m_angle);
m_matrix[1][0] = 0.0; m_matrix[1][1] = 1.0; m_matrix[1][2] = 0.0;
m_matrix[2][0] = std::sin(m_angle); m_matrix[2][1] = 0.0; m_matrix[2][2] = std::cos(m_angle);
}
GZRotation::GZRotation() :
m_angle(0)
{
GenerateMatrix();
}
GZRotation::GZRotation(double angle) :
m_angle(angle)
{
GenerateMatrix();
}
GZRotation::~GZRotation() {}
void GZRotation::GenerateMatrix() {
m_matrix[0][0] = std::cos(m_angle); m_matrix[0][1] = -std::sin(m_angle); m_matrix[0][2] = 0.0;
m_matrix[1][0] = std::sin(m_angle); m_matrix[1][1] = std::cos(m_angle); m_matrix[1][2] = 0.0;
m_matrix[2][0] = 0.0; m_matrix[2][1] = 0.0; m_matrix[2][2] = 1.0;
}

1
src/Plotter.cpp
View File

@ -1,4 +1,5 @@
#include "Plotter.h" #include "Plotter.h"
#include <iostream>
Plotter::Plotter() : Plotter::Plotter() :
table(new THashTable()) table(new THashTable())

762
src/SabreDetector.cpp
View File

@ -1,10 +1,6 @@
#include "SabreDetector.h" #include "SabreDetector.h"
#include <cmath>
#include <iostream> #include <iostream>
using namespace std;
/* /*
Distances in meters, angles in radians. Distances in meters, angles in radians.
@ -13,11 +9,11 @@ using namespace std;
0---------------------1 0---------------------1
| | | |
| | y | | x
| | ^ | | <-----
| | | | | |
| | | | | |
3---------------------2 -----> x 3---------------------2 y
(z is hence positive along beam direction) (z is hence positive along beam direction)
The channel numbers, also as looking back from target pov, are: The channel numbers, also as looking back from target pov, are:
@ -44,610 +40,182 @@ using namespace std;
| | | | | | | | | | | |
| | | | | | | | | | | |
***note that these are for ARRAY storage, and may not necessarily
correspond to the PHYSICAL channels; this will need to be taken
into account if used in actual data analysis
-- kgh, March 2020 >> Note that the detector starts centered on the x-axis (central phi = 0) untilted, and then is rotated to wherever the frick
it is supposed to go; phi = 90 is centered on y axis, pointing down towards the bottom of the scattering chamber
-- gwm, Dec 2020; based on the og code from kgh
*/ */
SabreDetGeometry::SabreDetGeometry(double iRinner_flat,
double iRouter_flat,
double ideltaPhi_flat,
double ibeamPhi_central,
double itiltFromVertical,
double idistFromTarget,
double xoffset,
double yoffset) :
NUMRINGS(16), NUMWEDGES(8) {
rbc = 0; //ring bottom channel SabreDetector::SabreDetector() :
rtc = NUMRINGS-1; //ring top channel m_Router(0.1351), m_Rinner(0.0326), m_deltaPhi_flat(54.4*deg2rad), m_phiCentral(0.0), m_tilt(0.0), m_translation(0.,0.,0.)
wrc = 0; //wedge right channel {
wlc = NUMWEDGES-1; //wedge left channel m_YRot.SetAngle(m_tilt);
m_ZRot.SetAngle(m_phiCentral);
Rinner_flat = iRinner_flat; m_ringCoords_flat.resize(m_nRings);
Router_flat = iRouter_flat; m_ringCoords_tilt.resize(m_nRings);
deltaR_flat = iRouter_flat - iRinner_flat; m_wedgeCoords_flat.resize(m_nWedges);
deltaR_flat_ring = deltaR_flat/NUMRINGS; m_wedgeCoords_tilt.resize(m_nWedges);
for(int i=0; i<m_nRings; i++) {
m_ringCoords_flat[i].resize(4);
m_ringCoords_tilt[i].resize(4);
}
for(int i=0; i<m_nWedges; i++) {
m_wedgeCoords_flat[i].resize(4);
m_wedgeCoords_tilt[i].resize(4);
}
deltaPhi_flat = ideltaPhi_flat; m_deltaR_flat = m_Router - m_Rinner;
deltaPhi_flat_wedge = deltaPhi_flat/NUMWEDGES; m_deltaR_flat_ring = m_deltaR_flat/m_nRings;
beamPhi_central = ibeamPhi_central;
tiltFromVertical = itiltFromVertical;
//this distance from target is from the CENTRAL point of
//the detector, i.e. the center of the "circle" which the
//detector forms. NOTE: these are all assumed to be negative!
ZdistFromTarget = idistFromTarget;
XdistFromTarget = xoffset;
YdistFromTarget = yoffset;
random = new TRandom3();
random->SetSeed();
ringch_flat_cart = new CartCoords*[NUMRINGS];
wedgech_flat_cart = new CartCoords*[NUMWEDGES];
ringch_tilted_cart = new CartCoords*[NUMRINGS];
wedgech_tilted_cart = new CartCoords*[NUMWEDGES];
for (int i=0; i<NUMRINGS; i++) {
ringch_flat_cart[i] = new CartCoords[4];
ringch_tilted_cart[i] = new CartCoords[4];
}
for (int i=0; i<NUMWEDGES; i++) {
wedgech_flat_cart[i] = new CartCoords[4];
wedgech_tilted_cart[i] = new CartCoords[4];
}
//Define rotation matricies
XRotMatrix = new double*[3];
ZRotMatrix = new double*[3];
XRotMatrixInv = new double*[3];
ZRotMatrixInv = new double*[3];
for (int i=0; i<3; i++) {
XRotMatrix[i] = new double[3];
ZRotMatrix[i] = new double[3];
XRotMatrixInv[i] = new double[3];
ZRotMatrixInv[i] = new double[3];
}
XRotMatrix[0][0] = 1;
XRotMatrix[0][1] = 0; XRotMatrix[0][2] = 0; XRotMatrix[1][0] = 0; XRotMatrix[2][0] = 0;
XRotMatrix[1][1] = cos(tiltFromVertical); XRotMatrix[1][2] = -sin(tiltFromVertical);
XRotMatrix[2][1] = sin(tiltFromVertical); XRotMatrix[2][2] = cos(tiltFromVertical);
ZRotMatrix[0][0] = cos(beamPhi_central); ZRotMatrix[0][1] = -sin(beamPhi_central);
ZRotMatrix[1][0] = sin(beamPhi_central); ZRotMatrix[1][1] = cos(beamPhi_central);
ZRotMatrix[2][2] = 1;
ZRotMatrix[2][0] = 0; ZRotMatrix[2][1] = 0; ZRotMatrix[0][2] = 0; ZRotMatrix[1][2] = 0;
//Inverse of a rotation is its transpose
XRotMatrixInv[0][0] = XRotMatrix[0][0]; XRotMatrixInv[0][1] = XRotMatrix[1][0]; XRotMatrixInv[0][2] = XRotMatrix[2][0];
XRotMatrixInv[1][0] = XRotMatrix[0][1]; XRotMatrixInv[1][1] = XRotMatrix[1][1]; XRotMatrixInv[1][2] = XRotMatrix[2][1];
XRotMatrixInv[2][0] = XRotMatrix[0][2]; XRotMatrixInv[2][1] = XRotMatrix[1][2]; XRotMatrixInv[2][2] = XRotMatrix[2][2];
ZRotMatrixInv[0][0] = ZRotMatrix[0][0]; ZRotMatrixInv[0][1] = ZRotMatrix[1][0]; ZRotMatrixInv[0][2] = ZRotMatrix[2][0];
ZRotMatrixInv[1][0] = ZRotMatrix[0][1]; ZRotMatrixInv[1][1] = ZRotMatrix[1][1]; ZRotMatrixInv[1][2] = ZRotMatrix[2][1];
ZRotMatrixInv[2][0] = ZRotMatrix[0][2]; ZRotMatrixInv[2][1] = ZRotMatrix[1][2]; ZRotMatrixInv[2][2] = ZRotMatrix[2][2];
//
/* SETTING FLAT GEOMETRY */
//RINGS:
for (int rch=0; rch<NUMRINGS; rch++) {
//top left:
ringch_flat_cart[rch][0].x = (Rinner_flat + (rch+1)*deltaR_flat_ring)*sin(deltaPhi_flat/2); //x
ringch_flat_cart[rch][0].y = (Rinner_flat + (rch+1)*deltaR_flat_ring)*cos(deltaPhi_flat/2); //y
ringch_flat_cart[rch][0].z = 0; //z
//top right:
ringch_flat_cart[rch][1].x = (Rinner_flat + (rch+1)*deltaR_flat_ring)*sin(-deltaPhi_flat/2); //x
ringch_flat_cart[rch][1].y = (Rinner_flat + (rch+1)*deltaR_flat_ring)*cos(-deltaPhi_flat/2); //y
ringch_flat_cart[rch][1].z = 0; //z
//bottom right:
ringch_flat_cart[rch][2].x = (Rinner_flat + rch*deltaR_flat_ring)*sin(-deltaPhi_flat/2); //x
ringch_flat_cart[rch][2].y = (Rinner_flat + rch*deltaR_flat_ring)*cos(-deltaPhi_flat/2); //y
ringch_flat_cart[rch][2].z = 0; //z
//bottom left:
ringch_flat_cart[rch][3].x = (Rinner_flat + rch*deltaR_flat_ring)*sin(deltaPhi_flat/2); //x
ringch_flat_cart[rch][3].y = (Rinner_flat + rch*deltaR_flat_ring)*cos(deltaPhi_flat/2); //y
ringch_flat_cart[rch][3].z = 0; //z
}
//WEDGES:
for (int wch=0; wch<NUMWEDGES; wch++) {
//top left:
wedgech_flat_cart[wch][0].x = Router_flat*sin(-deltaPhi_flat/2 + (wch+1)*deltaPhi_flat_wedge); //x
wedgech_flat_cart[wch][0].y = Router_flat*cos(-deltaPhi_flat/2 + (wch+1)*deltaPhi_flat_wedge); //y
wedgech_flat_cart[wch][0].z = 0; //z
//top right:
wedgech_flat_cart[wch][1].x = Router_flat*sin(-deltaPhi_flat/2 + wch*deltaPhi_flat_wedge); //x
wedgech_flat_cart[wch][1].y = Router_flat*cos(-deltaPhi_flat/2 + wch*deltaPhi_flat_wedge); //y
wedgech_flat_cart[wch][1].z = 0; //z
//bottom right:
wedgech_flat_cart[wch][2].x = Rinner_flat*sin(-deltaPhi_flat/2 + wch*deltaPhi_flat_wedge); //x
wedgech_flat_cart[wch][2].y = Rinner_flat*cos(-deltaPhi_flat/2 + wch*deltaPhi_flat_wedge); //y
wedgech_flat_cart[wch][2].z = 0; //z
//bottom left:
wedgech_flat_cart[wch][3].x = Rinner_flat*sin(-deltaPhi_flat/2 + (wch+1)*deltaPhi_flat_wedge); //x
wedgech_flat_cart[wch][3].y = Rinner_flat*cos(-deltaPhi_flat/2 + (wch+1)*deltaPhi_flat_wedge); //y
wedgech_flat_cart[wch][3].z = 0; //z
}
//now tilting, rotating, and translating
//RINGS:
for (int rch=0; rch<NUMRINGS; rch++) {
for (int c=0; c<4; c++) {
ringch_tilted_cart[rch][c] = TransformVector(ringch_flat_cart[rch][c]);
}
}
//WEDGES:
for (int wch=0; wch<NUMWEDGES; wch++) {
for (int c=0; c<4; c++) {
wedgech_tilted_cart[wch][c] = TransformVector(wedgech_flat_cart[wch][c]);
}
}
CalculateCorners();
} }
SabreDetGeometry::~SabreDetGeometry() { SabreDetector::SabreDetector(double Rin, double Rout, double deltaPhi_flat, double phiCentral, double tiltFromVert, double zdist, double xdist, double ydist) :
m_Router(Rout), m_Rinner(Rin), m_deltaPhi_flat(deltaPhi_flat), m_phiCentral(phiCentral), m_tilt(tiltFromVert), m_translation(xdist, ydist, zdist)
{
std::cout<<m_tilt<<std::endl;
m_YRot.SetAngle(m_tilt);
m_ZRot.SetAngle(m_phiCentral);
for (int i=0; i<NUMRINGS; i++) { m_ringCoords_flat.resize(m_nRings);
delete [] ringch_flat_cart[i]; m_ringCoords_tilt.resize(m_nRings);
delete [] ringch_tilted_cart[i]; m_wedgeCoords_flat.resize(m_nWedges);
} m_wedgeCoords_tilt.resize(m_nWedges);
delete [] ringch_flat_cart; for(int i=0; i<m_nRings; i++) {
delete [] ringch_tilted_cart; m_ringCoords_flat[i].resize(4);
m_ringCoords_tilt[i].resize(4);
for (int i=0; i<NUMWEDGES; i++) { }
delete [] wedgech_flat_cart[i]; for(int i=0; i<m_nWedges; i++) {
delete [] wedgech_tilted_cart[i]; m_wedgeCoords_flat[i].resize(4);
} m_wedgeCoords_tilt[i].resize(4);
delete [] wedgech_flat_cart; }
delete [] wedgech_tilted_cart;
for(int i=0; i<3; i++) {
delete [] XRotMatrix[i];
delete [] ZRotMatrix[i];
delete [] XRotMatrixInv[i];
delete [] ZRotMatrixInv[i];
}
delete [] XRotMatrix;
delete [] ZRotMatrix;
delete random;
CalculateCorners();
} }
double SabreDetGeometry::Get_Ring_Flat_X(int ch, int corner) { SabreDetector::~SabreDetector() {}
return CheckBothRing(ch,corner) ? Get_Cart(0,0,ch,corner,0) : 0.;
} void SabreDetector::CalculateCorners() {
double SabreDetGeometry::Get_Ring_Flat_Y(int ch, int corner) {
return CheckBothRing(ch,corner) ? Get_Cart(0,0,ch,corner,1) : 0.; double deltaPhi_per_wedge = m_deltaPhi_flat/double(m_nWedges);
} double deltaR_per_ring = (m_Router - m_Rinner)/double(m_nRings);
double SabreDetGeometry::Get_Ring_Flat_Z(int ch, int corner) {
return CheckBothRing(ch,corner) ? Get_Cart(0,0,ch,corner,2) : 0.; double x0, x1, x2, x3;
} double y0, y1, y2, y3;
double SabreDetGeometry::Get_Ring_Flat_R(int ch, int corner) { double z0, z1, z2, z3;
return CheckBothRing(ch,corner) ?
sqrt(pow(Get_Cart(0,0,ch,corner,0),2) + //Generate flat rings
pow(Get_Cart(0,0,ch,corner,1),2) + for(int i=0; i<m_nRings; i++) {
pow(Get_Cart(0,0,ch,corner,2),2)) x0 = (m_Rinner + deltaR_per_ring*(i+1))*std::cos(-m_deltaPhi_flat/2.0);
: 0.; y0 = (m_Rinner + deltaR_per_ring*(i+1))*std::sin(-m_deltaPhi_flat/2.0);
} z0 = 0.0;
double SabreDetGeometry::Get_Ring_Flat_Theta(int ch, int corner) { m_ringCoords_flat[i][0].SetVectorCartesian(x0, y0, z0);
return CheckBothRing(ch,corner) ?
acos(Get_Ring_Flat_Z(ch,corner)/Get_Ring_Flat_R(ch,corner)) x1 = (m_Rinner + deltaR_per_ring*(i))*std::cos(-m_deltaPhi_flat/2.0);
: 0.; y1 = (m_Rinner + deltaR_per_ring*(i))*std::sin(-m_deltaPhi_flat/2.0);
} z1 = 0.0;
//GWM: Note this now returns in normal phi from x-axis format m_ringCoords_flat[i][1].SetVectorCartesian(x1, y1, z1);
double SabreDetGeometry::Get_Ring_Flat_Phi(int ch, int corner) {
if (!(CheckBothRing(ch,corner))) return 0.; x2 = (m_Rinner + deltaR_per_ring*(i))*std::cos(m_deltaPhi_flat/2.0);
double x = Get_Ring_Flat_X(ch,corner); y2 = (m_Rinner + deltaR_per_ring*(i))*std::sin(m_deltaPhi_flat/2.0);
double y = Get_Ring_Flat_Y(ch,corner); z2 = 0.0;
double phi = std::atan2(y,x); m_ringCoords_flat[i][2].SetVectorCartesian(x2, y2, z2);
if (phi<0) phi += M_PI*2.0;
return phi; x3 = (m_Rinner + deltaR_per_ring*(i+1))*std::cos(m_deltaPhi_flat/2.0);
y3 = (m_Rinner + deltaR_per_ring*(i+1))*std::sin(m_deltaPhi_flat/2.0);
z3 = 0.0;
m_ringCoords_flat[i][3].SetVectorCartesian(x3, y3, z3);
}
//Generate flat wedges
for(int i=0; i<m_nWedges; i++) {
x0 = m_Router * std::cos(-m_deltaPhi_flat/2.0 + i*deltaPhi_per_wedge);
y0 = m_Router * std::sin(-m_deltaPhi_flat/2.0 + i*deltaPhi_per_wedge);
z0 = 0.0;
m_wedgeCoords_flat[i][0].SetVectorCartesian(x0, y0, z0);
x1 = m_Rinner * std::cos(-m_deltaPhi_flat/2.0 + i*deltaPhi_per_wedge);
y1 = m_Rinner * std::sin(-m_deltaPhi_flat/2.0 + i*deltaPhi_per_wedge);
z1 = 0.0;
m_wedgeCoords_flat[i][1].SetVectorCartesian(x1, y1, z1);
x2 = m_Rinner * std::cos(-m_deltaPhi_flat/2.0 + (i+1)*deltaPhi_per_wedge);
y2 = m_Rinner * std::sin(-m_deltaPhi_flat/2.0 + (i+1)*deltaPhi_per_wedge);
z2 = 0.0;
m_wedgeCoords_flat[i][2].SetVectorCartesian(x2, y2, z2);
x3 = m_Router * std::cos(-m_deltaPhi_flat/2.0 + (i+1)*deltaPhi_per_wedge);
y3 = m_Router * std::sin(-m_deltaPhi_flat/2.0 + (i+1)*deltaPhi_per_wedge);
z3 = 0.0;
m_wedgeCoords_flat[i][3].SetVectorCartesian(x3, y3, z3);
}
//Generate tilted rings
for(int i=0; i<m_nRings; i++) {
for(int j=0; j<4; j++) {
m_ringCoords_tilt[i][j] = TransformToTiltedFrame(m_ringCoords_flat[i][j]);
}
}
//Generate tilted wedges
for(int i=0; i<m_nWedges; i++) {
for(int j=0; j<4; j++) {
m_wedgeCoords_tilt[i][j] = TransformToTiltedFrame(m_wedgeCoords_flat[i][j]);
}
}
} }
double SabreDetGeometry::Get_Wedge_Flat_X(int ch, int corner) { //Note: float used for calculations due to lack of precision from sin, cos, and tangent functions
return CheckBothWedge(ch,corner) ? Get_Cart(0,1,ch,corner,0) : 0.; G3Vec SabreDetector::GetTrajectoryCoordinates(double theta, double phi) {
} if(m_translation.GetX() != 0.0 || m_translation.GetY() != 0.0) {
double SabreDetGeometry::Get_Wedge_Flat_Y(int ch, int corner) { return G3Vec();
return CheckBothWedge(ch,corner) ? Get_Cart(0,1,ch,corner,1) : 0.; }
}
double SabreDetGeometry::Get_Wedge_Flat_Z(int ch, int corner) { //Calculate the *potential* phi in the flat detector
return CheckBothWedge(ch,corner) ? Get_Cart(0,1,ch,corner,2) : 0.; double phi_numerator = std::cos(m_tilt)*(std::sin(phi)*std::cos(m_phiCentral) - std::sin(m_phiCentral)*std::cos(phi));
} double phi_denominator = std::cos(m_phiCentral)*std::cos(phi) + std::sin(m_phiCentral)*std::sin(phi);
double SabreDetGeometry::Get_Wedge_Flat_R(int ch, int corner) { double phi_flat = std::atan2(phi_numerator, phi_denominator);
return CheckBothWedge(ch,corner) ? if(phi_flat < 0) phi_flat += M_PI*2.0;
sqrt(pow(Get_Cart(0,1,ch,corner,0),2) +
pow(Get_Cart(0,1,ch,corner,1),2) + //Calculate the *potential* R in the flat detector
pow(Get_Cart(0,1,ch,corner,2),2)) double r_numerator = m_translation.GetZ()*std::cos(phi)*std::sin(theta);
: 0.; double r_denominator = std::cos(phi_flat)*std::cos(m_phiCentral)*std::cos(m_tilt)*std::cos(theta) - std::sin(phi_flat)*std::sin(m_phiCentral)*std::cos(theta) - std::cos(phi_flat)*std::sin(m_tilt)*std::cos(phi)*std::sin(theta);
} double r_flat = r_numerator/r_denominator;
double SabreDetGeometry::Get_Wedge_Flat_Theta(int ch, int corner) {
return CheckBothWedge(ch,corner) ? //Calculate the distance from the origin to the hit on the detector
acos(Get_Wedge_Flat_Z(ch,corner)/Get_Wedge_Flat_R(ch,corner)) double R_to_detector = (r_flat*std::cos(phi_flat)*std::sin(m_tilt) + m_translation.GetZ())/std::cos(theta);
: 0.; double xhit = R_to_detector*std::sin(theta)*std::cos(phi);
} double yhit = R_to_detector*std::sin(theta)*std::sin(phi);
//GWM: Note this now returns in normal phi from x-axis format double zhit = R_to_detector*std::cos(theta);
double SabreDetGeometry::Get_Wedge_Flat_Phi(int ch, int corner) {
if (!(CheckBothWedge(ch,corner))) return false;
double x = Get_Wedge_Flat_X(ch,corner); //Check to see if our flat coords fall inside the flat detector
double y = Get_Wedge_Flat_Y(ch,corner); if(IsInside(r_flat, phi_flat)) {
double phi = std::atan2(y,x); return G3Vec(xhit, yhit, zhit);
if(phi<0) phi += M_PI*2.0; } else {
return phi; return G3Vec();
}
} }
double SabreDetGeometry::Get_Ring_Tilted_X(int ch, int corner) { G3Vec SabreDetector::GetHitCoordinates(int ringch, int wedgech) {
return CheckBothRing(ch,corner) ? Get_Cart(1,0,ch,corner,0) : 0.; if(!CheckRingChannel(ringch) || !CheckWedgeChannel(wedgech)) {
return G3Vec();
}
double deltaR_per_ring = (m_Router - m_Rinner)/double(m_nRings);
double deltaPhi_per_wedge = (m_deltaPhi_flat)/double(m_nWedges);
double r_center = m_Rinner + (0.5+ringch)*deltaR_per_ring;
double phi_center = -m_deltaPhi_flat/2.0 + (0.5+wedgech)*deltaPhi_per_wedge;
double x = r_center*std::cos(phi_center);
double y = r_center*std::sin(phi_center);
double z = 0;
G3Vec hit(x, y, z);
return TransformToTiltedFrame(hit);
} }
double SabreDetGeometry::Get_Ring_Tilted_Y(int ch, int corner) {
return CheckBothRing(ch,corner) ? Get_Cart(1,0,ch,corner,1) : 0.;
}
double SabreDetGeometry::Get_Ring_Tilted_Z(int ch, int corner) {
return CheckBothRing(ch,corner) ? Get_Cart(1,0,ch,corner,2) : 0.;
}
double SabreDetGeometry::Get_Ring_Tilted_R(int ch, int corner) {
return CheckBothRing(ch,corner) ?
sqrt(pow(Get_Cart(1,0,ch,corner,0),2) +
pow(Get_Cart(1,0,ch,corner,1),2) +
pow(Get_Cart(1,0,ch,corner,2),2))
: 0.;
}
double SabreDetGeometry::Get_Ring_Tilted_Theta(int ch, int corner) {
return CheckBothRing(ch,corner) ?
acos(Get_Ring_Tilted_Z(ch,corner)/Get_Ring_Tilted_R(ch,corner))
: 0.;
}
//GWM: Note this now returns in normal phi from x-axis format
double SabreDetGeometry::Get_Ring_Tilted_Phi(int ch, int corner) {
if (!(CheckBothRing(ch,corner))) return false;
double x = Get_Ring_Tilted_X(ch,corner);
double y = Get_Ring_Tilted_Y(ch,corner);
double phi = std::atan2(y,x);
if(phi<0) phi += M_PI*2.0;
return phi;
}
double SabreDetGeometry::Get_Wedge_Tilted_X(int ch, int corner) {
return CheckBothWedge(ch,corner) ? Get_Cart(1,1,ch,corner,0) : 0.;
}
double SabreDetGeometry::Get_Wedge_Tilted_Y(int ch, int corner) {
return CheckBothWedge(ch,corner) ? Get_Cart(1,1,ch,corner,1) : 0.;
}
double SabreDetGeometry::Get_Wedge_Tilted_Z(int ch, int corner) {
return CheckBothWedge(ch,corner) ? Get_Cart(1,1,ch,corner,2) : 0.;
}
double SabreDetGeometry::Get_Wedge_Tilted_R(int ch, int corner) {
return CheckBothWedge(ch,corner) ?
sqrt(pow(Get_Cart(1,1,ch,corner,0),2) +
pow(Get_Cart(1,1,ch,corner,1),2) +
pow(Get_Cart(1,1,ch,corner,2),2))
: 0.;
}
double SabreDetGeometry::Get_Wedge_Tilted_Theta(int ch, int corner) {
return CheckBothWedge(ch,corner) ?
acos(Get_Wedge_Tilted_Z(ch,corner)/Get_Wedge_Tilted_R(ch,corner))
: 0.;
}
//GWM: Note this now returns in normal phi from x-axis format
double SabreDetGeometry::Get_Wedge_Tilted_Phi(int ch, int corner) {
if (!(CheckBothWedge(ch,corner))) return false;
double x = Get_Wedge_Tilted_X(ch,corner);
double y = Get_Wedge_Tilted_Y(ch,corner);
double phi = std::atan2(y, x);
if(phi<0) phi += 2.0*M_PI;
return phi;
}
int SabreDetGeometry::NumRings() {return NUMRINGS;}
int SabreDetGeometry::NumWedges() {return NUMWEDGES;}
double SabreDetGeometry::Get_Cart(int fot, int row, int ch, int corner, int cart) {
//fot = flat (0) or tilted (1)
//row = ring (0) or wedge (1)
int marker = 10*fot + row;
switch (marker) {
case (00): return ringch_flat_cart[ch][corner][cart];
case (01): return wedgech_flat_cart[ch][corner][cart];
case (10): return ringch_tilted_cart[ch][corner][cart];
case (11): return wedgech_tilted_cart[ch][corner][cart];
default: return 0;
}
}
bool SabreDetGeometry::CheckRingChannel(int ich) {
return ((ich >= 0 && ich < NUMRINGS) ? true : false);
}
bool SabreDetGeometry::CheckWedgeChannel(int ich) {
return ((ich >= 0 && ich < NUMWEDGES) ? true: false);
}
bool SabreDetGeometry::CheckCorner(int icn) {
return ((icn >= 0 && icn < 4) ? true : false);
}
bool SabreDetGeometry::CheckBothRing(int ich, int icn) {
return (CheckRingChannel(ich)&&CheckCorner(icn));
}
bool SabreDetGeometry::CheckBothWedge(int ich, int icn) {
return (CheckWedgeChannel(ich)&&CheckCorner(icn));
}
/*Method by which the coordinates of a hit in a wedge/ring pixel are calculated.
*Currently, takes the center of the pixel as the value of a pixel hit, could be altered
*to take a uniformly sampled random point with in the pixel
*/
CartCoords SabreDetGeometry::GetCoordinates(int ringch, int wedgech) {
if(!CheckRingChannel(ringch) || !CheckWedgeChannel(wedgech)) {
CartCoords temp;
temp.x = 0;
temp.y = 0;
temp.z = 0;
return temp;
}
//define pixel by its center (half way between top and bottom radius, halfway between left and right phi)
//EDIT GWM: July 2020 change to randomize the location within the sabre pixel
double rcenter = Rinner_flat + deltaR_flat_ring*(ringch + random->Uniform(0.0, 1.0));
double phi_center = deltaPhi_flat/2.0-(wedgech + random->Uniform(0.0, 1.0))*deltaPhi_flat_wedge;
CartCoords PixelCenter;
PixelCenter.x = rcenter*sin(phi_center);
PixelCenter.y = rcenter*cos(phi_center);
PixelCenter.z = 0;
//return coords in final orientation
return TransformVector(PixelCenter);
}
CartCoords SabreDetGeometry::TransformVector(CartCoords vector) {
CartCoords xrot_vector, xzrot_vector, xzrot_t_vector;
std::cout<<"Starting coords -- x: "<<vector.x<<" y: "<<vector.y<<" z: "<<vector.z<<" r: "<<vector.GetR()<<std::endl;
//rotation about x-axis
xrot_vector.x = vector.x;
xrot_vector.y = XRotMatrix[1][1]*vector.y+XRotMatrix[1][2]*vector.z;
xrot_vector.z = XRotMatrix[2][1]*vector.y+XRotMatrix[2][2]*vector.z;
std::cout<<"after x-rotation -- x: "<<xrot_vector.x<<" y: "<<xrot_vector.y<<" z: "<<xrot_vector.z<<" r: "<<vector.GetR()<<std::endl;
//rotation about z-axis
xzrot_vector.x = ZRotMatrix[0][0]*xrot_vector.x+ZRotMatrix[0][1]*xrot_vector.y;
xzrot_vector.y = ZRotMatrix[1][0]*xrot_vector.x+ZRotMatrix[1][1]*xrot_vector.y;
xzrot_vector.z = ZRotMatrix[2][2]*xrot_vector.z;
std::cout<<"after z-rotation -- x: "<<xzrot_vector.x<<" y: "<<xzrot_vector.y<<" z: "<<xzrot_vector.z<<" r: "<<vector.GetR()<<std::endl;
//translate as needed
xzrot_t_vector.x = xzrot_vector.x-XdistFromTarget;
xzrot_t_vector.y = xzrot_vector.y-YdistFromTarget;
xzrot_t_vector.z = xzrot_vector.z-ZdistFromTarget;
std::cout<<"after z-translation -- x: "<<xzrot_t_vector.x<<" y: "<<xzrot_t_vector.y<<" z: "<<xzrot_t_vector.z<<std::endl;
return xzrot_t_vector;
}
CartCoords SabreDetGeometry::InverseTransformVector(CartCoords vector) {
CartCoords it_vector, izrot_vector, ixzrot_vector;
//Inverse of translation
/*it_vector.x = vector.x + XdistFromTarget;
it_vector.y = vector.y + YdistFromTarget;
it_vector.z = vector.z + ZdistFromTarget;*/
//Inverse of rotation about z-axis
izrot_vector.x = vector.x*ZRotMatrixInv[0][0] + vector.y*ZRotMatrixInv[0][1] + vector.z*ZRotMatrixInv[0][2];
izrot_vector.y = vector.x*ZRotMatrixInv[1][0] + vector.y*ZRotMatrixInv[1][1] + vector.z*ZRotMatrixInv[1][2];
izrot_vector.z = vector.x*ZRotMatrixInv[2][0] + vector.y*ZRotMatrixInv[2][1] + vector.z*ZRotMatrixInv[2][2];
//Inverse of rotation about x-axis
ixzrot_vector.x = izrot_vector.x*XRotMatrixInv[0][0] + izrot_vector.y*XRotMatrixInv[0][1] + izrot_vector.z*XRotMatrixInv[0][2];
ixzrot_vector.y = izrot_vector.x*XRotMatrixInv[1][0] + izrot_vector.y*XRotMatrixInv[1][1] + izrot_vector.z*XRotMatrixInv[1][2];
ixzrot_vector.z = izrot_vector.x*XRotMatrixInv[2][0] + izrot_vector.y*XRotMatrixInv[2][1] + izrot_vector.z*XRotMatrixInv[2][2];
return ixzrot_vector;
}
void SabreDetGeometry::Recenter(double x, double y) {
XdistFromTarget = x;
YdistFromTarget = y;
}
double SabreDetGeometry::GetScatteringAngle(int ringch, int wedgech) {
if(!CheckRingChannel(ringch) || !CheckWedgeChannel(wedgech)) return 0;
CartCoords these_coords = GetCoordinates(ringch, wedgech);
double r = sqrt(pow(these_coords.x, 2.) + pow(these_coords.y, 2.) + pow(these_coords.z, 2.));
return acos(these_coords.z/r);
}
/*
*/
bool SabreDetGeometry::IsInside(double theta, double phi) {
//Make unit vector
CartCoords unit, origin;
origin.x = 0; origin.y = 0; origin.z = ZdistFromTarget;
origin = InverseTransformVector(origin);
unit.x = sin(theta)*cos(phi); unit.y = sin(theta)*sin(phi); unit.z = cos(theta);
unit = InverseTransformVector(unit);
double s = origin.y/unit.y;
double px = origin.x + s*unit.x;
double py = origin.z + s*unit.z;
double radius = std::sqrt(px*px + py*py);
std::cout<<"radius: "<<radius<<std::endl;
double phi_in_detector = unit.GetPhi();
double theta_in_detector = unit.GetTheta();
double xy_radius = std::sqrt(unit.x*unit.x + unit.y*unit.y);
double xz_radius = std::sqrt(unit.x*unit.x + unit.z*unit.z);
double val = std::atan2(py, px);
if(val < 0 ) val += M_PI*2.0;
std::cout<<"val: "<<val*180.0/M_PI<<std::endl;
std::cout<<"phi in det: "<<phi_in_detector*180.0/M_PI<<std::endl;
std::cout<<"theta in det: "<<theta_in_detector*180.0/M_PI<<std::endl;
std::cout<<"xy radius: "<<xy_radius<<std::endl;
std::cout<<"xz radius: "<<xz_radius<<std::endl;
std::cout<<"unit -- x: "<<unit.x<<" y: "<<unit.y<<" z: "<<unit.z<<std::endl;
std::cout<<" phi range tilt -- min: "<<Get_Wedge_Tilted_Phi(0,1)*180.0/M_PI<<" max: "<<Get_Wedge_Tilted_Phi(7,0)*180.0/M_PI<<std::endl;
std::cout<<" phi range flat -- min: "<<Get_Wedge_Flat_Phi(0,1)*180.0/M_PI<<" max: "<<Get_Wedge_Flat_Phi(7,0)*180.0/M_PI<<std::endl;
std::cout<<" wedge 0 1 theta -- "<<Get_Wedge_Tilted_Theta(0,1)*180.0/M_PI<<std::endl;
std::cout<<" ring 0 0 x: "<<Get_Ring_Flat_X(0,0)<<std::endl;
std::cout<<" ring 0 1 x: "<<Get_Ring_Flat_X(0,1)<<std::endl;
std::cout<<" ring 0 0 y: "<<Get_Ring_Flat_Y(0,0)<<std::endl;
std::cout<<" ring 0 1 y: "<<Get_Ring_Flat_Y(0,1)<<std::endl;
std::cout<<" t ring 0 0 x: "<<Get_Ring_Tilted_X(0,0)<<std::endl;
std::cout<<" t ring 0 1 x: "<<Get_Ring_Tilted_X(0,1)<<std::endl;
std::cout<<" t ring 0 0 y: "<<Get_Ring_Tilted_Y(0,0)<<std::endl;
std::cout<<" t ring 0 1 y: "<<Get_Ring_Tilted_Y(0,1)<<std::endl;
std::cout<<"Test inversion"<<std::endl;
CartCoords test = InverseTransformVector(wedgech_tilted_cart[0][1]);
std::cout<<"Inverted x: "<<test.x<<" y: "<<test.y<<" z: "<<test.z<<std::endl;
std::cout<<"Given x: "<<wedgech_flat_cart[0][1].x<<" y: "<<wedgech_flat_cart[0][1].y<<" z: "<<wedgech_flat_cart[0][1].z<<std::endl;
//Check phi condition
//double phi_in_detector = phi - beamPhi_central;
double phi_min = Get_Wedge_Tilted_Phi(0,1);
double phi_max = Get_Wedge_Tilted_Phi(7,0);
bool passed;
if(phi > phi_min && phi < phi_max) {
passed = true;
} else {
return false;
}
std::cout<<" passed phi "<<std::endl;
if(xy_radius>Router_flat || xy_radius<Rinner_flat) passed = false;
return passed;
}
CartCoords SabreDetGeometry::GetTrajectoryCoordinates(double theta, double phi) {
CartCoords temp;
temp.x = 0; temp.y = 0; temp.z = 0;
int ringchan = -1, wedgechan = -1;
double phi_in_detector = phi - beamPhi_central;
double phi_min, phi_max;
for(int i=0; i<NUMWEDGES; i++) {
phi_min = wedgech_flat_cart[i][0].GetPhi();
phi_max = wedgech_flat_cart[i][1].GetPhi();
if(phi_in_detector > phi_min && phi_in_detector < phi_max) {
wedgechan = i;
break;
}
}
if(wedgechan == -1) return temp;
CartCoords innerPosition_flat, innerPosition_tilt;
CartCoords outerPosition_flat, outerPosition_tilt;
double r_inner, r_outer;
double theta_min, theta_max;
//Since we passed phi, use the phi coordinate to calculate the point at the "top" and
//bottom of the detector at the given phi, then transform to the lampshade frame.
for(int i=0; i<NUMRINGS; i++) {
r_inner = Rinner_flat + deltaR_flat_ring*(i);
r_outer = Rinner_flat + deltaR_flat_ring*(i+1);
innerPosition_flat.x = r_inner*sin(phi_in_detector);
innerPosition_flat.y = r_inner*cos(phi_in_detector);
innerPosition_flat.z = 0;
innerPosition_tilt = TransformVector(innerPosition_flat);
outerPosition_flat.x = r_outer*sin(phi_in_detector);
outerPosition_flat.y = r_outer*cos(phi_in_detector);
outerPosition_flat.z = 0;
outerPosition_tilt = TransformVector(outerPosition_flat);
//Get the theta values
theta_max = innerPosition_tilt.GetTheta();
theta_min = outerPosition_tilt.GetTheta();
if(theta<=theta_max && theta>=theta_min) {
ringchan = i;
break;
}
}
if(ringchan == -1) return temp;
return GetCoordinates(ringchan, wedgechan);
}
//Coordinate functions
double CartCoords::GetTheta() {
double r = std::sqrt(std::pow(x, 2.) + std::pow(y, 2.) + std::pow(z, 2.));
return std::acos(z/r);
}
double CartCoords::GetR() {
return std::sqrt(std::pow(x, 2.) + std::pow(y, 2.) + std::pow(z, 2.));
}
double CartCoords::GetPhi() {
/*double r = std::sqrt(std::pow(x, 2.) + std::pow(y, 2.));
if((x >= 0 && y >= 0) || (x <= 0 && y >= 0)) return std::acos(x/r);
else if((x <= 0 && y <= 0) || (x >= 0 && y <= 0)) return (2.0*M_PI - std::acos(x/r));
else return 0.0;*/
double phi = std::atan2(y, x);
if(phi<0) phi += M_PI*2.0;
return phi;
}
double CartCoords::GetDetectorPhi() {
/*double r = std::sqrt(std::pow(x, 2.) + std::pow(y, 2.));
if((x >= 0 && y >= 0) || (x <= 0 && y >= 0)) return std::acos(x/r);
else if((x <= 0 && y <= 0) || (x >= 0 && y <= 0)) return (2.0*M_PI - std::acos(x/r));
else return 0.0;*/
double phi = std::atan2(x, y);
if(phi<0) phi += M_PI*2.0;
return phi;
}

89
src/SabreEfficiency.cpp
View File

@ -17,6 +17,29 @@ SabreEfficiency::SabreEfficiency() :
detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI2*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI2*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
G3Vec coords;
for(int i=0; i<5; i++) {
for(int j=0; j<detectors[i].GetNumberOfRings(); j++) {
for(int k=0; k<4; k++) {
coords = detectors[i].GetRingTiltCoords(j, k);
ringxs.push_back(coords.GetX());
ringys.push_back(coords.GetY());
ringzs.push_back(coords.GetZ());
}
}
}
for(int i=0; i<5; i++) {
for(int j=0; j<detectors[i].GetNumberOfWedges(); j++) {
for(int k=0; k<4; k++) {
coords = detectors[i].GetWedgeTiltCoords(j, k);
wedgexs.push_back(coords.GetX());
wedgeys.push_back(coords.GetY());
wedgezs.push_back(coords.GetZ());
}
}
}
} }
SabreEfficiency::~SabreEfficiency() {} SabreEfficiency::~SabreEfficiency() {}
@ -60,15 +83,17 @@ void SabreEfficiency::RunDecay(const char* file) {
tree->SetBranchAddress("residual", &resid); tree->SetBranchAddress("residual", &resid);
double nevents = tree->GetEntries(); double nevents = tree->GetEntries();
std::vector<double> resid_thetas, eject_thetas; std::vector<double> resid_xs, eject_xs;
std::vector<double> resid_phis, eject_phis; std::vector<double> resid_ys, eject_ys;
std::vector<double> resid_kes, eject_kes; std::vector<double> resid_zs, eject_zs;
//Progress tracking //Progress tracking
int percent5 = nevents*0.05; int percent5 = nevents*0.05;
int count = 0; int count = 0;
int npercent = 0; int npercent = 0;
G3Vec coordinates;
for(int i=0; i<tree->GetEntries(); i++) { for(int i=0; i<tree->GetEntries(); i++) {
if(++count == percent5) {//Show progress every 5% if(++count == percent5) {//Show progress every 5%
npercent++; npercent++;
@ -80,10 +105,11 @@ void SabreEfficiency::RunDecay(const char* file) {
if(eject->KE >= ENERGY_THRESHOLD) { if(eject->KE >= ENERGY_THRESHOLD) {
for(auto& det : detectors) { for(auto& det : detectors) {
if(det.IsInside(eject->theta, eject->phi)) { coordinates = det.GetTrajectoryCoordinates(eject->theta, eject->phi);
eject_thetas.push_back(eject->theta); if(coordinates.GetX() != 0) {
eject_phis.push_back(eject->phi); eject_xs.push_back(coordinates.GetX());
eject_kes.push_back(eject->KE); eject_ys.push_back(coordinates.GetY());
eject_zs.push_back(coordinates.GetZ());
break; break;
} }
} }
@ -91,10 +117,10 @@ void SabreEfficiency::RunDecay(const char* file) {
if(resid->KE > ENERGY_THRESHOLD) { if(resid->KE > ENERGY_THRESHOLD) {
for(auto& det : detectors) { for(auto& det : detectors) {
if(det.IsInside(resid->theta, resid->phi)) { if(det.GetTrajectoryCoordinates(resid->theta, resid->phi).GetX() != 0) {
resid_thetas.push_back(resid->theta); resid_xs.push_back(coordinates.GetX());
resid_phis.push_back(resid->phi); resid_ys.push_back(coordinates.GetY());
resid_kes.push_back(resid->KE); resid_zs.push_back(coordinates.GetZ());
break; break;
} }
} }
@ -102,15 +128,42 @@ void SabreEfficiency::RunDecay(const char* file) {
} }
double ejecteff = ((double) eject_thetas.size())/nevents; double ejecteff = ((double) eject_xs.size())/nevents;
double resideff = ((double) resid_thetas.size())/nevents; double resideff = ((double) resid_xs.size())/nevents;
TParameter<double> eject_eff("Light Breakup Efficiency", ejecteff); TParameter<double> eject_eff("Light Breakup Efficiency", ejecteff);
TParameter<double> resid_eff("Heavy Breakup Efficiency", resideff); TParameter<double> resid_eff("Heavy Breakup Efficiency", resideff);
TGraph2D* gde = new TGraph2D(eject_xs.size(), &(eject_xs[0]), &(eject_ys[0]), &(eject_zs[0]));
gde->SetName("detected_eject_points");
gde->SetMarkerStyle(2);
gde->SetMarkerColor(2);
TGraph2D* gr = new TGraph2D(ringxs.size(), &(ringxs[0]), &(ringys[0]), &(ringzs[0]));
gr->SetName("ring_detector_edges");
gr->SetTitle("SABRE Detector; x(m); y(m); z(m)");
gr->SetMarkerStyle(2);
TGraph2D* gw = new TGraph2D(wedgexs.size(), &(wedgexs[0]), &(wedgeys[0]), &(wedgezs[0]));
gw->SetName("wedge_detector_edges");
gw->SetTitle("SABRE Detector Wedges; x(m); y(m); z(m)");
gw->SetMarkerStyle(2);
TCanvas* canvas = new TCanvas();
canvas->SetName("detectors_and_particles");
canvas->cd();
gr->Draw("AP");
gw->Draw("same P");
gde->Draw("same P");
input->cd(); input->cd();
eject_eff.Write(); eject_eff.Write();
resid_eff.Write(); resid_eff.Write();
gr->Write();
gw->Write();
gde->Write();
canvas->Write();
input->Close(); input->Close();
} }
void SabreEfficiency::Run2Step(const char* file) { void SabreEfficiency::Run2Step(const char* file) {
@ -146,7 +199,7 @@ void SabreEfficiency::Run2Step(const char* file) {
if(break1->KE >= ENERGY_THRESHOLD) { if(break1->KE >= ENERGY_THRESHOLD) {
for(auto& det : detectors) { for(auto& det : detectors) {
if(det.IsInside(break1->theta, break1->phi)) { if(det.GetTrajectoryCoordinates(break1->theta, break1->phi).GetX() != 0) {
b1_thetas.push_back(break1->theta); b1_thetas.push_back(break1->theta);
b1_phis.push_back(break1->phi); b1_phis.push_back(break1->phi);
b1_kes.push_back(break1->KE); b1_kes.push_back(break1->KE);
@ -157,7 +210,7 @@ void SabreEfficiency::Run2Step(const char* file) {
if(break2->KE > ENERGY_THRESHOLD) { if(break2->KE > ENERGY_THRESHOLD) {
for(auto& det : detectors) { for(auto& det : detectors) {
if(det.IsInside(break2->theta, break2->phi)) { if(det.GetTrajectoryCoordinates(break2->theta, break2->phi).GetX() != 0) {
b2_thetas.push_back(break2->theta); b2_thetas.push_back(break2->theta);
b2_phis.push_back(break2->phi); b2_phis.push_back(break2->phi);
b2_kes.push_back(break2->KE); b2_kes.push_back(break2->KE);
@ -214,7 +267,7 @@ void SabreEfficiency::Run3Step(const char* file) {
if(break1->KE > ENERGY_THRESHOLD) { if(break1->KE > ENERGY_THRESHOLD) {
for(auto& det : detectors) { for(auto& det : detectors) {
if(det.IsInside(break1->theta, break1->phi)) { if(det.GetTrajectoryCoordinates(break1->theta, break1->phi).GetX() != 0) {
b1_thetas.push_back(break1->theta); b1_thetas.push_back(break1->theta);
b1_phis.push_back(break1->phi); b1_phis.push_back(break1->phi);
b1_kes.push_back(break1->KE); b1_kes.push_back(break1->KE);
@ -226,7 +279,7 @@ void SabreEfficiency::Run3Step(const char* file) {
if(break3->KE > ENERGY_THRESHOLD) { if(break3->KE > ENERGY_THRESHOLD) {
for(auto& det : detectors) { for(auto& det : detectors) {
if(det.IsInside(break3->theta, break3->phi)) { if(det.GetTrajectoryCoordinates(break3->theta, break3->phi).GetX() != 0) {
b3_thetas.push_back(break3->theta); b3_thetas.push_back(break3->theta);
b3_phis.push_back(break3->phi); b3_phis.push_back(break3->phi);
b3_kes.push_back(break3->KE); b3_kes.push_back(break3->KE);
@ -237,7 +290,7 @@ void SabreEfficiency::Run3Step(const char* file) {
if(break4->KE > ENERGY_THRESHOLD) { if(break4->KE > ENERGY_THRESHOLD) {
for(auto& det : detectors) { for(auto& det : detectors) {
if(det.IsInside(break4->theta, break4->phi)) { if(det.GetTrajectoryCoordinates(break4->theta, break4->phi).GetX() != 0) {
b4_thetas.push_back(break4->theta); b4_thetas.push_back(break4->theta);
b4_phis.push_back(break4->phi); b4_phis.push_back(break4->phi);
b4_kes.push_back(break4->KE); b4_kes.push_back(break4->KE);

1
src/TwoStepSystem.cpp
View File

@ -87,4 +87,5 @@ void TwoStepSystem::RunSystem() {
step2.TurnOnResidualEloss(); step2.TurnOnResidualEloss();
step2.Calculate(); step2.Calculate();
} }

4
src/kinematics_dict.cxx
View File

@ -115,8 +115,8 @@ namespace {
0 0
}; };
static const char* includePaths[] = { static const char* includePaths[] = {
"/usr/local/root/include/", "/home/gordon/cern/root-6.22.02/root-install/include/",
"/Users/gordonmccann/Desktop/kinematics/", "/home/gordon/Kinematics/",
0 0
}; };
static const char* fwdDeclCode = R"DICTFWDDCLS( static const char* fwdDeclCode = R"DICTFWDDCLS(

49
src/main.cpp
View File

@ -1,21 +1,16 @@
#include <iostream> #include <iostream>
#include "Kinematics.h" #include "Kinematics.h"
#include "SabreEfficiency.h" #include "SabreEfficiency.h"
#include "SabreDetector.h"
#include "KinematicsExceptions.h" #include "KinematicsExceptions.h"
#include <TGraph2D.h>
#include <TGraph.h>
#include <TFile.h>
#include <TAxis.h>
#include <TCanvas.h>
int main(int argc, char** argv) { int main(int argc, char** argv) {
if(argc<2) { if(argc<2) {
std::cerr<<"Incorrect number of arguments!"<<std::endl; std::cerr<<"Incorrect number of arguments!"<<std::endl;
return 1; return 1;
} }
/*Kinematics calculator; Kinematics calculator;
try { try {
if(!calculator.LoadConfig(argv[1])) { if(!calculator.LoadConfig(argv[1])) {
return 1; return 1;
@ -28,19 +23,19 @@ int main(int argc, char** argv) {
} }
SabreEfficiency sabre; SabreEfficiency sabre;
sabre.SetReactionType(calculator.GetReactionType()); sabre.SetReactionType(calculator.GetReactionType());
sabre.CalculateEfficiency(calculator.GetOutputName());*/ sabre.CalculateEfficiency(calculator.GetOutputName());
std::vector<SabreDetGeometry> detectors; /*std::vector<SabreDetector> detectors;
const double INNER_R = 0.0326; const double INNER_R = 0.0326;
const double OUTER_R = 0.1351; const double OUTER_R = 0.1351;
const double TILT = 40.0; const double TILT = 40.0;
const double DIST_2_TARG = 0.1245; const double DIST_2_TARG = -0.1245;
const double PHI_COVERAGE = 54.4; //delta phi for each det const double PHI_COVERAGE = 54.4; //delta phi for each det
const double PHI0 = 36.0; //center phi values for each det in array const double PHI0 = 234.0; //center phi values for each det in array
const double PHI1 = 108.0; //# is equal to detID in channel map const double PHI1 = 162.0; //# is equal to detID in channel map
const double PHI2 = 324.0; const double PHI2 = 306.0;
const double PHI3 = 252.0; const double PHI3 = 18.0;
const double PHI4 = 180.0; const double PHI4 = 90.0;
const double DEG2RAD = M_PI/180.0; const double DEG2RAD = M_PI/180.0;
detectors.reserve(5); detectors.reserve(5);
@ -50,15 +45,23 @@ int main(int argc, char** argv) {
detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG);
double theta = 145*M_PI/180.0; double phi = 92.1428*M_PI/180.0; double theta, phi, expected_flat_t, expected_flat_p;
std::cout<<"theta: "<<theta<<" phi: "<<phi<<std::endl; for(int h=0; h<5; h++) {
for(int i=0; i<5; i++) { for(int j=0; j<16; j++) {
if(detectors[i].IsInside(theta, phi)) { for(int k=0; k<4; k ++) {
std::cout<<"Caught it in detector: "<<i<<std::endl; theta = detectors[h].GetRingTiltCoords(j, k).GetTheta();
phi = detectors[h].GetRingTiltCoords(j, k).GetPhi();
expected_flat_p = detectors[h].GetRingFlatCoords(j, k).GetPhi();
for(int i=0; i<5; i++) {
if(detectors[i].GetTrajectoryCoordinates(theta, phi).GetX() != 0) {
break;
} else if(i == 4) {
std::cout<<" Not found! detector: "<<h<<" ring: "<<j<<" corner: "<<k<<" theta: "<<theta/DEG2RAD<<" phi: "<<phi/DEG2RAD<<" flat_p: "<<expected_flat_p/DEG2RAD<<std::endl;
}
}
}
} }
break; }*/
}
return 0; return 0;

396
testplots/testplots.cpp
View File

@ -1,396 +0,0 @@
#include <iostream>
#include <cmath>
#include <fstream>
#include <TGraph.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TGraph2D.h>
#include <TLatex.h>
#include <TRandom3.h>
#include "SabreDetector.h"
using namespace std;
int main() {
const int NUMDETS = 5;
SabreDetGeometry **sabdet = new SabreDetGeometry*[NUMDETS];
//dimensions common to all dets:
const double INNER_RADIUS = 0.0326; //meters
const double OUTER_RADIUS = 0.1351; //meters
const double TILT_FROM_VERTICAL = 40.0; //degrees
const double DIST_FROM_TARGET = 0.1245; //meters
sabdet[0] = new SabreDetGeometry(INNER_RADIUS,
OUTER_RADIUS,
54.4*M_PI/180., //delta phi
36.*M_PI/180., //center phi
TILT_FROM_VERTICAL*M_PI/180.,
DIST_FROM_TARGET);
sabdet[1] = new SabreDetGeometry(INNER_RADIUS,
OUTER_RADIUS,
54.4*M_PI/180.,
108.*M_PI/180.,
TILT_FROM_VERTICAL*M_PI/180.,
DIST_FROM_TARGET);
sabdet[2] = new SabreDetGeometry(INNER_RADIUS,
OUTER_RADIUS,
54.4*M_PI/180.,
180.*M_PI/180.,
TILT_FROM_VERTICAL*M_PI/180.,
DIST_FROM_TARGET);
sabdet[3] = new SabreDetGeometry(INNER_RADIUS,
OUTER_RADIUS,
54.4*M_PI/180.,
252.*M_PI/180.,
TILT_FROM_VERTICAL*M_PI/180.,
DIST_FROM_TARGET);
sabdet[4] = new SabreDetGeometry(INNER_RADIUS,
OUTER_RADIUS,
54.4*M_PI/180.,
324.*M_PI/180.,
TILT_FROM_VERTICAL*M_PI/180.,
DIST_FROM_TARGET);
const int NUMCORNERS = 4;
const int NUMRINGS = sabdet[0]->NumRings();
const int NUMWEDGES = sabdet[0]->NumWedges();
const int NUMRINGPTS = NUMDETS*NUMRINGS*NUMCORNERS;
const int NUMWEDGEPTS = NUMDETS*NUMWEDGES*NUMCORNERS;
//const int NUMMCCOUNTS = 1E8;
double ring_xpts_flat[NUMRINGPTS];
double ring_ypts_flat[NUMRINGPTS];
double ring_zpts_flat[NUMRINGPTS];
double ring_rpts_flat[NUMRINGPTS];
double ring_tpts_flat[NUMRINGPTS];
double ring_ppts_flat[NUMRINGPTS];
double wedge_xpts_flat[NUMWEDGEPTS];
double wedge_ypts_flat[NUMWEDGEPTS];
double wedge_zpts_flat[NUMWEDGEPTS];
double wedge_rpts_flat[NUMWEDGEPTS];
double wedge_tpts_flat[NUMWEDGEPTS];
double wedge_ppts_flat[NUMWEDGEPTS];
double ring_xpts_tilted[NUMRINGPTS];
double ring_ypts_tilted[NUMRINGPTS];
double ring_zpts_tilted[NUMRINGPTS];
double ring_rpts_tilted[NUMRINGPTS];
double ring_tpts_tilted[NUMRINGPTS];
double ring_ppts_tilted[NUMRINGPTS];
double wedge_xpts_tilted[NUMWEDGEPTS];
double wedge_ypts_tilted[NUMWEDGEPTS];
double wedge_zpts_tilted[NUMWEDGEPTS];
double wedge_rpts_tilted[NUMWEDGEPTS];
double wedge_tpts_tilted[NUMWEDGEPTS];
double wedge_ppts_tilted[NUMWEDGEPTS];
ofstream checkfile;
checkfile.open("checkfile.txt");
checkfile << "RINGS\n"
<< "det ch corner x_fl y_fl z_fl r_fl t_fl p_fl x_ti y_ti z_ti r_ti t_ti p_ti\n";
for (int d=0; d<NUMDETS; d++) {
for (int i=0; i<NUMRINGS; i++) {
for (int j=0; j<NUMCORNERS; j++) {
int index = NUMCORNERS*NUMRINGS*d + NUMCORNERS*i + j;
ring_xpts_flat[index] = sabdet[d]->Get_Ring_Flat_X(i,j)*100;
ring_ypts_flat[index] = sabdet[d]->Get_Ring_Flat_Y(i,j)*100;
ring_zpts_flat[index] = sabdet[d]->Get_Ring_Flat_Z(i,j)*100;
ring_rpts_flat[index] = sabdet[d]->Get_Ring_Flat_R(i,j)*100;
ring_tpts_flat[index] = sabdet[d]->Get_Ring_Flat_Theta(i,j)*180./M_PI;
ring_ppts_flat[index] = sabdet[d]->Get_Ring_Flat_Phi(i,j)*180./M_PI;
ring_xpts_tilted[index] = sabdet[d]->Get_Ring_Tilted_X(i,j)*100;
ring_ypts_tilted[index] = sabdet[d]->Get_Ring_Tilted_Y(i,j)*100;
ring_zpts_tilted[index] = sabdet[d]->Get_Ring_Tilted_Z(i,j)*100;
ring_rpts_tilted[index] = sabdet[d]->Get_Ring_Tilted_R(i,j)*100;
ring_tpts_tilted[index] = sabdet[d]->Get_Ring_Tilted_Theta(i,j)*180./M_PI;
ring_ppts_tilted[index] = sabdet[d]->Get_Ring_Tilted_Phi(i,j)*180./M_PI;
checkfile << d << " " << i << " " << j << " "
<< ring_xpts_flat[index] << " "
<< ring_ypts_flat[index] << " "
<< ring_zpts_flat[index] << " "
<< ring_rpts_flat[index] << " "
<< ring_tpts_flat[index] << " "
<< ring_ppts_flat[index] << " "
<< ring_xpts_tilted[index] << " "
<< ring_ypts_tilted[index] << " "
<< ring_zpts_tilted[index] << " "
<< ring_rpts_tilted[index] << " "
<< ring_tpts_tilted[index] << " "
<< ring_ppts_tilted[index] << endl;
}
}
}
checkfile << "WEDGES\n"
<< "det ch corner x_fl y_fl z_fl r_fl t_fl p_fl x_ti y_ti z_ti r_ti t_ti p_ti\n";
for (int d=0; d<NUMDETS; d++) {
for (int i=0; i<NUMWEDGES; i++) {
for (int j=0; j<NUMCORNERS; j++) {
int index = NUMCORNERS*NUMWEDGES*d + NUMCORNERS*i + j;
wedge_xpts_flat[index] = sabdet[d]->Get_Wedge_Flat_X(i,j)*100;
wedge_ypts_flat[index] = sabdet[d]->Get_Wedge_Flat_Y(i,j)*100;
wedge_zpts_flat[index] = sabdet[d]->Get_Wedge_Flat_Z(i,j)*100;
wedge_rpts_flat[index] = sabdet[d]->Get_Wedge_Flat_R(i,j)*100;
wedge_tpts_flat[index] = sabdet[d]->Get_Wedge_Flat_Theta(i,j)*180./M_PI;
wedge_ppts_flat[index] = sabdet[d]->Get_Wedge_Flat_Phi(i,j)*180./M_PI;
wedge_xpts_tilted[index] = sabdet[d]->Get_Wedge_Tilted_X(i,j)*100;
wedge_ypts_tilted[index] = sabdet[d]->Get_Wedge_Tilted_Y(i,j)*100;
wedge_zpts_tilted[index] = sabdet[d]->Get_Wedge_Tilted_Z(i,j)*100;
wedge_rpts_tilted[index] = sabdet[d]->Get_Wedge_Tilted_R(i,j)*100;
wedge_tpts_tilted[index] = sabdet[d]->Get_Wedge_Tilted_Theta(i,j)*180./M_PI;
wedge_ppts_tilted[index] = sabdet[d]->Get_Wedge_Tilted_Phi(i,j)*180./M_PI;
checkfile << d << " " << i << " " << j << " "
<< wedge_xpts_flat[index] << " "
<< wedge_ypts_flat[index] << " "
<< wedge_zpts_flat[index] << " "
<< wedge_rpts_flat[index] << " "
<< wedge_tpts_flat[index] << " "
<< wedge_ppts_flat[index] << " "
<< wedge_xpts_tilted[index] << " "
<< wedge_ypts_tilted[index] << " "
<< wedge_zpts_tilted[index] << " "
<< wedge_rpts_tilted[index] << " "
<< wedge_tpts_tilted[index] << " "
<< wedge_ppts_tilted[index] << endl;
}
}
}
checkfile.close();
/*
UNDER CONSTRUCTION: MC efficiency integration
TRandom3 *rand = new TRandom3();
rand->SetSeed();
int hits=0;
bool didhit=false;
double randx=0;
double randy=0;
double randz=0;
double randr=0;
double randt=0;
double randp=0;
for (int mc=0; mc<NUMMCCOUNTS; mc++) {
didhit = false;
randx = rand->Uniform(0,1);
randy = rand->Uniform(0,1);
randz = rand->Uniform(0,1);
randr = sqrt(pow(randx,2)+pow(randy,2)+pow(randz,2));
randt = acos(randz/randr)*180./M_PI;
randp = atan(randy/randx)*180./M_PI;
if (randx <= 0 && randy >= 0) randp += 90;
else if (randx <= 0 && randy <= 0) randp += 90;
else if (randx >= 0 && randy <= 0) randp += 270;
else if (randx >= 0 && randy >= 0) randp += 270;
for (int d=0; d<NUMDETS; d++) {
if (didhit) break;
for (int f=0; f<NUMRINGS; f++) {
if (didhit) break;
for (int b=0; b<NUMWEDGES; b++) {
if (randt <= sabdet[d]->Get_Ring_Tilted_Theta(f,)) {
hit++;
didhit=true;
break;
}
}
}
}
}
cout << "Geom eff = " << 1.*hits/NUMMCCOUNTS*100 << "%" << endl;
*/
TGraph ring_flat_xy_gr(NUMRINGPTS,ring_xpts_flat,ring_ypts_flat);
ring_flat_xy_gr.SetMarkerStyle(2);
ring_flat_xy_gr.SetName("ring_flat_xy_gr");
TGraph wedge_flat_xy_gr(NUMWEDGEPTS,wedge_xpts_flat,wedge_ypts_flat);
wedge_flat_xy_gr.SetMarkerStyle(2);
wedge_flat_xy_gr.SetName("wedge_flat_xy_gr");
wedge_flat_xy_gr.SetTitle("Flat (Untilted);x (cm);y (cm)");
TGraph ring_flat_zy_gr(NUMRINGPTS,ring_zpts_flat,ring_ypts_flat);
ring_flat_zy_gr.SetMarkerStyle(2);
ring_flat_zy_gr.SetName("ring_flat_zy_gr");
TGraph wedge_flat_zy_gr(NUMWEDGEPTS,wedge_zpts_flat,wedge_ypts_flat);
wedge_flat_zy_gr.SetMarkerStyle(2);
wedge_flat_zy_gr.SetName("wedge_flat_zy_gr");
wedge_flat_zy_gr.SetTitle("Flat (Untilted);z (cm);y (cm)");
TGraph ring_flat_pt_gr(NUMRINGPTS,ring_ppts_flat,ring_tpts_flat);
ring_flat_pt_gr.SetMarkerStyle(2);
ring_flat_pt_gr.SetName("ring_flat_pt_gr");
TGraph wedge_flat_pt_gr(NUMWEDGEPTS,wedge_ppts_flat,wedge_tpts_flat);
wedge_flat_pt_gr.SetMarkerStyle(2);
wedge_flat_pt_gr.SetName("wedge_flat_pt_gr");
wedge_flat_pt_gr.SetTitle("Flat (Untilted);#phi (deg.);#theta (deg.)");
TGraph2D ring_flat_2d_gr(NUMRINGPTS,ring_xpts_flat,ring_zpts_flat,ring_ypts_flat);
ring_flat_2d_gr.SetMarkerStyle(2);
ring_flat_2d_gr.SetName("ring_flat_2d_gr");
TGraph2D wedge_flat_2d_gr(NUMWEDGEPTS,wedge_xpts_flat,wedge_zpts_flat,wedge_ypts_flat);
wedge_flat_2d_gr.SetMarkerStyle(2);
wedge_flat_2d_gr.SetName("wedge_flat_2d_gr");
wedge_flat_2d_gr.SetTitle("Flat (Untilted);x (cm); z (cm); y (cm)");
TGraph ring_tilted_xy_gr(NUMRINGPTS,ring_xpts_tilted,ring_ypts_tilted);
ring_tilted_xy_gr.SetMarkerStyle(2);
ring_tilted_xy_gr.SetName("ring_tilted_xy_gr");
TGraph wedge_tilted_xy_gr(NUMWEDGEPTS,wedge_xpts_tilted,wedge_ypts_tilted);
wedge_tilted_xy_gr.SetMarkerStyle(2);
wedge_tilted_xy_gr.SetName("wedge_tilted_xy_gr");
wedge_tilted_xy_gr.SetTitle("Tilted;x (cm);y (cm)");
TGraph ring_tilted_zy_gr(NUMRINGPTS,ring_zpts_tilted,ring_ypts_tilted);
ring_tilted_zy_gr.SetMarkerStyle(2);
ring_tilted_zy_gr.SetName("ring_tilted_zy_gr");
TGraph wedge_tilted_zy_gr(NUMWEDGEPTS,wedge_zpts_tilted,wedge_ypts_tilted);
wedge_tilted_zy_gr.SetMarkerStyle(2);
wedge_tilted_zy_gr.SetName("wedge_tilted_zy_gr");
wedge_tilted_zy_gr.SetTitle("Tilted;z (cm);y (cm)");
TGraph ring_tilted_pt_gr(NUMRINGPTS,ring_ppts_tilted,ring_tpts_tilted);
ring_tilted_pt_gr.SetMarkerStyle(2);
ring_tilted_pt_gr.SetName("ring_tilted_pt_gr");
TGraph wedge_tilted_pt_gr(NUMWEDGEPTS,wedge_ppts_tilted,wedge_tpts_tilted);
wedge_tilted_pt_gr.SetMarkerStyle(2);
wedge_tilted_pt_gr.SetName("wedge_tilted_pt_gr");
wedge_tilted_pt_gr.SetTitle("Tilted;#phi (deg.);#theta (deg.)");
TGraph2D ring_tilted_2d_gr(NUMRINGPTS,ring_xpts_tilted,ring_zpts_tilted,ring_ypts_tilted);
ring_tilted_2d_gr.SetMarkerStyle(2);
ring_tilted_2d_gr.SetName("ring_tilted_2d_gr");
TGraph2D wedge_tilted_2d_gr(NUMWEDGEPTS,wedge_xpts_tilted,wedge_zpts_tilted,wedge_ypts_tilted);
wedge_tilted_2d_gr.SetMarkerStyle(2);
wedge_tilted_2d_gr.SetName("wedge_tilted_2d_gr");
wedge_tilted_2d_gr.SetTitle("Tilted;x (cm); z (cm); y (cm)");
int nCenters = 5*16*8;
int nCent_det = 16*8;
double centerx[nCenters], centery[nCenters], centerz[nCenters];
double **det_centx, **det_centy, **det_centz;
det_centx = new double*[5]; det_centy = new double*[5]; det_centz = new double*[5];
for(int i=0; i<5; i++) {
det_centx[i] = new double[nCent_det];
det_centy[i] = new double[nCent_det];
det_centz[i] = new double[nCent_det];
}
CartCoords point;
int index = 0;
int small_index;
for(int d=0; d<5; d++) {
small_index = 0;
for(int r=0; r<16; r++) {
for(int w=0; w<8; w++) {
point = sabdet[d]->GetCoordinates(r,w);
centerx[index] = point.x*100;
centery[index] = point.y*100;
centerz[index] = point.z*100;
det_centx[d][small_index] = point.x*100;
det_centy[d][small_index] = point.y*100;
det_centz[d][small_index] = point.z*100;
index++; small_index++;
}
}
}
TGraph2D centerpoint(nCenters, &centerx[0], &centerz[0], &centery[0]);
centerpoint.SetName("centerpoint");
centerpoint.SetTitle("centerpoint;x;z;y");
centerpoint.SetMarkerStyle(2);
centerpoint.SetMarkerColor(2);
TGraph2D centerpoint_d1(nCent_det, &(det_centx[0][0]), &(det_centz[0][0]), &(det_centy[0][0]));
centerpoint_d1.SetName("centerpoint_d1");
centerpoint_d1.SetMarkerStyle(2);
centerpoint_d1.SetMarkerColor(1);
TGraph2D centerpoint_d2(nCent_det, &(det_centx[1][0]), &(det_centz[1][0]), &(det_centy[1][0]));
centerpoint_d2.SetName("centerpoint_d2");
centerpoint_d2.SetMarkerStyle(2);
centerpoint_d2.SetMarkerColor(2);
TGraph2D centerpoint_d3(nCent_det, &(det_centx[2][0]), &(det_centz[2][0]), &(det_centy[2][0]));
centerpoint_d3.SetName("centerpoint_d3");
centerpoint_d3.SetMarkerStyle(2);
centerpoint_d3.SetMarkerColor(3);
TGraph2D centerpoint_d4(nCent_det, &(det_centx[3][0]), &(det_centz[3][0]), &(det_centy[3][0]));
centerpoint_d4.SetName("centerpoint_d4");
centerpoint_d4.SetMarkerStyle(2);
centerpoint_d4.SetMarkerColor(4);
TGraph2D centerpoint_d5(nCent_det, &(det_centx[4][0]), &(det_centz[4][0]), &(det_centy[4][0]));
centerpoint_d5.SetName("centerpoint_d5");
centerpoint_d5.SetMarkerStyle(2);
centerpoint_d5.SetMarkerColor(5);
TCanvas flat_xy_canv("flat_xy_canv");
wedge_flat_xy_gr.Draw("AP");
ring_flat_xy_gr.Draw("same P");
TCanvas flat_zy_canv("flat_zy_canv");
wedge_flat_zy_gr.Draw("AP");
ring_flat_zy_gr.Draw("same P");
TCanvas flat_pt_canv("flat_pt_canv");
wedge_flat_pt_gr.Draw("AP");
ring_flat_pt_gr.Draw("same P");
TCanvas flat_2d_canv("flat_2d_canv");
wedge_flat_2d_gr.Draw("P");
ring_flat_2d_gr.Draw("same P");
TCanvas tilted_xy_canv("tilted_xy_canv");
wedge_tilted_xy_gr.Draw("AP");
ring_tilted_xy_gr.Draw("same P");
TCanvas tilted_zy_canv("tilted_zy_canv");
wedge_tilted_zy_gr.Draw("AP");
ring_tilted_zy_gr.Draw("same P");
TCanvas tilted_pt_canv("tilted_pt_canv");
wedge_tilted_pt_gr.Draw("AP");
ring_tilted_pt_gr.Draw("same P");
TCanvas tilted_2d_canv("tilted_2d_canv");
wedge_tilted_2d_gr.Draw("P");
ring_tilted_2d_gr.Draw("same P");
centerpoint.Draw("same P");
TCanvas tilted_point("tilted_point");
centerpoint.Draw("P");
TCanvas tilted_point_det("tilted_point_det");
wedge_tilted_2d_gr.Draw("P");
ring_tilted_2d_gr.Draw("same P");
centerpoint_d1.Draw("P same");
centerpoint_d2.Draw("P same");
centerpoint_d3.Draw("P same");
centerpoint_d4.Draw("P same");
centerpoint_d5.Draw("P same");
TFile outfile("outfile.root","RECREATE");
ring_flat_xy_gr.Write();
wedge_flat_xy_gr.Write();
flat_xy_canv.Write();
ring_tilted_xy_gr.Write();
wedge_tilted_xy_gr.Write();
tilted_xy_canv.Write();
ring_flat_zy_gr.Write();
wedge_flat_zy_gr.Write();
flat_zy_canv.Write();
ring_tilted_zy_gr.Write();
wedge_tilted_zy_gr.Write();
tilted_zy_canv.Write();
ring_flat_pt_gr.Write();
wedge_flat_pt_gr.Write();
flat_pt_canv.Write();
ring_tilted_pt_gr.Write();
wedge_tilted_pt_gr.Write();
tilted_pt_canv.Write();
ring_flat_2d_gr.Write();
wedge_flat_2d_gr.Write();
flat_2d_canv.Write();
ring_tilted_2d_gr.Write();
wedge_tilted_2d_gr.Write();
tilted_2d_canv.Write();
tilted_point.Write();
tilted_point_det.Write();
outfile.Close();
delete [] sabdet;
return 0;
}

BIN
tp
View File

Binary file not shown.