Skip to content

Commit

Permalink
Merge pull request #56 from NESTCollaboration/test
Browse files Browse the repository at this point in the history 
Sync with NEST v2.2.1
  • Loading branch information
@sophiafarrell
sophiafarrell authored Mar 1, 2021
2 parents 788f279 + e455835 commit 11b1710
Show file tree
Hide file tree
Showing 16 changed files with 477 additions and 422 deletions.
4 changes: 2 additions & 2 deletions MANIFEST.in
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
recursive-include src *.cpp
recursive-include src *.hh
recursive-include include *.h
recursive-include src *.py
include *.md
include *.txt
global-include *
include *.txt
314 changes: 157 additions & 157 deletions src/nestpy/GammaHandler.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,157 +1,157 @@
#include <exception>
#include "GammaHandler.hh"
#include <stdexcept>
#include "TestSpectra.hh"

using namespace std;

double yMax = 1.0; //arbitrary y max, might need to change
double brThresh = 0.1;

const vector<vector<double>>& GammaHandler::sourceLookupTable(const std::string& source)
{
// energy container vector orginized as {energy, branch ratio, PE mass attenuation coef, Compton coef, Pair Production
// coef}
typedef vector<vector<double>> LookupTable;
static const LookupTable co57Info {
{ 122.0, 0.856, 1.793, 0.1081, 0.00 },
{ 136.0, 0.1068, 0.5651, 0.1019, 0.00 },
{ 14.0, 0.0916, 55.5, 0.0744, 0.00 },
};
static const LookupTable co60Info {
{ 1332.0, 0.9998, 0.001991, 0.04244, 0.0008853 },
{ 1173.0, 0.9985, 0.004126, 0.05156, 0.00 },
};
static const LookupTable cs137Info {
{ 662.0, 0.851, 0.01338, 0.06559, 0.00 },
{ 284.0, 0.0006, 0.08009, 0.08495, 0.00 },
};
if (source == "Co57")
{
return co57Info;
}
else if (source == "Co60")
{
return co60Info;
}
else if (source == "Cs137")
{
return cs137Info;
}
cerr << source << " Is not a valid option!" << endl;
throw std::invalid_argument { source + " is not a valid option!" };
return co57Info;
}

const vector<double> GammaHandler::combineSpectra(double emin, double emax, string source) {
double brSum = 0.0;
double fValue = 0.0;
double pe = 0.0;
double compton = 0.0;
double pp = 0.0;

vector<vector<double>> sourceInfo = GammaHandler::sourceLookupTable(source);

vector<double> xyTry = {
emin + (emax - emin) * RandomGen::rndm()->rand_uniform(),
yMax * RandomGen::rndm()->rand_uniform(), 1.};

while(xyTry[2] > 0.) {
pe = GammaHandler::photoIonization(sourceInfo, xyTry);
compton = GammaHandler::compton(sourceInfo, xyTry);
pp = GammaHandler::pairProduction(sourceInfo, xyTry);
fValue = pe + compton + pp;
xyTry = RandomGen::rndm()->VonNeumann(emin, emax, 0., yMax, xyTry[0],
xyTry[1], fValue);
}

vector<double> keV_vec = {-1, -1, -1};

if( pe > 0.0) {
keV_vec[0] = xyTry[0];
}
if(compton > 0.0) {
keV_vec[1] = xyTry[0];
}
if(pp > 0.0) {
keV_vec[2] = xyTry[0];
}
return keV_vec;
}

double GammaHandler::photoIonization(const vector<vector<double>>& sourceInfo, const vector<double>& xyTry)
{
// implement simple delta function to the spectrum
std::size_t index { 0 };
bool found = false;
for (int i = 0; i < sourceInfo.size(); i++)
{
double initialEnergy = sourceInfo[i][0];
if (abs(xyTry[0] - initialEnergy) < brThresh)
{
index = i;
found = true;
}
}
double br = sourceInfo[index][1];
double pe = sourceInfo[index][2];
double co = sourceInfo[index][3];
double pp = sourceInfo[index][4];
if(found) return yMax * br * (pe / (pe + co + pp));
return 0.0;
}

