diff --git a/include/AngularDistribution.h b/include/AngularDistribution.h index 393699f..05a7767 100644 --- a/include/AngularDistribution.h +++ b/include/AngularDistribution.h @@ -5,35 +5,39 @@ #include #include -class AngularDistribution { -public: - AngularDistribution(); - AngularDistribution(const std::string& file); - ~AngularDistribution(); - void ReadDistributionFile(const std::string& file); - void AttachRandomNumberGenerator(std::mt19937* random) { generator = random; }; - double GetRandomCosTheta(); - int GetL() { return L; }; - double GetBranchingRatio() { return branchingRatio; }; +namespace Mask { -private: - bool IsIsotropic(); - bool IsGeneratorSet() { - if(generator) { - return true; - } else { - return false; + class AngularDistribution { + public: + AngularDistribution(); + AngularDistribution(const std::string& file); + ~AngularDistribution(); + void ReadDistributionFile(const std::string& file); + inline void AttachRandomNumberGenerator(std::mt19937* random) { generator = random; } + double GetRandomCosTheta(); + inline int GetL() { return L; } + inline double GetBranchingRatio() { return branchingRatio; } + + private: + inline bool IsIsotropic() { return isoFlag; } + inline bool IsGeneratorSet() { + if(generator) { + return true; + } else { + return false; + } } - } + + std::mt19937* generator; //NOT OWNED BY ANGULAR DISTRIBUTION + std::uniform_real_distribution uniform_cosine_dist; + std::uniform_real_distribution uniform_prob_dist; + + double branchingRatio; + int L; + std::vector constants; + bool isoFlag; + }; - std::mt19937* generator; //NOT OWNED BY ANGULAR DISTRIBUTION - std::uniform_real_distribution uniform_cosine_dist; - std::uniform_real_distribution uniform_prob_dist; - - double branchingRatio; - int L; - std::vector constants; - bool isoFlag; -}; +} #endif \ No newline at end of file diff --git a/include/Eloss_Tables.h b/include/Eloss_Tables.h index fc12245..67b5abe 100644 --- a/include/Eloss_Tables.h +++ b/include/Eloss_Tables.h @@ -1,130 +1,133 @@ #ifndef ELOSS_TABLES_H #define ELOSS_TABLES_H + +namespace Mask { -#define MAX_Z 93 //Maximum number of elements for which we have hydrogen coefficients + #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} -}; + /*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 index fc0283f..bb0bad9 100644 --- a/include/EnergyLoss.h +++ b/include/EnergyLoss.h @@ -19,39 +19,40 @@ Written by G.W. McCann Aug. 2020 #include #include "MassLookup.h" -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(std::vector& Zt, std::vector& At, 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 PI = 3.14159265358979323846; - static constexpr double H_RESTMASS = 938.27231; //MeV, for beta calc -}; +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(std::vector& Zt, std::vector& At, 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/LayeredTarget.h b/include/LayeredTarget.h index 84da673..16206e4 100644 --- a/include/LayeredTarget.h +++ b/include/LayeredTarget.h @@ -1,7 +1,7 @@ /* LayeredTarget.h -Functional unit for targets in the SPANCRedux environment. Contains a +Functional unit for targets in the Mask environment. Contains a set (read: vector) of Targets for use in reaction calculations. In this way handles situations such as carbon backed targets @@ -15,27 +15,30 @@ Written by G.W. McCann Aug. 2020 #include #include -#include #include "Target.h" -class LayeredTarget { - - public: - LayeredTarget(); - ~LayeredTarget(); - void AddLayer(std::vector& Z, std::vector& A, std::vector& stoich, double thickness); - double GetProjectileEnergyLoss(int zp, int ap, double startEnergy, int rxnLayer, double angle); - double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle); - double GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle); - int GetNumberOfLayers(); - int FindLayerContaining(int Z, int A); - void SetName(std::string& n); - Target& GetLayerInfo(int index); - std::string& GetName(); +namespace Mask { - private: - std::vector layers; - std::string name; -}; + class LayeredTarget { + + public: + LayeredTarget(); + ~LayeredTarget(); + void AddLayer(std::vector& Z, std::vector& A, std::vector& stoich, double thickness); + double GetProjectileEnergyLoss(int zp, int ap, double startEnergy, int rxnLayer, double angle); + double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle); + double GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, int rxnLayer, double angle); + int FindLayerContaining(int Z, int A); + inline int GetNumberOfLayers() { return layers.size(); } + inline void SetName(std::string& n) { name = n; } + inline const Target& GetLayerInfo(int index) { return layers[index]; } + inline const std::string& GetName() { return name; } + + private: + std::vector layers; + std::string name; + }; + +} #endif diff --git a/include/LegendrePoly.h b/include/LegendrePoly.h index 270d42e..6e9f04b 100644 --- a/include/LegendrePoly.h +++ b/include/LegendrePoly.h @@ -1,12 +1,16 @@ #ifndef LEGENDREPOLY_H #define LEGENDREPOLY_H -double P_l(int l, double x); -double P_l_ROOT(double* x, double* pars); -double Normed_P_l_sq(int l, double x); +namespace Mask { -double P_0(double x); -double P_1(double x); -double P_2(double x); + double P_l(int l, double x); + double P_l_ROOT(double* x, double* pars); + double Normed_P_l_sq(int l, double x); + + double P_0(double x); + double P_1(double x); + double P_2(double x); + +} #endif \ No newline at end of file diff --git a/include/MassLookup.h b/include/MassLookup.h index ba7bd52..0a6f255 100644 --- a/include/MassLookup.h +++ b/include/MassLookup.h @@ -12,36 +12,38 @@ Converted to true singleton to simplify usage -- Aug. 2021 GWM #ifndef MASS_LOOKUP_H #define MASS_LOOKUP_H -#include #include #include #include -#include -class MassLookup { +namespace Mask { - public: - ~MassLookup(); - double FindMass(int Z, int A); - std::string FindSymbol(int Z, int A); - static MassLookup* GetInstance() { - if(s_instance == nullptr) { - s_instance = new MassLookup(); - } - return s_instance; - } + class MassLookup { + public: + ~MassLookup(); + double FindMass(int Z, int A); + std::string FindSymbol(int Z, int A); + + static MassLookup* GetInstance() { + if(s_instance == nullptr) { + s_instance = new MassLookup(); + } + return s_instance; + } + + private: + MassLookup(); + std::unordered_map massTable; + std::unordered_map elementTable; + + static MassLookup* s_instance; + + //constants + static constexpr double u_to_mev = 931.4940954; + static constexpr double electron_mass = 0.000548579909; + + }; - private: - MassLookup(); - std::unordered_map massTable; - std::unordered_map elementTable; - - static MassLookup* s_instance; - - //constants - static constexpr double u_to_mev = 931.4940954; - static constexpr double electron_mass = 0.000548579909; - -}; +} #endif diff --git a/include/Reaction.h b/include/Reaction.h index 1c3b940..396c935 100644 --- a/include/Reaction.h +++ b/include/Reaction.h @@ -14,66 +14,65 @@ namespace Mask { -class Reaction { -public: - Reaction(); - Reaction(int zt, int at, int zp, int ap, int ze, int ae); - ~Reaction(); - bool Calculate(); - void SetNuclei(int zt, int at, int zp, int ap, int ze, int ae); - void SetNuclei(const Nucleus* nucs); - void SetBeamKE(double bke); - void SetEjectileThetaType(int type); + class Reaction { + public: + Reaction(); + Reaction(int zt, int at, int zp, int ap, int ze, int ae); + ~Reaction(); + bool Calculate(); + void SetNuclei(int zt, int at, int zp, int ap, int ze, int ae); + void SetNuclei(const Nucleus* nucs); + void SetBeamKE(double bke); + void SetEjectileThetaType(int type); + + inline void SetLayeredTarget(LayeredTarget* targ) { target = targ; }; + inline void SetPolarRxnAngle(double theta) { m_theta = theta; }; + inline void SetAzimRxnAngle(double phi) { m_phi = phi; }; + inline void SetExcitation(double ex) { m_ex = ex; }; + inline void SetTarget(const Nucleus& nuc) { reactants[0] = nuc; }; + inline void SetTarget(int z, int a) { reactants[0] = Nucleus(z, a); }; + inline void SetProjectile(const Nucleus& nuc) { reactants[1] = nuc; }; + inline void SetProjectile(int z, int a) { reactants[1] = Nucleus(z, a); }; + inline void SetEjectile(const Nucleus& nuc) { reactants[2] = nuc; }; + inline void SetEjectile(int z, int a) { reactants[2] = Nucleus(z, a); }; + inline void SetResidual(const Nucleus& nuc) { reactants[3] = nuc; }; + inline void SetResidual(int z, int a) { reactants[3] = Nucleus(z, a); }; + inline void SetRxnLayer(int layer) { rxnLayer = layer; }; + inline void TurnOffResidualEloss() { resid_elossFlag = false; }; + inline void TurnOnResidualEloss() { resid_elossFlag = true; }; + inline bool IsDecay() { return decayFlag; }; + inline const Nucleus* GetNuclei() const { return &(reactants[0]); }; + inline const Nucleus& GetProjectile() const { return reactants[1]; }; + inline const Nucleus& GetTarget() const { return reactants[0]; }; + inline const Nucleus& GetEjectile() const { return reactants[2]; }; + inline const Nucleus& GetResidual() const { return reactants[3]; }; + inline int GetRxnLayer() { return rxnLayer; }; + inline void ResetTarget() { reactants[0].SetVectorCartesian(0,0,0, reactants[0].GetGroundStateMass()); } + inline void ResetProjectile() { reactants[1].SetVectorCartesian(0,0,0, reactants[1].GetGroundStateMass()); } + inline void ResetEjectile() { reactants[2].SetVectorCartesian(0,0,0, reactants[2].GetGroundStateMass()); } + inline void ResetResidual() { reactants[3].SetVectorCartesian(0,0,0, reactants[3].GetGroundStateMass()); } + + private: + void CalculateDecay(); //target -> light_decay (eject) + heavy_decay(resid) + void CalculateReaction(); //target + project -> eject + resid + void CalculateReactionThetaLab(); + void CalculateReactionThetaCM(); + + Nucleus reactants[4]; //0=target, 1=projectile, 2=ejectile, 3=residual + LayeredTarget* target; //not owned by Reaction + + double m_bke, m_theta, m_phi, m_ex; + + int rxnLayer; + int m_eject_theta_type; + + bool decayFlag, nuc_initFlag, resid_elossFlag; + + static constexpr int lab = 0; + static constexpr int center_of_mass = 1; + + }; - /*Setters and getters*/ - inline void SetLayeredTarget(LayeredTarget* targ) { target = targ; }; - inline void SetPolarRxnAngle(double theta) { m_theta = theta; }; - inline void SetAzimRxnAngle(double phi) { m_phi = phi; }; - inline void SetExcitation(double ex) { m_ex = ex; }; - inline void SetTarget(const Nucleus& nuc) { reactants[0] = nuc; }; - inline void SetTarget(int z, int a) { reactants[0] = Nucleus(z, a); }; - inline void SetProjectile(const Nucleus& nuc) { reactants[1] = nuc; }; - inline void SetProjectile(int z, int a) { reactants[1] = Nucleus(z, a); }; - inline void SetEjectile(const Nucleus& nuc) { reactants[2] = nuc; }; - inline void SetEjectile(int z, int a) { reactants[2] = Nucleus(z, a); }; - inline void SetResidual(const Nucleus& nuc) { reactants[3] = nuc; }; - inline void SetResidual(int z, int a) { reactants[3] = Nucleus(z, a); }; - inline void SetRxnLayer(int layer) { rxnLayer = layer; }; - inline void TurnOffResidualEloss() { resid_elossFlag = false; }; - inline void TurnOnResidualEloss() { resid_elossFlag = true; }; - inline bool IsDecay() { return decayFlag; }; - inline const Nucleus* GetNuclei() const { return &(reactants[0]); }; - inline const Nucleus& GetProjectile() const { return reactants[1]; }; - inline const Nucleus& GetTarget() const { return reactants[0]; }; - inline const Nucleus& GetEjectile() const { return reactants[2]; }; - inline const Nucleus& GetResidual() const { return reactants[3]; }; - inline void ResetTarget() { reactants[0].SetVectorCartesian(0,0,0,0); }; - inline void ResetProjectile() { reactants[1].SetVectorCartesian(0,0,0,0); }; - inline void ResetEjectile() { reactants[2].SetVectorCartesian(0,0,0,0); }; - inline void ResetResidual() { reactants[3].SetVectorCartesian(0,0,0,0); }; - inline int GetRxnLayer() { return rxnLayer; }; - -private: - void CalculateReaction(); //target + project -> eject + resid - void CalculateReactionThetaLab(); - void CalculateReactionThetaCM(); - void CalculateDecay(); //target -> light_decay (eject) + heavy_decay(resid) - - Nucleus reactants[4]; //0=target, 1=projectile, 2=ejectile, 3=residual - LayeredTarget* target; //not owned by Reaction - - double m_bke, m_theta, m_phi, m_ex; - - int rxnLayer; - int m_eject_theta_type; - - bool decayFlag, nuc_initFlag, resid_elossFlag; - - static constexpr int lab = 0; - static constexpr int center_of_mass = 1; - -}; - -}; +} #endif \ No newline at end of file diff --git a/include/Rotation.h b/include/Rotation.h index 7c3aab3..9b8ba3c 100644 --- a/include/Rotation.h +++ b/include/Rotation.h @@ -10,69 +10,69 @@ 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); + 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]; }; - -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); + + 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]; }; - -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); + + 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]; }; - -private: - void GenerateMatrix(); - double m_angle; - double m_matrix[3][3]; -}; - -}; + +} #endif \ No newline at end of file diff --git a/include/SabreDetector.h b/include/SabreDetector.h index 434f7fa..13cb241 100644 --- a/include/SabreDetector.h +++ b/include/SabreDetector.h @@ -171,4 +171,5 @@ private: }; + #endif diff --git a/include/SabreEfficiency.h b/include/SabreEfficiency.h index b795acf..7f0c619 100644 --- a/include/SabreEfficiency.h +++ b/include/SabreEfficiency.h @@ -21,8 +21,8 @@ private: std::vector detectors; - Target deadlayer; - Target sabre_eloss; + Mask::Target deadlayer; + Mask::Target sabre_eloss; DeadChannelMap dmap; //Sabre constants diff --git a/include/Target.h b/include/Target.h index 84dcce3..93c977b 100644 --- a/include/Target.h +++ b/include/Target.h @@ -1,7 +1,7 @@ /* Target.h -A basic target unit for use in the SPANCRedux environment. A target +A basic target unit for use in the Mask environment. A target is defined as a single compound with elements Z,A of a given stoichiometry Holds an energy loss class @@ -15,32 +15,35 @@ Written by G.W. McCann Aug. 2020 #include #include +#include #include "EnergyLoss.h" -class Target { +namespace Mask { - public: - Target(double thick); - ~Target(); - void SetElements(std::vector& z, std::vector& a, std::vector& stoich); - bool ContainsElement(int z, int a); - double getEnergyLossTotal(int zp, int ap, double startEnergy, double angle); - double getEnergyLossHalf(int zp, int ap, double startEnergy, double angle); - double getReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle); - double getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double angle); - double& GetThickness(); - int GetNumberOfElements(); - int GetElementZ(int index); - int GetElementA(int index); - int GetElementStoich(int index); - + class Target { + + public: + Target(double thick); + ~Target(); + void SetElements(std::vector& z, std::vector& a, std::vector& stoich); + bool ContainsElement(int z, int a); + double getEnergyLossTotal(int zp, int ap, double startEnergy, double angle); + double getEnergyLossHalf(int zp, int ap, double startEnergy, double angle); + double getReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle); + double getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double angle); + inline const double& GetThickness() { return thickness; } + inline int GetNumberOfElements() { return Z.size(); } + inline int GetElementZ(int index) { return Z[index]; } + inline int GetElementA(int index) { return A[index]; } + inline int GetElementStoich(int index) { return Stoich[index]; } + + private: + EnergyLoss eloss; + double thickness; + std::vector Z, A, Stoich; + + }; - private: - EnergyLoss eloss; - double thickness; - std::vector Z, A, Stoich; - static constexpr double PI = 3.14159265358979323846; - -}; +} #endif diff --git a/include/Vec3.h b/include/Vec3.h index 9be16c4..dcd1ad8 100644 --- a/include/Vec3.h +++ b/include/Vec3.h @@ -11,51 +11,51 @@ 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; + 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]; + }; - 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 index 5364bdd..15ca5b3 100644 --- a/include/Vec4.h +++ b/include/Vec4.h @@ -11,59 +11,59 @@ 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());}; + 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]; - 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 \ No newline at end of file diff --git a/input.txt b/input.txt index bc9234c..9311917 100644 --- a/input.txt +++ b/input.txt @@ -1,5 +1,5 @@ ----------Data Information---------- -OutputFile: /data1/gwm17/mask_tests/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.mask +OutputFile: test/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.mask SaveTree: yes SavePlots: yes ----------Reaction Information---------- diff --git a/src/AngularDistribution.cpp b/src/AngularDistribution.cpp index 1065c71..1919270 100644 --- a/src/AngularDistribution.cpp +++ b/src/AngularDistribution.cpp @@ -4,101 +4,105 @@ #include #include "LegendrePoly.h" -AngularDistribution::AngularDistribution() : - generator(nullptr), uniform_cosine_dist(-1.0, 1.0), uniform_prob_dist(0.0, 1.0), branchingRatio(1.0), L(0), isoFlag(true) -{ -} +namespace Mask { -AngularDistribution::AngularDistribution(const std::string& file) : - generator(nullptr), branchingRatio(1.0), L(0), isoFlag(true) -{ - ReadDistributionFile(file); -} - -AngularDistribution::~AngularDistribution() {} - -void AngularDistribution::ReadDistributionFile(const std::string& file) { - - if(file == "none" || file == "") { - L=0; - branchingRatio=1.0; - constants.clear(); - constants.push_back(0.5); - isoFlag = true; - return; + AngularDistribution::AngularDistribution() : + generator(nullptr), uniform_cosine_dist(-1.0, 1.0), uniform_prob_dist(0.0, 1.0), branchingRatio(1.0), L(0), isoFlag(true) + { } - - std::ifstream input(file); - std::string junk; - int l; - double par; - - if(!input.is_open()) { - std::cerr<<"Unable to open distribution file. All values reset to default."<>junk>>l; - while(input>>junk) { - input>>par; - constants.push_back(par); + + AngularDistribution::~AngularDistribution() {} + + void AngularDistribution::ReadDistributionFile(const std::string& file) { + + if(file == "none" || file == "") { + L=0; + branchingRatio=1.0; + constants.clear(); + constants.push_back(0.5); + isoFlag = true; + return; + } + + std::ifstream input(file); + std::string junk; + int l; + double par; + + if(!input.is_open()) { + std::cerr<<"Unable to open distribution file. All values reset to default."<>junk>>l; + while(input>>junk) { + input>>par; + constants.push_back(par); + } + input.close(); + + if(constants.size() != ((unsigned int) l+1)) { + std::cerr<<"Unexpected number of constants for given angular momentum! Expected "< branchingRatio) return -10; - - do { - probability = 0.0; - costheta = uniform_cosine_dist(*generator); + + double AngularDistribution::GetRandomCosTheta() { + + if(!IsGeneratorSet()) { + std::cerr<<"Random number generator is not set in AngularDistribution! Returning default value of 0"< probability); + if(test > branchingRatio) return -10; + + do { + probability = 0.0; + costheta = uniform_cosine_dist(*generator); + test = uniform_prob_dist(*generator); + for(unsigned int i=0; i probability); + + return costheta; + } - return costheta; } \ No newline at end of file diff --git a/src/DecaySystem.cpp b/src/DecaySystem.cpp index eaa0ba0..56308fd 100644 --- a/src/DecaySystem.cpp +++ b/src/DecaySystem.cpp @@ -59,17 +59,14 @@ namespace Mask { LinkTarget(); } - if(step1.IsDecay()) { - double rxnTheta = std::acos(decay1dist.GetRandomCosTheta()); - double rxnPhi = (*m_phi1Range)(*generator); - step1.SetPolarRxnAngle(rxnTheta); - step1.SetAzimRxnAngle(rxnPhi); + double rxnTheta = std::acos(decay1dist.GetRandomCosTheta()); + double rxnPhi = (*m_phi1Range)(*generator); + step1.SetPolarRxnAngle(rxnTheta); + step1.SetAzimRxnAngle(rxnPhi); - step1.TurnOnResidualEloss(); - step1.Calculate(); - } else { - return; - } + step1.TurnOnResidualEloss(); + step1.Calculate(); + } } \ No newline at end of file diff --git a/src/Detectors/SabreDetector.cpp b/src/Detectors/SabreDetector.cpp index 4e9a81d..c98fd0b 100644 --- a/src/Detectors/SabreDetector.cpp +++ b/src/Detectors/SabreDetector.cpp @@ -49,6 +49,7 @@ */ + #include "SabreDetector.h" SabreDetector::SabreDetector() : diff --git a/src/EnergyLoss.cpp b/src/EnergyLoss.cpp index c332ecf..f6b8467 100644 --- a/src/EnergyLoss.cpp +++ b/src/EnergyLoss.cpp @@ -14,231 +14,230 @@ Written by G.W. McCann Aug. 2020 #include "KinematicsExceptions.h" #include -EnergyLoss::EnergyLoss() { - comp_denom = 0; - ZP = -1; +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(std::vector& Zt, std::vector& At, 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; iFindMass(ZP, AP)*MEV2U; + } + + double e_final = e_initial; + double x_traversed = 0; + double x_step = 0.25*thickness; //initial step in x + double e_step = GetTotalStoppingPower(e_final)*x_step/1000.0; //initial step in e + double e_threshold = 0.05*e_initial; + + int depth=0; + + + if(thickness == 0.0 || e_initial == 0.0) + return 0; + + bool go = true; + while(go) { + //If intial guess of step size is too large, shrink until in range + if(e_step/e_final > 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"); + } + + 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& Zt, std::vector& At, 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; iFindMass(ZP, AP)*MEV2U; - } - - double e_final = e_initial; - double x_traversed = 0; - double x_step = 0.25*thickness; //initial step in x - double e_step = GetTotalStoppingPower(e_final)*x_step/1000.0; //initial step in e - double e_threshold = 0.05*e_initial; - - int depth=0; - - - if(thickness == 0.0) return 0; - else if(e_initial == 0.0) return 0; - - bool go = true; - while(go) { - //If intial guess of step size is too large, shrink until in range - if(e_step/e_final > 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(depth > 20)std::cout<<"depth: "<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"); - } - - 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 -LayeredTarget::LayeredTarget() {} +namespace Mask { -LayeredTarget::~LayeredTarget() {} - -/*Add in a Target made of a compound defined by a set of Zs, As, Ss, and a thickness*/ -void LayeredTarget::AddLayer(std::vector& Z, std::vector& A, std::vector& stoich, double thickness) { - Target t(thickness); - t.SetElements(Z, A, stoich); - layers.push_back(t); -} - -/* - Here projectile refers to the incoming reactant particle (i.e. the beam) - Calculates energy loss assuming that the reaction occurs in the middle of the target - Note that the layer order can matter! -*/ -double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy, int rxnLayer, double angle) { - - if(rxnLayer < 0 || ((unsigned int) rxnLayer) > layers.size()) { - std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"< layers.size()) { - std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"< layers.size()) { - std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"<=rxnLayer; i--) { - if(i == rxnLayer) { - eloss += layers[i].getReverseEnergyLossHalf(ze, ae, newEnergy, angle); - newEnergy = startEnergy + eloss; - } else { - eloss += layers[i].getReverseEnergyLossTotal(ze, ae, newEnergy, angle); - newEnergy = startEnergy + eloss; - } - } - - return eloss; -} - -/*Getters and Setters*/ - -int LayeredTarget::GetNumberOfLayers() { - return layers.size(); -} - -int LayeredTarget::FindLayerContaining(int Z, int A) { - for(unsigned int i=0; i& Z, std::vector& A, std::vector& stoich, double thickness) { + Target t(thickness); + t.SetElements(Z, A, stoich); + layers.push_back(t); + } + + /* + Here projectile refers to the incoming reactant particle (i.e. the beam) + Calculates energy loss assuming that the reaction occurs in the middle of the target + Note that the layer order can matter! + */ + double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy, int rxnLayer, double angle) { + + if(rxnLayer < 0 || ((unsigned int) rxnLayer) > layers.size()) { + std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"< layers.size()) { + std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"< layers.size()) { + std::cerr<<"Reaction layer in eloss calculation is not in range! Returning 0"<=rxnLayer; i--) { + if(i == rxnLayer) { + eloss += layers[i].getReverseEnergyLossHalf(ze, ae, newEnergy, angle); + newEnergy = startEnergy + eloss; + } else { + eloss += layers[i].getReverseEnergyLossTotal(ze, ae, newEnergy, angle); + newEnergy = startEnergy + eloss; + } + } + + return eloss; + } + + int LayeredTarget::FindLayerContaining(int Z, int A) { + for(unsigned int i=0; i -double P_l(int l, double x) { - if(l == 0) { - return 1.0; - } else if (l == 1) { - return x; - } else { - return (2.0*l - 1.0)/l*x*P_l(l-1, x) - (l-1.0)/l*P_l(l-2, x); +namespace Mask { + + double P_l(int l, double x) { + if(l == 0) { + return 1.0; + } else if (l == 1) { + return x; + } else { + return (2.0*l - 1.0)/l*x*P_l(l-1, x) - (l-1.0)/l*P_l(l-2, x); + } + } + + double Normed_P_l_sq(int l, double x) { + return (2.0*l+1.0)/2.0*std::pow(P_l(l, x), 2.0); + } + + double P_0(double x) { + return 1.0; + } + + double P_1(double x) { + return x; + } + + double P_2(double x) { + return 0.5*(3.0*x*x -1.0); + } + + double P_l_ROOT(double* x, double* pars) { + return P_l(pars[0], x[0]); } -} -double Normed_P_l_sq(int l, double x) { - return (2.0*l+1.0)/2.0*std::pow(P_l(l, x), 2.0); -} - -double P_0(double x) { - return 1.0; -} - -double P_1(double x) { - return x; -} - -double P_2(double x) { - return 0.5*(3.0*x*x -1.0); -} - -double P_l_ROOT(double* x, double* pars) { - return P_l(pars[0], x[0]); } \ No newline at end of file diff --git a/src/MassLookup.cpp b/src/MassLookup.cpp index 14e3171..3e7a396 100644 --- a/src/MassLookup.cpp +++ b/src/MassLookup.cpp @@ -11,52 +11,51 @@ Written by G.W. McCann Aug. 2020 #include "MassLookup.h" #include "KinematicsExceptions.h" +namespace Mask { -/* - Read in AMDC mass file, preformated to remove excess info. Here assumes that by default - the file is in a local directory etc/ -*/ - -MassLookup* MassLookup::s_instance = nullptr; - -MassLookup::MassLookup() { - std::ifstream massfile("./etc/mass.txt"); - if(massfile.is_open()) { - std::string junk, A, element; - int Z; - double atomicMassBig, atomicMassSmall, isotopicMass; - getline(massfile,junk); - getline(massfile,junk); - while(massfile>>junk) { - massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall; - isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev; - std::string key = "("+std::to_string(Z)+","+A+")"; - massTable[key] = isotopicMass; - elementTable[Z] = element; + MassLookup* MassLookup::s_instance = nullptr; + + MassLookup::MassLookup() { + + std::ifstream massfile("etc/mass.txt"); + if(massfile.is_open()) { + std::string junk, A, element; + int Z; + double atomicMassBig, atomicMassSmall, isotopicMass; + getline(massfile,junk); + getline(massfile,junk); + while(massfile>>junk) { + massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall; + isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev; + std::string key = "("+std::to_string(Z)+","+A+")"; + massTable[key] = isotopicMass; + elementTable[Z] = element; + } + } else { + throw MassFileException(); } - } else { - throw MassFileException(); } -} + + MassLookup::~MassLookup() {} + + //Returns nuclear mass in MeV + double MassLookup::FindMass(int Z, int A) { + std::string key = "("+std::to_string(Z)+","+std::to_string(A)+")"; + auto data = massTable.find(key); + if(data == massTable.end()) + throw MassException(); + + return data->second; + } + + //returns element symbol + std::string MassLookup::FindSymbol(int Z, int A) { + auto data = elementTable.find(Z); + if(data == elementTable.end()) + throw MassException(); + + std::string fullsymbol = std::to_string(A) + data->second; + return fullsymbol; + } -MassLookup::~MassLookup() {} - -//Returns nuclear mass in MeV -double MassLookup::FindMass(int Z, int A) { - std::string key = "("+std::to_string(Z)+","+std::to_string(A)+")"; - auto data = massTable.find(key); - if(data == massTable.end()) { - throw MassException(); - } - return data->second; -} - -//returns element symbol -std::string MassLookup::FindSymbol(int Z, int A) { - auto data = elementTable.find(Z); - if(data == elementTable.end()) { - throw MassException(); - } - std::string fullsymbol = std::to_string(A) + data->second; - return fullsymbol; -} +} \ No newline at end of file diff --git a/src/Nucleus.cpp b/src/Nucleus.cpp index d5d621f..153431a 100644 --- a/src/Nucleus.cpp +++ b/src/Nucleus.cpp @@ -11,37 +11,37 @@ namespace Mask { -Nucleus::Nucleus () : - Vec4(), m_z(0), m_a(0), m_gs_mass(0), m_theta_cm(0), m_symbol(""), m_detectFlag(false) -{ -} - -Nucleus::Nucleus(int Z, int A) : - Vec4(), m_z(Z), m_a(A), m_theta_cm(0), m_detectFlag(false) -{ - m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A); - m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A); - SetVectorCartesian(0,0,0,m_gs_mass); //by defualt a nucleus has mass given by the g.s. -} - -Nucleus::Nucleus(int Z, int A, double px, double py, double pz, double E) : - Vec4(px, py, pz, E), m_z(Z), m_a(A) -{ - m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A); - m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A); -} - -Nucleus::~Nucleus() {} - -bool Nucleus::SetIsotope(int Z, int A) { - if(Z>A) return false; + Nucleus::Nucleus () : + Vec4(), m_z(0), m_a(0), m_gs_mass(0), m_theta_cm(0), m_symbol(""), m_detectFlag(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; -} + Nucleus::Nucleus(int Z, int A) : + Vec4(), m_z(Z), m_a(A), m_theta_cm(0), m_detectFlag(false) + { + m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A); + m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A); + SetVectorCartesian(0,0,0,m_gs_mass); //by defualt a nucleus has mass given by the g.s. + } + + Nucleus::Nucleus(int Z, int A, double px, double py, double pz, double E) : + Vec4(px, py, pz, E), m_z(Z), m_a(A) + { + m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A); + m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A); + } + + Nucleus::~Nucleus() {} + + bool Nucleus::SetIsotope(int Z, int A) { + if(Z>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 +} \ No newline at end of file diff --git a/src/OneStepSystem.cpp b/src/OneStepSystem.cpp index ad2e257..6d0b053 100644 --- a/src/OneStepSystem.cpp +++ b/src/OneStepSystem.cpp @@ -55,23 +55,19 @@ namespace Mask { LinkTarget(); } - if(!step1.IsDecay()) { - //Sample parameters - double bke = (*m_beamDist)(*generator); - double rxnTheta = std::acos((*m_theta1Range)(*generator)); - double rxnPhi = (*m_phi1Range)(*generator); - double residEx = (*m_exDist)(*generator); + //Sample parameters + double bke = (*m_beamDist)(*generator); + double rxnTheta = std::acos((*m_theta1Range)(*generator)); + double rxnPhi = (*m_phi1Range)(*generator); + double residEx = (*m_exDist)(*generator); - step1.SetBeamKE(bke); - step1.SetPolarRxnAngle(rxnTheta); - step1.SetAzimRxnAngle(rxnPhi); - step1.SetExcitation(residEx); + step1.SetBeamKE(bke); + step1.SetPolarRxnAngle(rxnTheta); + step1.SetAzimRxnAngle(rxnPhi); + step1.SetExcitation(residEx); - step1.TurnOnResidualEloss(); - step1.Calculate(); - } else { - return; - } + step1.TurnOnResidualEloss(); + step1.Calculate(); } } \ No newline at end of file diff --git a/src/Reaction.cpp b/src/Reaction.cpp index 3aac0cf..bf9e26e 100644 --- a/src/Reaction.cpp +++ b/src/Reaction.cpp @@ -11,234 +11,236 @@ namespace Mask { -Reaction::Reaction() : - target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), nuc_initFlag(false), resid_elossFlag(false) -{ -} - -Reaction::Reaction(int zt, int at, int zp, int ap, int ze, int ae) : - target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), resid_elossFlag(false) -{ - SetNuclei(zt, at, zp, ap, ze, ae); -} - -Reaction::~Reaction() -{ -} - -bool Reaction::Calculate() { - - if(!nuc_initFlag) return false; - - if(decayFlag) { - CalculateDecay(); - return true; - } else { - CalculateReaction(); - return true; - } -} - -//Deep copy of nucleus array -void Reaction::SetNuclei(const Nucleus* nucs) { - reactants[0] = nucs[0]; - reactants[1] = nucs[1]; - reactants[2] = nucs[2]; - reactants[3] = nucs[3]; - nuc_initFlag = true; -} - -void Reaction::SetNuclei(int zt, int at, int zp, int ap, int ze, int ae) { - int zr, ar; - reactants[0] = Nucleus(zt, at); - reactants[2] = Nucleus(ze, ae); - if(ap == 0) { - decayFlag = true; - zr = zt - ze; - ar = at - ae; - } else { - reactants[1] = Nucleus(zp, ap); - decayFlag = false; - zr = zt + zp - ze; - ar = at + ap - ae; - } - - if(zr < 0 || ar <= 0) { - nuc_initFlag = false; - } else { - reactants[3] = Nucleus(zr, ar); - nuc_initFlag = true; - } -} - -void Reaction::SetBeamKE(double bke) { - if(!nuc_initFlag) return; - else if(decayFlag) return; - m_bke = bke - target->GetProjectileEnergyLoss(reactants[1].GetZ(), reactants[1].GetA(), bke, rxnLayer, 0); -}; - -void Reaction::SetEjectileThetaType(int type) { - if(decayFlag) return; - if(type != center_of_mass && type != lab) return; - - m_eject_theta_type = type; -} - -//Methods given by Iliadis in Nuclear Physics of Stars, Appendix C -//For use with lab frame restricted angles. May not give appropriate disribution for ejectile -void Reaction::CalculateReactionThetaLab() { - reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass()); - double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass())); - double beam_E = m_bke + reactants[1].GetGroundStateMass(); - reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E); - - double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex); - - double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()); - if(m_bke < Ethresh) { - throw EnergyThresholdException(); - } - - double term1 = sqrt(reactants[1].GetGroundStateMass()*reactants[2].GetGroundStateMass()*m_bke)/(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())*cos(m_theta); - double term2 = (m_bke*(reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()) + reactants[3].GetGroundStateMass()*Q)/(reactants[3].GetGroundStateMass() + reactants[2].GetGroundStateMass()); - double sqrt_pos_ejectKE = term1 + sqrt(term1*term1 + term2); - double sqrt_neg_ejectKE = term1 - sqrt(term1*term1 + term2); - double ejectKE; - if(sqrt_pos_ejectKE > 0) { - ejectKE = sqrt_pos_ejectKE*sqrt_pos_ejectKE; - } else { - ejectKE = sqrt_neg_ejectKE*sqrt_neg_ejectKE; - } - double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass())); - double ejectE = ejectKE + reactants[2].GetGroundStateMass(); - - reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP, ejectE); - - reactants[3] = reactants[0] + reactants[1] - reactants[2]; - - ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, m_theta); - ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass())); - ejectE = ejectKE + reactants[2].GetGroundStateMass(); - reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP, ejectE); - - if(resid_elossFlag) { - double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta()); - double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass())); - double residE = residKE + reactants[3].GetInvMass(); - reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE); - } -} - -//Methods from original ANASEN. Gives proper distribution for inverse kinematics. -void Reaction::CalculateReactionThetaCM() { - //Target assumed at rest, with 0 excitation energy - reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass()); - double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass())); - double beam_E = m_bke + reactants[1].GetGroundStateMass(); - reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E); - - - double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex); - - double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()); - if(m_bke < Ethresh) { - throw EnergyThresholdException(); + Reaction::Reaction() : + target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), nuc_initFlag(false), resid_elossFlag(false) + { } - auto parent = reactants[0] + reactants[1]; - double boost2lab[3]; - double boost2cm[3]; - const double* boost = parent.GetBoost(); - for(int i=0; i<3; i++) { - boost2lab[i] = boost[i]; - boost2cm[i] = boost[i]*(-1.0); + Reaction::Reaction(int zt, int at, int zp, int ap, int ze, int ae) : + target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), resid_elossFlag(false) + { + SetNuclei(zt, at, zp, ap, ze, ae); } - parent.ApplyBoost(boost2cm); - double ejectE_cm = (std::pow(reactants[2].GetGroundStateMass(), 2.0) - std::pow(reactants[3].GetGroundStateMass() + m_ex, 2.0) + std::pow(parent.GetE(),2.0))/(2.0*parent.GetE()); - double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - std::pow(reactants[2].GetGroundStateMass(), 2.0)); - reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP_cm, ejectE_cm); - reactants[2].ApplyBoost(boost2lab); - reactants[3] = reactants[0] + reactants[1] - reactants[2]; - - double ejectKE = reactants[2].GetKE(); - double ejectP = reactants[2].GetP(); - double ejectE = reactants[2].GetE(); - //energy loss for ejectile (after reaction!) - ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, reactants[2].GetTheta()); - ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass())); - ejectE = ejectKE + reactants[2].GetGroundStateMass(); - reactants[2].SetVectorSpherical(reactants[2].GetTheta(), reactants[2].GetPhi(), ejectP, ejectE); - - //if on, get eloss for residual (after reaction!) - if(resid_elossFlag) { - double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta()); - double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass())); - double residE = residKE + reactants[3].GetInvMass(); - reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE); + + Reaction::~Reaction() + { } -} - -void Reaction::CalculateReaction() { - switch(m_eject_theta_type) { - case center_of_mass: - { - CalculateReactionThetaCM(); - break; - } - case lab: - { - CalculateReactionThetaLab(); - break; + + bool Reaction::Calculate() { + + if(!nuc_initFlag) + return false; + + if(decayFlag) { + CalculateDecay(); + return true; + } else { + CalculateReaction(); + return true; } } + + //Deep copy of nucleus array + void Reaction::SetNuclei(const Nucleus* nucs) { + reactants[0] = nucs[0]; + reactants[1] = nucs[1]; + reactants[2] = nucs[2]; + reactants[3] = nucs[3]; + nuc_initFlag = true; + } + + void Reaction::SetNuclei(int zt, int at, int zp, int ap, int ze, int ae) { + int zr, ar; + reactants[0] = Nucleus(zt, at); + reactants[2] = Nucleus(ze, ae); + if(ap == 0) { + decayFlag = true; + zr = zt - ze; + ar = at - ae; + } else { + reactants[1] = Nucleus(zp, ap); + decayFlag = false; + zr = zt + zp - ze; + ar = at + ap - ae; + } + + if(zr < 0 || ar <= 0) { + nuc_initFlag = false; + } else { + reactants[3] = Nucleus(zr, ar); + nuc_initFlag = true; + } + } + + void Reaction::SetBeamKE(double bke) { + if(!nuc_initFlag || decayFlag) + return; + + m_bke = bke - target->GetProjectileEnergyLoss(reactants[1].GetZ(), reactants[1].GetA(), bke, rxnLayer, 0); + } + + void Reaction::SetEjectileThetaType(int type) { + if(decayFlag) return; + if(type != center_of_mass && type != lab) return; + + m_eject_theta_type = type; + } + + //Methods given by Iliadis in Nuclear Physics of Stars, Appendix C + //For use with lab frame restricted angles. May not give appropriate disribution for ejectile + void Reaction::CalculateReactionThetaLab() { + reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass()); + double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass())); + double beam_E = m_bke + reactants[1].GetGroundStateMass(); + reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E); + + double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex); + + double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()); + if(m_bke < Ethresh) { + throw EnergyThresholdException(); + } + + double term1 = sqrt(reactants[1].GetGroundStateMass()*reactants[2].GetGroundStateMass()*m_bke)/(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())*cos(m_theta); + double term2 = (m_bke*(reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()) + reactants[3].GetGroundStateMass()*Q)/(reactants[3].GetGroundStateMass() + reactants[2].GetGroundStateMass()); + double sqrt_pos_ejectKE = term1 + sqrt(term1*term1 + term2); + double sqrt_neg_ejectKE = term1 - sqrt(term1*term1 + term2); + double ejectKE; + if(sqrt_pos_ejectKE > 0) { + ejectKE = sqrt_pos_ejectKE*sqrt_pos_ejectKE; + } else { + ejectKE = sqrt_neg_ejectKE*sqrt_neg_ejectKE; + } + double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass())); + double ejectE = ejectKE + reactants[2].GetGroundStateMass(); + + reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP, ejectE); + + reactants[3] = reactants[0] + reactants[1] - reactants[2]; + + ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, m_theta); + ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass())); + ejectE = ejectKE + reactants[2].GetGroundStateMass(); + reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP, ejectE); + + if(resid_elossFlag) { + double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta()); + double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass())); + double residE = residKE + reactants[3].GetInvMass(); + reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE); + } + } + + //Methods from original ANASEN. Gives proper distribution for inverse kinematics. + void Reaction::CalculateReactionThetaCM() { + //Target assumed at rest, with 0 excitation energy + reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass()); + double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass())); + double beam_E = m_bke + reactants[1].GetGroundStateMass(); + reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E); + + + double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex); + + double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()); + if(m_bke < Ethresh) { + throw EnergyThresholdException(); + } + + auto parent = reactants[0] + reactants[1]; + double boost2lab[3]; + double boost2cm[3]; + const double* boost = parent.GetBoost(); + for(int i=0; i<3; i++) { + boost2lab[i] = boost[i]; + boost2cm[i] = boost[i]*(-1.0); + } + parent.ApplyBoost(boost2cm); + double ejectE_cm = (std::pow(reactants[2].GetGroundStateMass(), 2.0) - std::pow(reactants[3].GetGroundStateMass() + m_ex, 2.0) + std::pow(parent.GetE(),2.0))/(2.0*parent.GetE()); + double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - std::pow(reactants[2].GetGroundStateMass(), 2.0)); + reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP_cm, ejectE_cm); + reactants[2].ApplyBoost(boost2lab); + reactants[3] = reactants[0] + reactants[1] - reactants[2]; + + double ejectKE = reactants[2].GetKE(); + double ejectP = reactants[2].GetP(); + double ejectE = reactants[2].GetE(); + //energy loss for ejectile (after reaction!) + ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, reactants[2].GetTheta()); + ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass())); + ejectE = ejectKE + reactants[2].GetGroundStateMass(); + reactants[2].SetVectorSpherical(reactants[2].GetTheta(), reactants[2].GetPhi(), ejectP, ejectE); + + //if on, get eloss for residual (after reaction!) + if(resid_elossFlag) { + double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta()); + double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass())); + double residE = residKE + reactants[3].GetInvMass(); + reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE); + } + } + + void Reaction::CalculateReaction() { + switch(m_eject_theta_type) { + case center_of_mass: + { + CalculateReactionThetaCM(); + break; + } + case lab: + { + CalculateReactionThetaLab(); + break; + } + } + } + + //Calculate in CM, where decay is isotropic + void Reaction::CalculateDecay() { + + double Q = reactants[0].GetInvMass() - reactants[2].GetGroundStateMass() - reactants[3].GetGroundStateMass(); + if(Q < 0) { + throw QValueException(); + } + + const double* boost = reactants[0].GetBoost(); + double boost2cm[3]; + double boost2lab[3]; + for(int i=0; i<3; i++) { + boost2lab[i] = boost[i]; + boost2cm[i] = boost[i]*(-1.0); + } + + reactants[0].ApplyBoost(&(boost2cm[0])); + double ejectE_cm = (reactants[2].GetGroundStateMass()*reactants[2].GetGroundStateMass() - reactants[3].GetGroundStateMass()*reactants[3].GetGroundStateMass() + reactants[0].GetE()*reactants[0].GetE())/ + (2.0*reactants[0].GetE()); + double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - reactants[2].GetGroundStateMass()*reactants[2].GetGroundStateMass()); + + reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP_cm, ejectE_cm); + reactants[2].SetThetaCM(m_theta); + + reactants[0].ApplyBoost(boost2lab); + reactants[2].ApplyBoost(boost2lab); + + reactants[3] = reactants[0] - reactants[2]; + + //energy loss for the *light* break up nucleus + double ejectKE = reactants[2].GetKE() - target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), reactants[2].GetKE(), rxnLayer, reactants[2].GetTheta()); + double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0*reactants[2].GetGroundStateMass())); + double ejectE = ejectKE + reactants[2].GetGroundStateMass(); + reactants[2].SetVectorSpherical(reactants[2].GetTheta(), reactants[2].GetPhi(), ejectP, ejectE); + + //if on, get eloss for *heavy* break up nucleus + if(resid_elossFlag) { + + double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta()); + double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass())); + double residE = residKE + reactants[3].GetInvMass(); + reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE); + } + } + } -//Calculate in CM, where decay is isotropic -void Reaction::CalculateDecay() { - - double Q = reactants[0].GetInvMass() - reactants[2].GetGroundStateMass() - reactants[3].GetGroundStateMass(); - if(Q < 0) { - throw QValueException(); - } - - const double* boost = reactants[0].GetBoost(); - double boost2cm[3]; - double boost2lab[3]; - for(int i=0; i<3; i++) { - boost2lab[i] = boost[i]; - boost2cm[i] = boost[i]*(-1.0); - } - - reactants[0].ApplyBoost(&(boost2cm[0])); - double ejectE_cm = (reactants[2].GetGroundStateMass()*reactants[2].GetGroundStateMass() - reactants[3].GetGroundStateMass()*reactants[3].GetGroundStateMass() + reactants[0].GetE()*reactants[0].GetE())/ - (2.0*reactants[0].GetE()); - double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - reactants[2].GetGroundStateMass()*reactants[2].GetGroundStateMass()); - - reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP_cm, ejectE_cm); - reactants[2].SetThetaCM(m_theta); - - reactants[0].ApplyBoost(boost2lab); - reactants[2].ApplyBoost(boost2lab); - - reactants[3] = reactants[0] - reactants[2]; - - //energy loss for the *light* break up nucleus - double ejectKE = reactants[2].GetKE() - target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), reactants[2].GetKE(), rxnLayer, reactants[2].GetTheta()); - double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0*reactants[2].GetGroundStateMass())); - double ejectE = ejectKE + reactants[2].GetGroundStateMass(); - reactants[2].SetVectorSpherical(reactants[2].GetTheta(), reactants[2].GetPhi(), ejectP, ejectE); - - //if on, get eloss for *heavy* break up nucleus - if(resid_elossFlag) { - - double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta()); - double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass())); - double residE = residKE + reactants[3].GetInvMass(); - reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE); - } -} - -}; - diff --git a/src/Rotation.cpp b/src/Rotation.cpp index 8b40a30..aa07ad8 100644 --- a/src/Rotation.cpp +++ b/src/Rotation.cpp @@ -7,65 +7,65 @@ 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; -} + 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/Target.cpp b/src/Target.cpp index 24b5b35..83e3599 100644 --- a/src/Target.cpp +++ b/src/Target.cpp @@ -12,77 +12,70 @@ Written by G.W. McCann Aug. 2020 */ #include "Target.h" -/*Targets must be of known thickness*/ -Target::Target(double thick) { - thickness = thick; -} +namespace Mask { -Target::~Target() { -} + /*Targets must be of known thickness*/ + Target::Target(double thick) { + thickness = thick; + } + + Target::~Target() {} + + /*Set target elements of given Z, A, S*/ + void Target::SetElements(std::vector& z, std::vector& a, std::vector& stoich) { + Z = z; + A = a; + Stoich = stoich; + + eloss.SetTargetComponents(Z, A, Stoich); + } + + /*Element verification*/ + bool Target::ContainsElement(int z, int a) { + for(unsigned int i=0; i M_PI/2.) + theta = M_PI - theta; -/*Set target elements of given Z, A, S*/ -void Target::SetElements(std::vector& z, std::vector& a, std::vector& stoich) { - Z = z; - A = a; - Stoich = stoich; - - eloss.SetTargetComponents(Z, A, Stoich); -} + return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta))); + } + + /*Calculates energy loss for travelling halfway through the target*/ + double Target::getEnergyLossHalf(int zp, int ap, double startEnergy, double theta) { + if(theta == M_PI/2.) + return startEnergy; + else if (theta > M_PI/2.) + theta = M_PI - theta; -/*Element verification*/ -bool Target::ContainsElement(int z, int a) { - for(unsigned int i=0; i M_PI/2.) + theta = M_PI - theta; -/*Calculates energy loss for travelling all the way through the target*/ -double Target::getEnergyLossTotal(int zp, int ap, double startEnergy, double theta) { - if(theta == PI/2.) return startEnergy; - else if (theta > PI/2.) theta = PI - theta; - return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta))); -} + return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta))); + } + + /*Calculates reverse energy loss for travelling half way through the target*/ + double Target::getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double theta) { + if(theta == M_PI/2.) + return finalEnergy; + else if (theta > M_PI/2.) + theta = M_PI - theta; -/*Calculates energy loss for travelling halfway through the target*/ -double Target::getEnergyLossHalf(int zp, int ap, double startEnergy, double theta) { - if(theta == PI/2.) return startEnergy; - else if (theta > PI/2.) theta = PI - theta; - return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/(2.0*fabs(cos(theta)))); -} + return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta)))); + } -/*Calculates reverse energy loss for travelling all the way through the target*/ -double Target::getReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double theta) { - if(theta == PI/2.) return finalEnergy; - else if (theta > PI/2.) theta = PI - theta; - return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta))); -} - -/*Calculates reverse energy loss for travelling half way through the target*/ -double Target::getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double theta) { - if(theta == PI/2.) return finalEnergy; - else if (theta > PI/2.) theta = PI - theta; - return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta)))); -} - -/*Getter functions*/ - -double& Target::GetThickness() { - return thickness; -} - -int Target::GetNumberOfElements() { - return Z.size(); -} - -int Target::GetElementZ(int index) { - return Z[index]; -} - -int Target::GetElementA(int index) { - return A[index]; -} - -int Target::GetElementStoich(int index) { - return Stoich[index]; } diff --git a/src/Vec3.cpp b/src/Vec3.cpp index dde73be..42abdc9 100644 --- a/src/Vec3.cpp +++ b/src/Vec3.cpp @@ -8,41 +8,41 @@ namespace Mask { -Vec3::Vec3() { - m_data[0] = 0.; - m_data[1] = 0.; - m_data[2] = 0.; + 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); + } + } - -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/Vec4.cpp b/src/Vec4.cpp index 8115349..493ea3e 100644 --- a/src/Vec4.cpp +++ b/src/Vec4.cpp @@ -10,64 +10,64 @@ namespace Mask { -Vec4::Vec4() { - for(auto& val: m_data) - val = 0.0; - for(auto& val: m_boost) - val = 0.0; -} + 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(); + } -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 +} \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index 3fd5446..0574e2c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -17,6 +17,7 @@ int main(int argc, char** argv) { sw.Start(); try { if(!calculator.LoadConfig(argv[1])) { + std::cerr<<"Unable to read input file!"<