From 66900ea9b4bba282504ec3356a9894ebf6ef0a59 Mon Sep 17 00:00:00 2001 From: gwm17 Date: Thu, 18 Aug 2022 10:31:16 -0400 Subject: [PATCH] Update gitignore --- .gitignore | 5 +- bin/.gitignore | 5 - bin/PyPlotViewer | 5 - bin/PyPlotter | 6 - include/Eloss_Tables.h | 133 ----------- include/EnergyLoss.h | 58 ----- include/MaskFile.h | 77 ------ include/Rotation.h | 78 ------- include/Vec3.h | 61 ----- include/Vec4.h | 69 ------ src/Mask/EnergyLoss.cpp | 243 ------------------- src/Mask/MaskFile.cpp | 364 ----------------------------- src/Mask/Nucleus.cpp | 47 ---- src/Mask/Rotation.cpp | 71 ------ src/Mask/Vec3.cpp | 48 ---- src/Mask/Vec4.cpp | 73 ------ src/Plotters/Python/.gitignore | 3 - src/Plotters/Python/MaskFile.py | 80 ------- src/Plotters/Python/NucData.py | 70 ------ src/Plotters/Python/Nucleus.py | 106 --------- src/Plotters/Python/PyPlotter.py | 135 ----------- src/Plotters/Python/ViewPyPlots.py | 32 --- 22 files changed, 3 insertions(+), 1766 deletions(-) delete mode 100644 bin/.gitignore delete mode 100755 bin/PyPlotViewer delete mode 100755 bin/PyPlotter delete mode 100644 include/Eloss_Tables.h delete mode 100644 include/EnergyLoss.h delete mode 100644 include/MaskFile.h delete mode 100644 include/Rotation.h delete mode 100644 include/Vec3.h delete mode 100644 include/Vec4.h delete mode 100644 src/Mask/EnergyLoss.cpp delete mode 100644 src/Mask/MaskFile.cpp delete mode 100644 src/Mask/Nucleus.cpp delete mode 100644 src/Mask/Rotation.cpp delete mode 100644 src/Mask/Vec3.cpp delete mode 100644 src/Mask/Vec4.cpp delete mode 100644 src/Plotters/Python/.gitignore delete mode 100755 src/Plotters/Python/MaskFile.py delete mode 100755 src/Plotters/Python/NucData.py delete mode 100755 src/Plotters/Python/Nucleus.py delete mode 100755 src/Plotters/Python/PyPlotter.py delete mode 100755 src/Plotters/Python/ViewPyPlots.py diff --git a/.gitignore b/.gitignore index 4b337c1..f0eaf4d 100644 --- a/.gitignore +++ b/.gitignore @@ -9,14 +9,15 @@ *.pcm *.o *.a +*.so *.swp -*.pickle .DS_Store Makefile *.make -__pycache__ .vscode/ build/ +lib/ +bin/ ###Keep this file### diff --git a/bin/.gitignore b/bin/.gitignore deleted file mode 100644 index 888fdbf..0000000 --- a/bin/.gitignore +++ /dev/null @@ -1,5 +0,0 @@ -###keep only the directory, not the contents### -* -!.gitignore -!PyPlotter -!PyPlotViewer \ No newline at end of file diff --git a/bin/PyPlotViewer b/bin/PyPlotViewer deleted file mode 100755 index c612b15..0000000 --- a/bin/PyPlotViewer +++ /dev/null @@ -1,5 +0,0 @@ -#!/bin/bash - -FileName=$1; - -./src/Plotters/Python/ViewPyPlots.py ${FileName} \ No newline at end of file diff --git a/bin/PyPlotter b/bin/PyPlotter deleted file mode 100755 index 253ddd3..0000000 --- a/bin/PyPlotter +++ /dev/null @@ -1,6 +0,0 @@ -#!/bin/bash - -InputFile=$1; -OutputFile=$2; - -./src/Plotters/Python/PyPlotter.py ${InputFile} ${OutputFile} \ No newline at end of file diff --git a/include/Eloss_Tables.h b/include/Eloss_Tables.h deleted file mode 100644 index 67b5abe..0000000 --- a/include/Eloss_Tables.h +++ /dev/null @@ -1,133 +0,0 @@ -#ifndef ELOSS_TABLES_H -#define ELOSS_TABLES_H - -namespace Mask { - - #define MAX_Z 93 //Maximum number of elements for which we have hydrogen coefficients - - /*Atomic Masses for elements H through U. Taken from ELAST data*/ - static double NATURAL_MASS[MAX_Z] = - {0, - 1.00797, 4.0026, 6.939, 9.0122, 10.818, - 12.01115, 14.0067, 15.9994, 18.99984, 20.183, - 22.9898, 24.312, 26.9815, 28.086, 30.9738, - 32.064, 35.453, 39.948, 39.102, 40.08, - 44.956, 47.90, 50.942, 51.996, 54.938, - 55.847, 58.933, 58.71, 63.54, 65.37, - 69.72, 72.59, 74.922, 78.96, 79.909, - 83.80, 85.47, 87.62, 88.909, 91.22, - 92.906, 95.94, 98., 101.07, 102.905, - 106.4, 107.87, 112.4, 114.82, 118.69, - 121.75, 127.60, 126.904, 131.3, 132.905, - 137.34, 138.91, 140.12, 140.907, 144.24, - 146., 150.35, 151.96, 157.25, 158.924, - 162.50, 164.93, 167.26, 168.934, 173.04, - 174.97, 178.49, 180.948, 183.85, 186.2, - 190.2, 192.2, 195.09, 196.967, 200.59, - 204.37, 207.19, 208.98, 209., 210., - 222., 223., 226., 227., 232.038, - 231., 238.03 }; - - /*Stopping power coefficients for hydrogen into Z - *Taken from Ziegler, Hydrogen Stopping Powers and Ranges in All Elements, Volume 3 of - *The Stopping and Ranges of Ions in Matter, 1977, Pergamon Press Inc. - *Expanded table for method given in O.G. Spanc - *Note: some elements represent extrapolations when there was no data to fit - */ - static double HYDROGEN_COEFF[MAX_Z][12] = { - /*Blank*/{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}, - /*H*/{1.262,1.44,242.6,1.2E4,0.1159,0.0005099,5.436E4,-5.052,2.049,-0.3044,0.01966,-0.0004659}, - /*He*/{1.229,1.397,484.5,5873,0.05225,0.00102,2.451E4,-2.158,0.8278,-0.1172,0.007259,-0.000166}, - /*Li*/{1.411,1.6,725.6,3013,0.04578,0.00153,2.147E4,-0.5831,0.562,-0.1183,0.009298,-0.0002498}, - /*Be*/{2.248,2.59,966,153.8,0.03475,0.002039,1.63E4,0.2779,0.1745,-0.05684,0.005155,-0.0001488}, - /*B*/{2.474,2.815,1206,1060,0.02855,0.002549,1.345E4,-2.445,1.283,-0.2205,0.0156,-0.000393}, - /*C*/{2.631,2.989,1445,957.2,0.02819,0.003059,1.322E4,-4.38,2.044,-0.3283,0.02221,-0.0005417}, - /*N*/{2.954,3.35,1683,1900,0.02513,0.003569,1.179E4,-5.054,2.325,-0.3713,0.02506,-0.0006109}, - /*O*/{2.652,3,1920,2000,0.0223,0.004079,1.046E4,-6.734,3.019,-0.4748,0.03171,-0.0007669}, - /*F*/{2.085,2.352,2157,2634,0.01816,0.004589,8517,-5.571,2.449,-0.3781,0.02483,-0.0005919}, - /*Ne*/{1.951,2.199,2393,2699,0.01568,0.005099,7353,-4.408,1.879,-0.2814,0.01796,-0.0004168}, - /*Na*/{2.542,2.869,2628,1854,0.01472,0.005609,6905,-4.959,2.073,-0.3054,0.01921,-0.0004403}, - /*Mg*/{3.792,4.293,2862,1009,0.01397,0.006118,6551,-5.51,2.266,-0.3295,0.02047,-0.0004637}, - /*Al*/{4.154,4.739,2766,164.5,0.02023,0.006628,6309,-6.061,2.46,-0.3535,0.02173,-0.0004871}, - /*Si*/{4.15,4.7,3329,550,0.01321,0.007138,6194,-6.294,2.538,-0.3628,0.0222,-0.0004956}, - /*P*/{3.232,3.647,3561,1560,0.01267,0.007648,5942,-6.527,2.616,-0.3721,0.02267,-0.000504}, - /*S*/{3.447,3.891,3792,1219,0.01211,0.008158,5678,-6.761,2.694,-0.3814,0.02314,-0.0005125}, - /*Cl*/{5.047,5.714,4023,878.6,0.01178,0.008668,5524,-6.994,2.773,-0.3907,0.02361,-0.0005209}, - /*Ar*/{5.731,6.5,4253,530,0.01123,0.009178,5268,-7.227,2.851,-0.4,0.02407,-0.0005294}, - /*K*/{5.151,5.833,4482,545.7,0.01129,0.009687,5295,-7.44,2.923,-0.4094,0.02462,-0.0005411}, - /*Ca*/{5.521,6.252,4710,553.3,0.01112,0.0102,5214,-7.653,2.995,-0.4187,0.02516,-0.0005529}, - /*Sc*/{5.201,5.884,4938,560.9,0.009995,0.01071,4688,-8.012,3.123,-0.435,0.02605,-0.0005707}, - /*Ti*/{4.862,5.496,5165,568.5,0.009474,0.01122,4443,-8.371,3.251,-0.4513,0.02694,-0.0005886}, - /*V*/{4.48,5.055,5391,952.3,0.009117,0.01173,4276,-8.731,3.379,-0.4676,0.02783,-0.0006064}, - /*Cr*/{3.983,4.489,5616,1336,0.008413,0.01224,3946,-9.09,3.507,-0.4838,0.02872,-0.0006243}, - /*Mn*/{3.469,3.907,5725,1461,0.008829,0.01275,3785,-9.449,3.635,-0.5001,0.02961,-0.0006421}, - /*Fe*/{3.519,3.963,6065,1243,0.007782,0.01326,3650,-9.809,3.763,-0.5164,0.0305,-0.00066}, - /*Co*/{3.14,3.535,6288,1372,0.007361,0.01377,3453,-10.17,3.891,-0.5327,0.03139,-0.0006779}, - /*Ni*/{3.553,4.004,6205,555.1,0.008763,0.01428,3297,-10.53,4.019,-0.549,0.03229,-0.0006957}, - /*Cu*/{3.696,4.175,4673,387.8,0.02188,0.01479,3174,-11.18,4.252,-0.5791,0.03399,-0.0007314}, - /*Zn*/{4.21,4.75,6953,295.2,0.006809,0.0153,3194,-11.57,4.394,-0.598,0.03506,-0.0007537}, - /*Ga*/{5.041,5.697,7173,202.6,0.006725,0.01581,3154,-11.95,4.537,-0.6169,0.03613,-0.0007759}, - /*Ge*/{5.554,6.3,6496,110,0.009689,0.01632,3097,-12.34,4.68,-0.6358,0.03721,-0.0007981}, - /*As*/{5.323,6.012,7611,292.5,0.006447,0.01683,3024,-12.72,4.823,-0.6547,0.03828,-0.0008203}, - /*Se*/{5.874,6.656,7395,117.5,0.007684,0.01734,3006,-13.11,4.965,-0.6735,0.03935,-0.0008425}, - /*Br*/{5.611,6.335,8046,365.2,0.006244,0.01785,2928,-13.4,5.083,-0.6906,0.04042,-0.0008675}, - /*Kr*/{6.411,7.25,8262,220,0.006087,0.01836,2855,-13.69,5.2,-0.7076,0.0415,-0.0008925}, - /*Rb*/{5.694,6.429,8478,292.9,0.006087,0.01886,2855,-13.92,5.266,-0.714,0.04173,-0.0008943}, - /*Sr*/{6.339,7.159,8693,330.3,0.006003,0.01937,2815,-14.14,5.331,-0.7205,0.04196,-0.0008962}, - /*Y*/{6.407,7.234,8907,367.8,0.005889,0.01988,2762,-14.36,5.397,-0.7269,0.04219,-0.000898}, - /*Zr*/{6.734,7.603,9120,405.2,0.005765,0.02039,2704,-14.59,5.463,-0.7333,0.04242,-0.0008998}, - /*Nb*/{6.902,7.791,9333,442.7,0.005587,0.0209,2621,-16.22,6.094,-0.8225,0.04791,-0.001024}, - /*Mo*/{6.425,7.248,9545,480.2,0.005367,0.02141,2517,-17.85,6.725,-0.9116,0.05339,-0.001148}, - /*Tc*/{6.799,7.671,9756,517.6,0.005315,0.02192,2493,-17.96,6.752,-0.9135,0.05341,-0.001147}, - /*Ru*/{6.108,6.887,9966,555.1,0.005151,0.02243,2416,-18.07,6.779,-0.9154,0.05342,-0.001145}, - /*Rh*/{5.924,6.677,1.018E4,592.5,0.004919,0.02294,2307,-18.18,6.806,-0.9173,0.05343,-0.001143}, - /*Pd*/{5.238,5.9,1.038E4,630,0.004758,0.02345,2231,-18.28,6.833,-0.9192,0.05345,-0.001142}, - /*Ag*/{5.623,6.354,7160,337.6,0.01394,0.02396,2193,-18.39,6.86,-0.9211,0.05346,-0.00114}, - /*Cd*/{5.814,6.554,1.08E4,355.5,0.004626,0.02447,2170,-18.62,6.915,-0.9243,0.0534,-0.001134}, - /*In*/{6.23,7.024,1.101E4,370.9,0.00454,0.02498,2129,-18.85,6.969,-0.9275,0.05335,-0.001127}, - /*Sn*/{6.41,7.227,1.121E4,386.4,0.004474,0.02549,2099,-19.07,7.024,-0.9308,0.05329,-0.001121}, - /*Sb*/{7.5,8.48,8608,348,0.009074,0.026,2069,-19.57,7.225,-0.9603,0.05518,-0.001165}, - /*Te*/{6.979,7.871,1.162E4,392.4,0.004402,0.02651,2065,-20.07,7.426,-0.9899,0.05707,-0.001209}, - /*I*/{7.725,8.716,1.183E4,394.8,0.004376,0.02702,2052,-20.56,7.627,-1.019,0.05896,-0.001254}, - /*Xe*/{8.231,9.289,1.203E4,397.3,0.004384,0.02753,2056,-21.06,7.828,-1.049,0.06085,-0.001298}, - /*Cs*/{7.287,8.218,1.223E4,399.7,0.004447,0.02804,2086,-20.4,7.54,-1.004,0.05782,-0.001224}, - /*Ba*/{7.899,8.911,1.243E4,402.1,0.004511,0.02855,2116,-19.74,7.252,-0.9588,0.05479,-0.001151}, - /*La*/{8.041,9.071,1.263E4,404.5,0.00454,0.02906,2129,-19.08,6.964,-0.9136,0.05176,-0.001077}, - /*Ce*/{7.489,8.444,1.283E4,406.9,0.00442,0.02957,2073,-18.43,6.677,-0.8684,0.04872,-0.001003}, - /*Pr*/{7.291,8.219,1.303E4,409.3,0.004298,0.03008,2016,-17.77,6.389,-0.8233,0.04569,-0.0009292}, - /*Nd*/{7.098,8,1.323E4,411.8,0.004182,0.03059,1962,-17.11,6.101,-0.7781,0.04266,-0.0008553}, - /*Pm*/{6.91,7.786,1.343E4,414.2,0.00405,0.0311,1903,-16.45,5.813,-0.733,0.03963,-0.0007815}, - /*Sm*/{6.728,7.58,1.362E4,416.6,0.003976,0.03161,1865,-15.79,5.526,-0.6878,0.0366,-0.0007077}, - /*Eu*/{6.551,7.38,1.382E4,419,0.003877,0.03212,1819,-15.13,5.238,-0.6426,0.03357,-0.0006339}, - /*Gd*/{6.739,7.592,1.402E4,421.4,0.003863,0.03263,1812,-14.47,4.95,-0.5975,0.03053,-0.0005601}, - /*Tb*/{6.212,6.996,1.421E4,423.9,0.003725,0.03314,1747,-14.56,4.984,-0.6022,0.03082,-0.0005668}, - /*Dy*/{5.517,6.21,1.44E4,426.3,0.003632,0.03365,1703,-14.65,5.018,-0.6069,0.03111,-0.0005734}, - /*Ho*/{5.219,5.874,1.46E4,428.7,0.003498,0.03416,1640,-14.74,5.051,-0.6117,0.03141,-0.0005801}, - /*Er*/{5.071,5.706,1.479E4,433,0.003405,0.03467,1597,-14.83,5.085,-0.6164,0.0317,-0.0005867}, - /*Tm*/{4.926,5.542,1.498E4,433.5,0.003342,0.03518,1567,-14.91,5.119,-0.6211,0.03199,-0.0005933}, - /*Yb*/{4.787, 5.386,1.517E4,435.9,0.003292,0.03569,1544,-15,5.153,-0.6258,0.03228,-0.0006}, - /*Lu*/{4.893, 5.505,1.536E4,438.4,0.003243,0.0362,1521,-15.09,5.186,-0.6305,0.03257,-0.0006066}, - /*Hf*/{5.028, 5.657,1.555E4,440.8,0.003195,0.03671,1499,-15.18,5.22,-0.6353,0.03286,-0.0006133}, - /*Ta*/{4.738, 5.329,1.574E4,443.2,0.003186,0.03722,1494,-15.27,5.254,-0.64,0.03315,-0.0006199}, - /*W*/{4.574, 5.144,1.593E4,442.4,0.003144,0.03773,1475,-15.67,5.392,-0.6577,0.03418,-0.0006426}, - /*Re*/{5.2, 5.851,1.612E4,441.6,0.003122,0.03824,1464,-16.07,5.529,-0.6755,0.03521,-0.0006654}, - /*Os*/{5.07, 5.704,1.63E4,440.9,0.003082,0.03875,1446,-16.47,5.667,-0.6932,0.03624,-0.0006881}, - /*Ir*/{4.945, 5.563,1.649E4,440.1,0.002965,0.03926,1390,-16.88,5.804,-0.711,0.03727,-0.0007109}, - /*Pt*/{4.476, 5.034,1.667E4,439.3,0.002871,0.03977,1347,-17.28,5.942,-0.7287,0.0383,-0.0007336}, - /*Au*/{4.856, 5.46,1.832E4,438.5,0.002542,0.04028,1354,-17.02,5.846,-0.7149,0.0374,-0.0007114}, - /*Hg*/{4.308, 4.843,1.704E4,487.8,0.002882,0.04079,1352,-17.84,6.183,-0.7659,0.04076,-0.0007925}, - /*Tl*/{4.723, 5.311,1.722E4,537,0.002913,0.0413,1366,-18.66,6.52,-0.8169,0.04411,-0.0008737}, - /*Pb*/{5.319, 5.982,1.74E4,586.3,0.002871,0.04181,1347,-19.48,6.857,-0.8678,0.04747,-0.0009548}, - /*Bi*/{5.956, 6.7,1.78E4,677,0.00266,0.04232,1336,-19.55,6.871,-0.8686,0.04748,-0.0009544}, - /*Po*/{6.158, 6.928,1.777E4,586.3,0.002812,0.04283,1319,-19.62,6.884,-0.8694,0.04748,-0.000954}, - /*At*/{6.204, 6.979,1.795E4,586.3,0.002776,0.04334,1302,-19.69,6.898,-0.8702,0.04749,-0.0009536}, - /*Rn*/{6.181, 6.954,1.812E4,586.3,0.002748,0.04385,1289,-19.76,6.912,-0.871,0.04749,-0.0009532}, - /*Fr*/{6.949, 7.82,1.83E4,586.3,0.002737,0.04436,1284,-19.83,6.926,-0.8718,0.0475,-0.0009528}, - /*Ra*/{7.506, 8.448,1.848E4,586.3,0.002727,0.04487,1279,-19.9,6.94,-0.8726,0.04751,-0.0009524}, - /*Ac*/{7.649, 8.609,1.866E4,586.3,0.002697,0.04538,1265,-19.97,6.953,-0.8733,0.04751,-0.000952}, - /*Th*/{7.71, 8.679,1.883E4,586.3,0.002641,0.04589,1239,-20.04,6.967,-0.8741,0.04752,-0.0009516}, - /*Pa*/{7.407, 8.336,1.901E4,586.3,0.002603,0.0464,1221,-20.11,6.981,-0.8749,0.04752,-0.0009512}, - /*U*/{7.29, 8.204,1.918E4,586.3,0.002573,0.04691,1207,-20.18,6.995,-0.8757,0.04753,-0.0009508} - }; -} -#endif diff --git a/include/EnergyLoss.h b/include/EnergyLoss.h deleted file mode 100644 index 969e524..0000000 --- a/include/EnergyLoss.h +++ /dev/null @@ -1,58 +0,0 @@ -/* - -EnergyLoss.h -Code for calculating the energy loss of a charged, massive particle through an arbitrary medium. -Based on code written by D.W. Visser at Yale for the original SPANC. Uses energy loss calulations -described by Ziegler in various SRIM textbooks. - -Written by G.W. McCann Aug. 2020 - -*/ - -#ifndef ENERGYLOSS_H -#define ENERGYLOSS_H - -#include -#include -#include -#include -#include -#include "MassLookup.h" - -namespace Mask { - class EnergyLoss { - - public: - EnergyLoss(); - ~EnergyLoss(); - double GetEnergyLoss(int zp, int ap, double e_initial, double thickness); - double GetReverseEnergyLoss(int zp, int ap, double e_final, double thickness); - double GetRange(double energy); - void SetTargetComponents(const std::vector& Zt, const std::vector& At, const std::vector& Stoich); - - private: - double GetElectronicStoppingPower(double energy); - double GetNuclearStoppingPower(double energy); - double GetTotalStoppingPower(double energy); - double Hydrogen_dEdx_Low(double e_per_u, int z); - double Hydrogen_dEdx_Med(double e_per_u, int z); - double Hydrogen_dEdx_High(double e_per_u, double energy, int z); - double CalculateEffectiveChargeRatio(double e_per_u, int z); - - int ZP, AP; - double MP; //units of u, isotopic - double comp_denom; - std::vector ZT, AT; - std::vector targ_composition; - - //constants for calculations - static constexpr double MAX_FRACTIONAL_STEP = 0.001; - static constexpr double MAX_DEPTH = 50; - static constexpr double MAX_H_E_PER_U = 100000.0; - static constexpr double AVOGADRO = 0.60221367; //N_A times 10^(-24) for converting - static constexpr double MEV2U = 1.0/931.4940954; - static constexpr double H_RESTMASS = 938.27231; //MeV, for beta calc - }; -} - -#endif diff --git a/include/MaskFile.h b/include/MaskFile.h deleted file mode 100644 index 42ffced..0000000 --- a/include/MaskFile.h +++ /dev/null @@ -1,77 +0,0 @@ -#ifndef MASKFILE_H -#define MASKFILE_H - -#include -#include -#include - -#include "Nucleus.h" -#include "RxnType.h" - -namespace Mask { - - struct MaskFileHeader { - RxnType rxn_type = RxnType::None; - int nsamples = -1; - }; - - struct MaskFileData { - std::vector E, KE, p, theta, phi; //ordered: target, (if not decay)projectile, ejectile, residual, break1... - std::vector Z, A; - std::vector detect_flag; - bool eof = false; //flag on end of file - }; - - class MaskFile { - public: - enum class FileType { - read, - write, - append, - none - }; - - MaskFile(); - MaskFile(const std::string& name, MaskFile::FileType type); - bool Open(const std::string& name, MaskFile::FileType type); - inline bool IsOpen() { return file.is_open(); } - void Close(); - - void WriteHeader(RxnType rxn_type, int nsamples); - void WriteData(const std::vector& data); - void WriteData(const MaskFileData& data); - MaskFileHeader ReadHeader(); - MaskFileData ReadData(); - - - - private: - - FileType file_type; - std::string filename; - uint64_t buffer_position; - uint64_t buffer_end; - uint32_t data_size; - RxnType m_rxn_type; - uint32_t buffersize_bytes; - - std::fstream file; - - std::vector data_buffer; - - static constexpr uint32_t onestep_rxn_n = 2; - static constexpr uint32_t twostep_rxn_n = 4; - static constexpr uint32_t threestep_rxn_n = 6; - - static constexpr uint64_t buffersize = 10000; //number of events - static constexpr int width = 0; - static constexpr int precision = 3; - - static constexpr std::size_t double_size = sizeof(double); - static constexpr std::size_t int_size = sizeof(uint32_t); - static constexpr std::size_t bool_size = sizeof(bool); - }; - -}; - -#endif \ No newline at end of file diff --git a/include/Rotation.h b/include/Rotation.h deleted file mode 100644 index 9b8ba3c..0000000 --- a/include/Rotation.h +++ /dev/null @@ -1,78 +0,0 @@ -/* - Classes which define rotations about the x, y, and z axes. Using these, - any arbitrary orientation can be described. Methods implemented for vector multiplication - as well as generating the inverse of the rotation. -*/ -#ifndef ROTATION_H -#define ROTATION_H - -#include "Vec3.h" - -namespace Mask { - - class XRotation { - public: - XRotation(); - XRotation(double ang); - ~XRotation(); - Vec3 Rotate(const Vec3& vector); - inline void SetAngle(double ang) { m_angle = ang; GenerateMatrix(); } - inline XRotation GetInverse() { return XRotation(-m_angle); } - inline Vec3 operator*(const Vec3& 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 Vec3(x, y, z); - } - - private: - void GenerateMatrix(); - double m_angle; - double m_matrix[3][3]; - }; - - class YRotation { - public: - YRotation(); - YRotation(double ang); - ~YRotation(); - Vec3 Rotate(const Vec3& vector); - inline void SetAngle(double ang) { m_angle = ang; GenerateMatrix(); } - inline YRotation GetInverse() { return YRotation(-m_angle); } - inline Vec3 operator*(const Vec3& 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 Vec3(x, y, z); - } - - private: - void GenerateMatrix(); - double m_angle; - double m_matrix[3][3]; - }; - - class ZRotation { - public: - ZRotation(); - ZRotation(double ang); - ~ZRotation(); - Vec3 Rotate(const Vec3& vector); - inline void SetAngle(double ang) { m_angle = ang; GenerateMatrix(); } - inline ZRotation GetInverse() { return ZRotation(-m_angle); } - inline Vec3 operator*(const Vec3& 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 Vec3(x, y, z); - } - - private: - void GenerateMatrix(); - double m_angle; - double m_matrix[3][3]; - }; - -} - -#endif \ No newline at end of file diff --git a/include/Vec3.h b/include/Vec3.h deleted file mode 100644 index dcd1ad8..0000000 --- a/include/Vec3.h +++ /dev/null @@ -1,61 +0,0 @@ -/* - Class to represent a 3-space vector in both cartesian and spherical coordinates. Can perform vector - addition, subtraction, and dot product. - - --GWM Dec 2020 -*/ -#ifndef VEC3_H -#define VEC3_H - -#include - -namespace Mask { - - class Vec3 { - public: - Vec3(); - Vec3(double x, double y, double z); - ~Vec3(); - - 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 Vec3& operator=(const Vec3& rhs) { SetVectorCartesian(rhs.GetX(), rhs.GetY(), rhs.GetZ()); return *this; } - inline Vec3 operator+(const Vec3& rhs) const { return Vec3(this->GetX()+rhs.GetX(), this->GetY()+rhs.GetY(), this->GetZ()+rhs.GetZ()); } - inline Vec3 operator-(const Vec3& rhs) const { return Vec3(this->GetX()-rhs.GetX(), this->GetY()-rhs.GetY(), this->GetZ()-rhs.GetZ()); } - - - double Dot(const Vec3& rhs) const; - Vec3 Cross(const Vec3& rhs) const; - - - - private: - - //Use instead of std::atan2. Better control over values close to x=0 - 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 \ No newline at end of file diff --git a/include/Vec4.h b/include/Vec4.h deleted file mode 100644 index 51eb879..0000000 --- a/include/Vec4.h +++ /dev/null @@ -1,69 +0,0 @@ -/* - Class which represents a 4-momentum vector. Can perform vector addition, subtraction, dot product - and generate a boost vector to its rest frame as well as apply a boost to itself. - - --GWM Dec 2020. -*/ -#ifndef VEC4_H -#define VEC4_H - -#include - -namespace Mask { - - class Vec4 { - public: - Vec4(); - Vec4(double px, double py, double pz, double E); - virtual ~Vec4(); - void SetVectorCartesian(double px, double py, double pz, double E); - void SetVectorSpherical(double theta, double phi, double p, double E); - - inline double GetE() const { return m_data[3]; } - inline double GetPx() const { return m_data[0]; } - inline double GetPy() const { return m_data[1]; } - 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 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 { - double phi = Atan2(GetPy(), GetPx()); - if(phi<0) phi += 2.0*M_PI; - return GetPx() == 0.0 && GetPy() == 0.0 ? 0.0 : phi; - } - - inline double GetInvMass() const { return std::sqrt(GetE()*GetE() - GetP()*GetP()); } - inline double GetKE() const { return GetE() - GetInvMass(); } - inline const double* GetBoost() const { return &m_boost[0]; } - - void ApplyBoost(const double* boost); - - //Only intended for use in looping access! - inline const double operator[] (int index) const { return index>3 || index < 0 ? 0.0 : m_data[index]; } - - inline Vec4& operator=(const Vec4& rhs) { SetVectorCartesian(rhs.GetPx(), rhs.GetPy(), rhs.GetPz(), rhs.GetE()); return *this; } - inline Vec4 operator+(const Vec4& rhs) const { return Vec4(rhs.GetPx()+GetPx(), rhs.GetPy()+GetPy(), rhs.GetPz()+GetPz(), rhs.GetE()+GetE()); } - inline Vec4 operator-(const Vec4& rhs) const { return Vec4(rhs.GetPx()-GetPx(), rhs.GetPy()-GetPy(), rhs.GetPz()-GetPz(), rhs.GetE()-GetE()); } - - double Dot(const Vec4& rhs) const; - Vec4 Cross(const Vec4& rhs) const; - - private: - void CalcBoostToCM(); - - //use instead of std::atan2. Better controll over x=0 - 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_boost[3]; - - }; - -} - -#endif diff --git a/src/Mask/EnergyLoss.cpp b/src/Mask/EnergyLoss.cpp deleted file mode 100644 index 2fdcbfb..0000000 --- a/src/Mask/EnergyLoss.cpp +++ /dev/null @@ -1,243 +0,0 @@ -/* - -EnergyLoss.cpp -Code for calculating the energy loss of a charged, massive particle through an arbitrary medium. -Based on code written by D.W. Visser at Yale for the original SPANC. Uses energy loss calulations -described by Ziegler in various SRIM textbooks. - -Written by G.W. McCann Aug. 2020 - -*/ - -#include "Eloss_Tables.h" -#include "EnergyLoss.h" -#include "KinematicsExceptions.h" -#include - -namespace Mask { - - EnergyLoss::EnergyLoss() : - ZP(-1), AP(-1), MP(-1.0), comp_denom(0) - { - } - - EnergyLoss::~EnergyLoss() {} - - /*Targets are defined by their atomic number, total number of nucleons, and their stoichiometry within the target compound*/ - void EnergyLoss::SetTargetComponents(const std::vector& Zt, const std::vector& At, const std::vector& Stoich) { - comp_denom = 0; - ZT = Zt; - AT = At; - for(unsigned int i=0; i MAX_Z) - throw ELossException("Maximum allowed target Z exceeded"); - } - targ_composition.resize(Stoich.size()); - - for(unsigned int i=0; i MAX_FRACTIONAL_STEP && depth < MAX_DEPTH) { - depth++; - x_step *= 0.5; - e_step = GetTotalStoppingPower(e_final)*x_step/1000.0; - } else if((x_step + x_traversed) >= thickness) { //last chunk - go = false; - x_step = thickness - x_traversed; //get valid portion of last chunk - e_final -= GetTotalStoppingPower(e_final)*x_step/1000.0; - if(e_final <= e_threshold) - return e_initial; - } else if(depth == MAX_DEPTH) { - return e_initial; - } else { - x_traversed += x_step; - e_step = GetTotalStoppingPower(e_final)*x_step/1000.0; - e_final -= e_step; - if(e_final <= e_threshold) - return e_initial; - } - } - return e_initial - e_final; - } - - /* - Returns units of MeV; thickness in ug/cm^2; e_final in units of MeV - Energy loss going through the target using energy of particle after traversal - */ - double EnergyLoss::GetReverseEnergyLoss(int zp, int ap, double e_final, double thickness) { - if( ZP != zp) { - ZP = zp; - AP = ap; - MP = MassLookup::GetInstance().FindMass(ZP, AP)*MEV2U; - - } - - double e_initial = e_final; - double x_traversed = 0.0; - double x_step = 0.25*thickness; //initial step in x - double e_step = GetTotalStoppingPower(e_initial)*x_step/1000.0; //initial step in E - - bool go = true; - while(go) { - if(e_step/e_initial > MAX_FRACTIONAL_STEP) { - x_step *= 0.5; - e_step = GetTotalStoppingPower(e_initial)*x_step/1000.0; - } else if (x_traversed+x_step > thickness) { - go = false; - x_step = thickness - x_traversed; - e_initial += GetTotalStoppingPower(e_initial)*x_step/1000.0; - } else { - x_traversed += x_step; - e_step = GetTotalStoppingPower(e_initial)*x_step/1000.0; - e_initial += e_step; - } - } - - return e_initial-e_final; - } - - /* - Returns units of keV/(ug/cm^2) - Calculates Electronic stopping power by first calculating SE for hydrogen through the target and then using - corrections to calculate SE for projectile of interest - */ - double EnergyLoss::GetElectronicStoppingPower(double energy) { - //Wants in units of keV - energy *= 1000.0; - double e_per_u = energy/MP; - std::vector values; - if(e_per_u > MAX_H_E_PER_U) { - throw ELossException("Exceeded maximum allowed energy per nucleon"); - } else if (e_per_u > 1000.0) { - for(auto& z: ZT) - values.push_back(Hydrogen_dEdx_High(e_per_u, energy, z)); - } else if (e_per_u > 10.0) { - for(auto& z: ZT) - values.push_back(Hydrogen_dEdx_Med(e_per_u, z)); - } else if (e_per_u > 0.0) { - for(auto& z: ZT) - values.push_back(Hydrogen_dEdx_Low(e_per_u, z)); - } else { - throw ELossException("Negative energy per nucleon; given energy: "+std::to_string(energy)); - } - - if(values.size() == 0) - throw ELossException("Size of value array is 0. Unable to iterate over target components"); - - if(ZP > 1) { //not hydrogen, need to account for effective charge - for(unsigned int i=0; i -#include -#include - -/* - - FORMAT - - HEADER (contains rxntype & nuclei numbers & beam kinetic energy) (64 bits, 8 bytes) - NSAMPLES(32bit) RXNTYPE(32bit) - end HEADER - - There are NSAMPLES * (number of saved nuclei) data in the file. The number of nuclei saved is related to the - RXNTYPE. All nuclei (target, projectile, ejectile, residual, break1, etc...) are saved. A datum is as follows: - - DATUM (contains kinematic data for a nucleus) (384 bits, 48 bytes) - Z(32bit) A(32bit) DETECTFLAG(32bit) E(64bit) KE(64bit) P(64bit) THETA(64bit) PHI(64bit) - end DATUM - - -*/ - -namespace Mask { - - MaskFile::MaskFile() : - file_type(MaskFile::FileType::none), filename(""), buffer_position(0), buffer_end(0), data_size(0), buffersize_bytes(0), file() - { - } - - MaskFile::MaskFile(const std::string& name, MaskFile::FileType type) : - file_type(type), filename(name), buffer_position(0), buffer_end(0), data_size(0), buffersize_bytes(0), file() - { - Open(filename, type); - } - - bool MaskFile::Open(const std::string& name, MaskFile::FileType type) { - if(IsOpen()) { - std::cerr<<"Attempted to open file that is already open!"< 0 && buffer_position < buffersize_bytes && (file_type == FileType::write || file_type == FileType::append)) { - file.write(data_buffer.data(), buffer_position); - } - file.close(); - } - - void MaskFile::WriteHeader(RxnType rxn_type, int nsamples) { - - m_rxn_type = rxn_type; - switch(rxn_type) - { - case RxnType::PureDecay : - { - data_size = 3 * ( 5 * double_size + 2 * int_size + bool_size); - break; - } - case RxnType::OneStepRxn : - { - data_size = 4 * ( 5 * double_size + 2 * int_size + bool_size); - break; - } - case RxnType::TwoStepRxn : - { - data_size = 6 * ( 5 * double_size + 2 * int_size + bool_size); - break; - } - case RxnType::ThreeStepRxn : - { - data_size = 8 * ( 5 * double_size + 2 * int_size + bool_size); - break; - } - case RxnType::None : - { - std::cerr<<"Invalid RxnType at MaskFile::WriteHeader!"< temp_buffer(4); - file.read(temp_buffer.data(), 4); - header.nsamples = *(int*)(&temp_buffer[0]); - file.read(temp_buffer.data(), 4); - uint32_t type_value = *(uint32_t*)(&temp_buffer[0]); - m_rxn_type = GetRxnTypeFromInt(type_value); - - switch(m_rxn_type) - { - case RxnType::PureDecay: - { - data_size = 3 * ( 5 * double_size + 2 * int_size + bool_size); - break; - } - case RxnType::OneStepRxn: - { - data_size = 4 * ( 5 * double_size + 2 * int_size + bool_size); - break; - } - case RxnType::TwoStepRxn: - { - data_size = 6 * ( 5 * double_size + 2 * int_size + bool_size); - break; - } - case RxnType::ThreeStepRxn: - { - data_size = 8 * ( 5 * double_size + 2 * int_size + bool_size); - break; - } - case RxnType::None: - { - std::cerr<<"Unexpected reaction type at MaskFile::ReadHeader (rxn type = "<& data) { - - char* data_pointer; - double datum; - int number; - bool flag; - std::size_t j; - for(unsigned int i=0; i buffer_end) { - std::cerr<<"Attempting to read past end of file at MaskFile::ReadData! Returning empty"<A) return false; - - m_z = Z; - m_a = A; - m_gs_mass = MassLookup::GetInstance().FindMass(Z, A); - m_symbol = MassLookup::GetInstance().FindSymbol(Z, A); - SetVectorCartesian(0,0,0,m_gs_mass); - return true; - } - -} \ No newline at end of file diff --git a/src/Mask/Rotation.cpp b/src/Mask/Rotation.cpp deleted file mode 100644 index aa07ad8..0000000 --- a/src/Mask/Rotation.cpp +++ /dev/null @@ -1,71 +0,0 @@ -/* - Classes which define rotations about the x, y, and z axes. Using these, - any arbitrary orientation can be described. Methods implemented for vector multiplication - as well as generating the inverse of the rotation. -*/ -#include "Rotation.h" - -namespace Mask { - - XRotation::XRotation() : - m_angle(0) - { - GenerateMatrix(); - } - - XRotation::XRotation(double angle) : - m_angle(angle) - { - GenerateMatrix(); - } - - XRotation::~XRotation() {} - - void XRotation::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); - } - - YRotation::YRotation() : - m_angle(0) - { - GenerateMatrix(); - } - - YRotation::YRotation(double angle) : - m_angle(angle) - { - GenerateMatrix(); - } - - YRotation::~YRotation() {} - - void YRotation::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); - } - - - ZRotation::ZRotation() : - m_angle(0) - { - GenerateMatrix(); - } - - ZRotation::ZRotation(double angle) : - m_angle(angle) - { - GenerateMatrix(); - } - - ZRotation::~ZRotation() {} - - void ZRotation::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; - } - -}; \ No newline at end of file diff --git a/src/Mask/Vec3.cpp b/src/Mask/Vec3.cpp deleted file mode 100644 index 42abdc9..0000000 --- a/src/Mask/Vec3.cpp +++ /dev/null @@ -1,48 +0,0 @@ -/* - Class to represent a 3-space vector in both cartesian and spherical coordinates. Can perform vector - addition, subtraction, and dot product. - - --GWM Dec 2020 -*/ -#include "Vec3.h" - -namespace Mask { - - Vec3::Vec3() { - m_data[0] = 0.; - m_data[1] = 0.; - m_data[2] = 0.; - } - - Vec3::Vec3(double x, double y, double z) { - m_data[0] = x; - m_data[1] = y; - m_data[2] = z; - } - - Vec3::~Vec3() {} - - void Vec3::SetVectorCartesian(double x, double y, double z) { - m_data[0] = x; - m_data[1] = y; - m_data[2] = z; - } - - void Vec3::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 Vec3::Dot(const Vec3& rhs) const { - return GetX()*rhs.GetX() + GetY()*rhs.GetY() + GetZ()*rhs.GetZ(); - } - - Vec3 Vec3::Cross(const Vec3& rhs) const { - double x = GetY()*rhs.GetZ() - GetZ()*rhs.GetY(); - double y = GetZ()*rhs.GetX() - GetX()*rhs.GetZ(); - double z = GetX()*rhs.GetY() - GetY()*rhs.GetX(); - return Vec3(x,y,z); - } - -} diff --git a/src/Mask/Vec4.cpp b/src/Mask/Vec4.cpp deleted file mode 100644 index 493ea3e..0000000 --- a/src/Mask/Vec4.cpp +++ /dev/null @@ -1,73 +0,0 @@ -/* - Class which represents a 4-momentum vector. Can perform vector addition, subtraction, dot product - and generate a boost vector to its rest frame as well as apply a boost to itself. - - --GWM Dec 2020. - NOTE: uses (-,-,-,+) metric (same as ROOT convention) -*/ -#include "Vec4.h" - - -namespace Mask { - - Vec4::Vec4() { - for(auto& val: m_data) - val = 0.0; - for(auto& val: m_boost) - val = 0.0; - } - - Vec4::Vec4(double px, double py, double pz, double E) { - m_data[0] = px; - m_data[1] = py; - m_data[2] = pz; - m_data[3] = E; - CalcBoostToCM(); - } - - Vec4::~Vec4() {} - - void Vec4::SetVectorCartesian(double px, double py, double pz, double E) { - m_data[0] = px; - m_data[1] = py; - m_data[2] = pz; - m_data[3] = E; - - CalcBoostToCM(); - } - - void Vec4::SetVectorSpherical(double theta, double phi, double p, double E) { - m_data[0] = p*cos(phi)*sin(theta); - m_data[1] = p*sin(phi)*sin(theta); - m_data[2] = p*cos(theta); - m_data[3] = E; - CalcBoostToCM(); - } - - void Vec4::CalcBoostToCM() { - m_boost[0] = m_data[0]/m_data[3]; - m_boost[1] = m_data[1]/m_data[3]; - m_boost[2] = m_data[2]/m_data[3]; - } - - void Vec4::ApplyBoost(const double* beta) { - double beta2 = beta[0]*beta[0] + beta[1]*beta[1] + beta[2]*beta[2]; - double gamma = 1.0/std::sqrt(1.0 - beta2); - double bdotp = beta[0]*m_data[0] + beta[1]*m_data[1] + beta[2]*m_data[2]; - double gfactor = beta2>0.0 ? (gamma - 1.0)/beta2 : 0.0; - - SetVectorCartesian(GetPx()+gfactor*bdotp*beta[0]+gamma*beta[0]*GetE(), - GetPy()+gfactor*bdotp*beta[1]+gamma*beta[1]*GetE(), - GetPz()+gfactor*bdotp*beta[2]+gamma*beta[2]*GetE(), - gamma*(GetE() + bdotp)); - } - - double Vec4::Dot(const Vec4& rhs) const { - return GetE()*rhs.GetE() - GetPx()*rhs.GetPx() - GetPy()*rhs.GetPy() - GetPz()*rhs.GetPz(); - } - - Vec4 Vec4::Cross(const Vec4& rhs) const { - return Vec4(); - } - -} \ No newline at end of file diff --git a/src/Plotters/Python/.gitignore b/src/Plotters/Python/.gitignore deleted file mode 100644 index 231854d..0000000 --- a/src/Plotters/Python/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -__pycache__ - -!.gitignore \ No newline at end of file diff --git a/src/Plotters/Python/MaskFile.py b/src/Plotters/Python/MaskFile.py deleted file mode 100755 index 210d8f9..0000000 --- a/src/Plotters/Python/MaskFile.py +++ /dev/null @@ -1,80 +0,0 @@ -#!/usr/bin/env python3 - -import numpy as np -import struct - -class MaskFileData : - def __init__(self, n): - self.Z = np.zeros(n, dtype=int) - self.A = np.zeros(n, dtype=int) - self.dFlag = np.zeros(n, dtype=bool) - self.E = np.zeros(n) - self.KE = np.zeros(n) - self.p = np.zeros(n) - self.theta = np.zeros(n) - self.phi = np.zeros(n) - -class MaskFile: - int_size = 4 - double_size = 8 - bool_size = 1 - - def __init__(self, filename=""): - self.eofFlag = False - self.openFlag = False - if filename != "" : - self.Open(filename) - - def Open(self, filename): - self.filename = filename - self.file = open(self.filename, mode="rb") - if self.file.closed : - self.openFlag = False - else: - self.openFlag = True - - def ReadHeader(self): - data = self.file.read(2*self.int_size) - - (self.nsamples, self.rxntype) = struct.unpack("ii", data) - self.datasize = (5*self.double_size+2*self.int_size+self.bool_size) - self.datastr = "=ii?ddddd" - - if self.rxntype == 0: - self.N_nuclei = 3 - elif self.rxntype == 1: - self.N_nuclei = 4 - elif self.rxntype == 2: - self.N_nuclei = 6 - elif self.rxntype == 3: - self.N_nuclei = 8 - - def ReadData(self): - data = MaskFileData(self.N_nuclei) - for i in range(self.N_nuclei): - buffer = self.file.read(self.datasize) - (data.Z[i], data.A[i], data.dFlag[i], data.E[i], data.KE[i], data.p[i], data.theta[i], data.phi[i]) = struct.unpack(self.datastr, buffer) - if buffer == "": - self.eofFlag = True - - return data - - def Close(self): - self.file.close() - -def main() : - file = MaskFile(filename="/data1/gwm17/mask_tests/7Bedp_870keV_beam_50CD2.mask") - file.ReadHeader() - print("samples: ", file.nsamples, "rxntype:", file.rxntype, "datasize:", file.datasize) - count=0 - for i in range(file.nsamples): - file.ReadData() - count += 1 - - print("count:",count) - print("eofFlag:",file.eofFlag) - - file.Close() - -if __name__ == '__main__': - main() \ No newline at end of file diff --git a/src/Plotters/Python/NucData.py b/src/Plotters/Python/NucData.py deleted file mode 100755 index 7b4200d..0000000 --- a/src/Plotters/Python/NucData.py +++ /dev/null @@ -1,70 +0,0 @@ -#!/usr/bin/env python3 - -import numpy as np -import requests -import lxml.html as xhtml - - -class MassTable: - def __init__(self): - file = open("./etc/mass.txt","r") - self.mtable = {} - u2mev = 931.4940954 - me = 0.000548579909 #amu - self.etable = {} - - line = file.readline() - line = file.readline() - - for line in file: - entries = line.split() - n = entries[0] - z = entries[1] - a = entries[2] - element = entries[3] - massBig = float(entries[4]) - massSmall = float(entries[5]) - - key = '('+z+','+a+')' - value = ((massBig+massSmall*1e-6) - float(z)*me)*u2mev - self.mtable[key] = value - self.etable[key] = element - file.close() - - def GetMass(self, z, a): - key = '('+str(z)+','+str(a)+')' - if key in self.mtable: - return self.mtable[key] - else: - return 0 - - def GetSymbol(self, z, a): - key = '('+str(z)+','+str(a)+')' - if key in self.etable: - return str(a)+self.etable[key] - else: - return 'none' - -Masses = MassTable() - -def GetExcitations(symbol): - levels = np.array(np.empty(0)) - text = '' - - site = requests.get("https://www.nndc.bnl.gov/nudat2/getdatasetClassic.jsp?nucleus="+symbol+"&unc=nds") - contents = xhtml.fromstring(site.content) - tables = contents.xpath("//table") - rows = tables[2].xpath("./tr") - for row in rows[1:-2]: - entries = row.xpath("./td") - if len(entries) != 0: - entry = entries[0] - data = entry.xpath("./a") - if len(data) == 0: - text = entry.text - else: - text = data[0].text - text = text.replace('?', '') - text = text.replace('\xa0\xa0≈','') - levels = np.append(levels, float(text)/1000.0) - return levels diff --git a/src/Plotters/Python/Nucleus.py b/src/Plotters/Python/Nucleus.py deleted file mode 100755 index a77a087..0000000 --- a/src/Plotters/Python/Nucleus.py +++ /dev/null @@ -1,106 +0,0 @@ -#!/usr/bin/env python3 - -import numpy as np -from NucData import Masses - -class Nucleus: - deg2rad = np.pi/180.0 - def __init__(self, Z=0, A=0): - self.Z = Z - self.A = A - if Z != 0 and A != 0: - self.gsMass = Masses.GetMass(self.Z, self.A) - self.vec4 = np.zeros(4) - self.symbol = Masses.GetSymbol(self.Z, self.A) - self.vec4[3] = self.gsMass - - def SetIsotope(self, Z, A): - self.gsMass = Masses.GetMass(Z, A) - self.symbol = Masses.GetSymbol(Z, A) - self.Z = Z - self.A = A - self.vec4 = np.zeros(4) - self.vec4[3] = self.gsMass - - def SetVectorCart(self, px, py, pz, E): - self.vec4[0] = px - self.vec4[1] = py - self.vec4[2] = pz - self.vec4[3] = E - - def SetVectorSpher(self, theta, phi, p, E): - self.vec4[0] = p*np.sin(theta)*np.cos(phi) - self.vec4[1] = p*np.sin(theta)*np.sin(phi) - self.vec4[2] = p*np.cos(theta) - self.vec4[3] = E - - def __add__(self, other): - vec4 = self.vec4 + other.vec4 - newNuc = Nucleus(self.Z + other.Z, self.A + other.A) - newNuc.SetVectorCart(vec4[0], vec4[1], vec4[2], vec4[3]) - return newNuc - - def __sub__(self, other): - vec4 = self.vec4 - other.vec4 - newNuc = Nucleus(self.Z - other.Z, self.A - other.A) - newNuc.SetVectorCart(vec4[0], vec4[1], vec4[2], vec4[3]) - return newNuc - - def __str__(self): - return "Nucleus({0},{1}) with 4-vector({2})".format(self.Z, self.A, self.vec4) - - def GetP(self): - return np.sqrt(self.vec4[0]**2.0 + self.vec4[1]**2.0 + self.vec4[2]**2.0) - - def GetInvMass(self): - return np.sqrt(self.vec4[3]**2.0 - self.GetP()**2.0) - - def GetKE(self): - return self.vec4[3] - self.GetInvMass() - - def GetTheta(self): - return np.arccos(self.vec4[2]/self.GetP()) - - def GetPhi(self): - result = np.arctan2(self.vec4[1], self.vec4[0]) - if result < 0.0: - result += 2.0*np.pi - return result - - def GetExcitation(self): - return self.GetInvMass() - self.gsMass - - def GetBoostToCMFrame(self): - boost_vec = np.zeros(3) - boost_vec[0] = self.vec4[0]/self.vec4[3] - boost_vec[1] = self.vec4[1]/self.vec4[3] - boost_vec[2] = self.vec4[2]/self.vec4[3] - return boost_vec - - def ApplyBoost(self, boost_vec): - beta2 = np.linalg.norm(boost_vec)**2.0 - gamma = 1.0/np.sqrt(1.0 - beta2) - bdotp = boost_vec[0]*self.vec4[0] + boost_vec[1]*self.vec4[1] + boost_vec[2]*self.vec4[2] - gfactor = (gamma-1.0)/beta2 if beta2 > 0.0 else 0.0 - - px = self.vec4[0]+gfactor*bdotp*boost_vec[0]+gamma*boost_vec[0]*self.vec4[3] - py = self.vec4[1]+gfactor*bdotp*boost_vec[1]+gamma*boost_vec[1]*self.vec4[3] - pz = self.vec4[2]+gfactor*bdotp*boost_vec[2]+gamma*boost_vec[2]*self.vec4[3] - E = gamma*(self.vec4[3] + bdotp) - - self.SetVectorCart(px, py, pz, E); - - -def main(): - nuc = Nucleus(1,1) - print("First",nuc) - nuc2 = Nucleus(2,4) - print("Second", nuc2) - result = nuc + nuc2 - print("Addition", result) - result2 = nuc2 - nuc - print("Subtraction", result2) - - -if __name__ == '__main__': - main() \ No newline at end of file diff --git a/src/Plotters/Python/PyPlotter.py b/src/Plotters/Python/PyPlotter.py deleted file mode 100755 index deb7042..0000000 --- a/src/Plotters/Python/PyPlotter.py +++ /dev/null @@ -1,135 +0,0 @@ -#!/usr/bin/env python3 - -import numpy as np -import matplotlib.pyplot as plt -from Nucleus import Nucleus -from MaskFile import MaskFile, MaskFileData -from NucData import Masses -import pickle -import sys - -def PlotData(inputname, outputname): - rad2deg = 180.0/np.pi - - datafile = MaskFile(inputname) - datafile.ReadHeader() - - print("MaskFile opened -- rxntype:", datafile.rxntype, "number of samples:", datafile.nsamples) - - data = MaskFileData(datafile.N_nuclei) - - ke = np.zeros((datafile.N_nuclei, datafile.nsamples)) - ke_d = np.zeros((datafile.N_nuclei, datafile.nsamples)) - theta = np.zeros((datafile.N_nuclei, datafile.nsamples)) - theta_d = np.zeros((datafile.N_nuclei, datafile.nsamples)) - phi = np.zeros((datafile.N_nuclei, datafile.nsamples)) - phi_d = np.zeros((datafile.N_nuclei, datafile.nsamples)) - detect_mask = np.ones((datafile.N_nuclei, datafile.nsamples), dtype=bool) - names = [] - for i in range(datafile.N_nuclei): - names.append(" ") - - nuc = Nucleus() - - for i in range(datafile.nsamples): - data = datafile.ReadData() - - if i == 0: - for j in range(datafile.N_nuclei): - names[j] = Masses.GetSymbol(data.Z[j], data.A[j]) - - for j in range(datafile.N_nuclei): - nuc.SetIsotope(data.Z[j], data.A[j]) - nuc.SetVectorSpher(data.theta[j], data.phi[j], data.p[j], data.E[j]) - - ke[j][i] = nuc.GetKE() - theta[j][i] = data.theta[j]*rad2deg - phi[j][i] = data.phi[j]*rad2deg - if data.dFlag[j] == True: - ke_d[j][i] = data.KE[j] - theta_d[j][i] = data.theta[j]*rad2deg - phi_d[j][i] = data.phi[j]*rad2deg - else: - detect_mask[j][i] = False - - datafile.Close() - - #Remove empty values from detection arrays - final_theta_d = theta_d[detect_mask] - final_phi_d = phi_d[detect_mask] - final_ke_d = ke_d[detect_mask] - - #figs = {} - #axes = {} - - ''' - for i in range(len(names)): - figs[i], axes[i] = plt.subplots(2,2) - ''' - fig, axes = plt.subplots(len(names)-1,4) - fig.set_size_inches(12, 12) - - for i in range(1, len(names)): - ''' - axes[i][0][0].plot(theta[i], ke[i], marker=',', linestyle='None') - axes[i][0][0].set_title(names[i]+" KE vs. Theta") - axes[i][0][0].set_xlabel(r"$\theta$ (degrees)") - axes[i][0][0].set_ylabel("KE (MeV)") - - axes[i][0][1].plot(phi[i], ke[i], marker=",", linestyle='None') - axes[i][0][1].set_title(names[i]+" KE vs. Phi") - axes[i][0][1].set_xlabel(r"$\phi$ (degrees)") - axes[i][0][1].set_ylabel("KE (MeV)") - - axes[i][1][0].plot(theta_d[i], ke_d[i], marker=',', linestyle='None') - axes[i][1][0].set_title(names[i]+" KE vs. Theta -- Detected") - axes[i][1][0].set_xlabel(r"$\theta$ (degrees)") - axes[i][1][0].set_ylabel("KE (MeV)") - - axes[i][1][1].plot(phi_d[i], ke_d[i], marker=",", linestyle='None') - axes[i][1][1].set_title(names[i]+" KE vs. Phi -- Detected") - axes[i][1][1].set_xlabel(r"$\phi$ (degrees)") - axes[i][1][1].set_ylabel("KE (MeV)") - ''' - - axes[i-1][0].plot(theta[i], ke[i], marker=',', linestyle='None') - axes[i-1][0].set_title(names[i]+" KE vs. Theta") - axes[i-1][0].set_xlabel(r"$\theta$ (degrees)") - axes[i-1][0].set_ylabel("KE (MeV)") - - axes[i-1][1].plot(phi[i], ke[i], marker=",", linestyle='None') - axes[i-1][1].set_title(names[i]+" KE vs. Phi") - axes[i-1][1].set_xlabel(r"$\phi$ (degrees)") - axes[i-1][1].set_ylabel("KE (MeV)") - - axes[i-1][2].plot(theta_d[i], ke_d[i], marker=',', linestyle='None') - axes[i-1][2].set_title(names[i]+" KE vs. Theta -- Detected") - axes[i-1][2].set_xlabel(r"$\theta$ (degrees)") - axes[i-1][2].set_ylabel("KE (MeV)") - - axes[i-1][3].plot(phi_d[i], ke_d[i], marker=",", linestyle='None') - axes[i-1][3].set_title(names[i]+" KE vs. Phi -- Detected") - axes[i-1][3].set_xlabel(r"$\phi$ (degrees)") - axes[i-1][3].set_ylabel("KE (MeV)") - - plt.tight_layout() - plt.show(block=True) - - print("Writing figure to file:", outputname) - with open(outputname, "wb") as outfile: - pickle.dump(fig, outfile) - outfile.close() - - print("Finished.") - -def main(): - if len(sys.argv) == 3: - PlotData(sys.argv[1], sys.argv[2]) - else: - print("Unable to run PyPlotter, incorrect number of arguments! Need an input datafile name, and an output plot file name") - -if __name__ == '__main__': - main() - - - diff --git a/src/Plotters/Python/ViewPyPlots.py b/src/Plotters/Python/ViewPyPlots.py deleted file mode 100755 index 3ed324d..0000000 --- a/src/Plotters/Python/ViewPyPlots.py +++ /dev/null @@ -1,32 +0,0 @@ -#!/usr/bin/env python3 - -import matplotlib.pyplot as plt -import numpy as np -import pickle -import sys - -def SetManager(figure): - dummy = plt.figure() - manager = dummy.canvas.manager - manager.canvas.figure = figure - figure.set_canvas(manager.canvas) - -def ViewPyPlots(filename): - - figure = pickle.load(open(filename, "rb")) - SetManager(figure) - - figure.set_size_inches(12, 12) - plt.tight_layout() - - plt.show(block=True) - - -def main(): - if len(sys.argv) == 2: - ViewPyPlots(sys.argv[1]) - else: - print("Unable to run ViewPyPlots, incorrect number of commandline arguments -- requires an input pickle file") - -if __name__ == '__main__': - main() \ No newline at end of file