double GammaHandler::compton(const vector<vector<double>>& sourceInfo, const vector<double>& xyTry) {
double energyScaleFactor = ElectronRestMassEnergy; //mc^2 for electron mass in keV
double thetaMin = 0.0;
double thetaMax = M_PI;
int simpIterations = 100;
double simpStep = (thetaMax - thetaMin)/simpIterations;
double simpCurrentStep = thetaMin;
double shiftedEnergy, simpResult, kn, B, rY, rPsi, initialEnergy;
bool draw = true;
double a = 1.0/137.04;
double re = pow(0.38616, -12);

//loop over gamma energies
for(int i = 0; i < sourceInfo.size(); i++) {
double initialEnergy = sourceInfo[i][0];
double br = sourceInfo[i][1];
double pe = sourceInfo[i][2];
double co = sourceInfo[i][3];
double pp = sourceInfo[i][4];
//get shifted energy with MC
bool draw = true;
while(draw){
rPsi = M_PI * RandomGen::rndm()->rand_uniform();
rY = 10* RandomGen::rndm()->rand_uniform();

B = 1.0/(1.0+initialEnergy/energyScaleFactor*(1-cos(rPsi)));
kn = M_PI*pow(B,2)*(B+1.0/B-pow(sin(rPsi),2))*sin(rPsi); //Klein-Nishina
if(rY<kn) draw = false;
}
shiftedEnergy = initialEnergy * (1.0-1.0/(1.0+initialEnergy/energyScaleFactor*(1.0-cos(rPsi)))); //shifted energy formula
if(abs(xyTry[0]-shiftedEnergy) < brThresh) {
return kn*yMax*br*(co/(pe+co+pp));
}
}
return 0.0;
}

double GammaHandler::pairProduction(const vector<vector<double>>& sourceInfo, const vector<double>& xyTry) {
double energyScaleFactor = ElectronRestMassEnergy; //mc^2 for electron mass in keV
double initialEnergy, shiftedEnergy;
//loop over allowed gamma energies
for(int i = 0; i < sourceInfo.size(); i++) {
double initialEnergy = sourceInfo[i][0];
double br = sourceInfo[i][1];
double pe = sourceInfo[i][2];
double co = sourceInfo[i][3];
double pp = sourceInfo[i][4];
shiftedEnergy = (0.5)*(initialEnergy - 2*energyScaleFactor); //Pair production energy
if(abs(xyTry[0]-shiftedEnergy) < brThresh) {
return yMax*br*(pp/(pe+co+pp));
}
}
return 0.0;
}
#include <exception>
#include "GammaHandler.hh"
#include <stdexcept>
#include "TestSpectra.hh"

using namespace std;

double yMax = 1.0; //arbitrary y max, might need to change
double brThresh = 0.1;

const vector<vector<double>>& GammaHandler::sourceLookupTable(const std::string& source)
{
// energy container vector orginized as {energy, branch ratio, PE mass attenuation coef, Compton coef, Pair Production
// coef}
typedef vector<vector<double>> LookupTable;
static const LookupTable co57Info {
{ 122.0, 0.856, 1.793, 0.1081, 0.00 },
{ 136.0, 0.1068, 0.5651, 0.1019, 0.00 },
{ 14.0, 0.0916, 55.5, 0.0744, 0.00 },
};
static const LookupTable co60Info {
{ 1332.0, 0.9998, 0.001991, 0.04244, 0.0008853 },
{ 1173.0, 0.9985, 0.004126, 0.05156, 0.00 },
};
static const LookupTable cs137Info {
{ 662.0, 0.851, 0.01338, 0.06559, 0.00 },
{ 284.0, 0.0006, 0.08009, 0.08495, 0.00 },
};
if (source == "Co57")
{
return co57Info;
}
else if (source == "Co60")
{
return co60Info;
}
else if (source == "Cs137")
{
return cs137Info;
}
cerr << source << " Is not a valid option!" << endl;
throw std::invalid_argument { source + " is not a valid option!" };
return co57Info;
}

const vector<double> GammaHandler::combineSpectra(double emin, double emax, string source) {
double brSum = 0.0;
double fValue = 0.0;
double pe = 0.0;
double compton = 0.0;
double pp = 0.0;

vector<vector<double>> sourceInfo = GammaHandler::sourceLookupTable(source);

vector<double> xyTry = {
emin + (emax - emin) * RandomGen::rndm()->rand_uniform(),
yMax * RandomGen::rndm()->rand_uniform(), 1.};

while(xyTry[2] > 0.) {
pe = GammaHandler::photoIonization(sourceInfo, xyTry);
compton = GammaHandler::compton(sourceInfo, xyTry);
pp = GammaHandler::pairProduction(sourceInfo, xyTry);
fValue = pe + compton + pp;
xyTry = RandomGen::rndm()->VonNeumann(emin, emax, 0., yMax, xyTry[0],
xyTry[1], fValue);
}

vector<double> keV_vec = {-1, -1, -1};

if( pe > 0.0) {
keV_vec[0] = xyTry[0];
}
if(compton > 0.0) {
keV_vec[1] = xyTry[0];
}
if(pp > 0.0) {
keV_vec[2] = xyTry[0];
}
return keV_vec;
}

double GammaHandler::photoIonization(const vector<vector<double>>& sourceInfo, const vector<double>& xyTry)
{
// implement simple delta function to the spectrum
std::size_t index { 0 };
bool found = false;
for (int i = 0; i < sourceInfo.size(); i++)
{
double initialEnergy = sourceInfo[i][0];
if (abs(xyTry[0] - initialEnergy) < brThresh)
{
index = i;
found = true;
}
}
double br = sourceInfo[index][1];
double pe = sourceInfo[index][2];
double co = sourceInfo[index][3];
double pp = sourceInfo[index][4];
if(found) return yMax * br * (pe / (pe + co + pp));
return 0.0;
}

double GammaHandler::compton(const vector<vector<double>>& sourceInfo, const vector<double>& xyTry) {
double energyScaleFactor = ElectronRestMassEnergy; //mc^2 for electron mass in keV
double thetaMin = 0.0;
double thetaMax = M_PI;
int simpIterations = 100;
double simpStep = (thetaMax - thetaMin)/simpIterations;
double simpCurrentStep = thetaMin;
double shiftedEnergy, simpResult, kn, B, rY, rPsi, initialEnergy;
bool draw = true;
double a = 1.0/137.04;
double re = pow(0.38616, -12);

//loop over gamma energies
for(int i = 0; i < sourceInfo.size(); i++) {
double initialEnergy = sourceInfo[i][0];
double br = sourceInfo[i][1];
double pe = sourceInfo[i][2];
double co = sourceInfo[i][3];
double pp = sourceInfo[i][4];
//get shifted energy with MC
bool draw = true;
while(draw){
rPsi = M_PI * RandomGen::rndm()->rand_uniform();
rY = 10* RandomGen::rndm()->rand_uniform();

B = 1.0/(1.0+initialEnergy/energyScaleFactor*(1-cos(rPsi)));
kn = M_PI*pow(B,2)*(B+1.0/B-pow(sin(rPsi),2))*sin(rPsi); //Klein-Nishina
if(rY<kn) draw = false;
}
shiftedEnergy = initialEnergy * (1.0-1.0/(1.0+initialEnergy/energyScaleFactor*(1.0-cos(rPsi)))); //shifted energy formula
if(abs(xyTry[0]-shiftedEnergy) < brThresh) {
return kn*yMax*br*(co/(pe+co+pp));
}
}
return 0.0;
}

double GammaHandler::pairProduction(const vector<vector<double>>& sourceInfo, const vector<double>& xyTry) {
double energyScaleFactor = ElectronRestMassEnergy; //mc^2 for electron mass in keV
double initialEnergy, shiftedEnergy;
//loop over allowed gamma energies
for(int i = 0; i < sourceInfo.size(); i++) {
double initialEnergy = sourceInfo[i][0];
double br = sourceInfo[i][1];
double pe = sourceInfo[i][2];
double co = sourceInfo[i][3];
double pp = sourceInfo[i][4];
shiftedEnergy = (0.5)*(initialEnergy - 2*energyScaleFactor); //Pair production energy
if(abs(xyTry[0]-shiftedEnergy) < brThresh) {
return yMax*br*(pp/(pe+co+pp));
}
}
return 0.0;
}
Loading

0 comments on commit 11b1710

Please sign in to comment.