randMC.cpp

#define _USE_MATH_DEFINES
#include <cmath>
#include "RandMC.h"
#include "random.h"
#include "material.h"
#include "AtomicData.h"
#include "Spectrum.h"
#include <iostream>


cRandMC::cRandMC()
{
	pRand = NULL;
	minEnergy = 1; // in keV ;
	maxEnergy = 1000;  // in keV
	noOfEnergySteps = 100;
}


cRandMC::~cRandMC()
{
}

void cRandMC::setRandBase(cRandKiss &rand) {
	pRand = &rand;
}

void cRandMC::prepare(cMaterial &mat, double tubeVoltage) {
	cAtomicData data;
	data.prepare();

	if (pRand == NULL) {
		throw runtime_error("pRand is null");
	}
	else {
		preparePhotonEnergy(tubeVoltage);
		prepareFreePath(mat);
		prepareInteractingElement(mat);

		for (unsigned i = 0; i < mat.getNoOfElements(); i++) {

			prepareInteraction(mat.getAtomicNumber(i));
			prepareSubShellCS(mat.getAtomicNumber(i));
			prepareAugerChar(mat.getAtomicNumber(i));
			prepareComptonPolar(mat.getAtomicNumber(i));
			prepareRayleighPolar(mat.getAtomicNumber(i));
		}
	}
}
double cRandMC::getPhotonEnergy() {
	return randPhotonEnergy.rand();
}
double cRandMC::getIsotropicPolar() {

	return acos(1 - 2 * pRand->randU());
}
double cRandMC::getIsotropicPolar(double thetaMax) {

	return acos(1 -(1 - cos(thetaMax)) * pRand->randU());
}
void cRandMC::getIsotropicAzimuth(double u[3]) {

	double azimuth; // new azimuthal angle
	double a ; // length on xy−plane
	azimuth = 2 * M_PI * pRand->randU ();
	a = sqrt (1 - u[2] * u[2] ) ;
	u[0] = a * sin (azimuth );
	u[1] = a * cos (azimuth );

}
double cRandMC::getFreePath(double energy) {

	unsigned i = 0; // index for mfp
	unsigned step; // actual step size
					// find index i
	for (step = step0Mfp; step; step >>= 1) {
		if ((i | step) < noOfEnergySteps && energy > meanFreePath[0][i|step]) {
			i |= step;
		}

		// apply linear interpolation
		return (meanFreePath[1][i] + (meanFreePath[1][i + 1] - meanFreePath[1][i]) *(energy - meanFreePath[0][i]) / (meanFreePath[0][i + 1] - meanFreePath[0][i])) * randExp.rand();
	}

}
unsigned cRandMC::getInteractingElement(double energy) {

	return (unsigned)randChooseElement.rand(energy);
}
unsigned cRandMC::getInteraction(unsigned Z, double energy) {

	return (unsigned)randInteraction[Z-1].rand(energy);
}
unsigned cRandMC::getPhotoSubshell(unsigned Z, double energy) {

	if (Z <= 2) {
		return 1;
	}
	else {
		return (unsigned)randSubshell[Z-1].rand(energy);
	}

}
double cRandMC::getAugerChar(unsigned Z, unsigned vacancy, unsigned &newVac1, unsigned &newVac2) {

	cAtomicData data;
	unsigned randVal;

	if (randAugerChar[Z - 1]->isPrepared()) {
		randVal = (unsigned)randAugerChar[Z-1][vacancy -1].rand();
		newVac1 = randVal & 0x3F;
		newVac2 = randVal >> 6;

		return data.getBindingEnergy(Z, vacancy) - data.getBindingEnergy(Z, newVac1) - data.getBindingEnergy(Z, newVac2);
	}
	else {
		newVac1 = 0;
		newVac2 = 0;
		return data.getBindingEnergy(Z, vacancy);
	}
}
double cRandMC::getComptonPolar(unsigned Z, double energy) {

	return randComPolAngle[Z-1].rand(energy);
}
double cRandMC::getRayleighPolar(unsigned Z, double energy) {

	return randRayPolAngle[Z-1].rand(energy);
}
void cRandMC::preparePhotonEnergy(double tubeVoltage) {

	double minEnergy;
	string spectrumName;
	cSpectrum spectrum;
	spectrum.readSpectrum("SRO33100ROT350.dat", tubeVoltage, minEnergy, spectrumName);


	randPhotonEnergy.reset();
	randPhotonEnergy.setRandomBase(*pRand);

	double step = (tubeVoltage - minEnergy) / spectrum.size();

	for (unsigned i = 0; i < spectrum.size(); i++) {
		randPhotonEnergy.addDistPoint(i * step, spectrum[i]);
	}
	randPhotonEnergy.addDistPoint(tubeVoltage, 0);
	randPhotonEnergy.prepare();

}
void cRandMC::prepareFreePath(cMaterial &mat) {

	meanFreePath[0].resize(noOfEnergySteps + 1);
	meanFreePath[1].resize(noOfEnergySteps + 1);

	for (unsigned k = 0; k <= noOfEnergySteps; k++) {

		meanFreePath[0][k] = minEnergy * exp(k * log(maxEnergy / minEnergy) / noOfEnergySteps);

		meanFreePath[1][k] = 1/ mat.getAttCoeff(minEnergy);
	}

	randExp.reset();
	randExp.setRandomBase(*pRand);

	for (step0Mfp = 1 << 31; !(step0Mfp & (noOfEnergySteps -1)); step0Mfp >>= 1);

}
void cRandMC::prepareInteractingElement(cMaterial &mat) {

	vector<double> dist;
	dist.resize(mat.getNoOfElements());

	for (unsigned i = 0; i < mat.getNoOfElements(); i++) {

		dist[i] = mat.getAtomicNumber(i) + 0.5;
	}

	randChooseElement.reset();
	randChooseElement.setRandomBase(*pRand);
	randChooseElement.setRandomValues(dist);

	double energy;

	for (unsigned k = 0; k <= noOfEnergySteps; k++) {

		energy = minEnergy * exp(k * log(maxEnergy / minEnergy) / noOfEnergySteps);

		for (unsigned m = 0; m < mat.getNoOfElements(); m++) {

			dist[m] = (mat.getAttCoeff(energy) * mat.getFraction(m));

			randChooseElement.addProbabilities(dist, energy);
		}
	}

	randChooseElement.prepare();
}
void cRandMC::prepareInteraction(unsigned Z) {

	cAtomicData data;
	vector<double> dist(3);
	randInteraction[Z-1].reset();
	randInteraction[Z-1].setRandomBase(*pRand);

	dist[0] = 1.5;
	dist[1] = 2.5;
	dist[2] = 3.5;

	randInteraction->setRandomValues(dist);

	double energy;

	for (unsigned k = 0; k <= noOfEnergySteps; k++) {

		energy = minEnergy * exp(k * log(maxEnergy / minEnergy) / noOfEnergySteps);

		dist[0] = data.getPhotoTotal(Z, energy);
		dist[1] = data.getComptonTotal(Z, energy);
		dist[2] = data.getRayleighTotal(Z, energy);

		randInteraction[Z - 1].addProbabilities(dist, energy);
	}
	randInteraction->prepare();
}
void cRandMC::prepareSubShellCS(unsigned Z) {

	vector<unsigned> subshell;
	vector<double> dist;

	if (Z > 2) {

		randSubshell[Z-1].reset();
		randSubshell[Z-1].setRandomBase(*pRand);

		cAtomicData data;
		data.prepare();
		data.getPhotoSubshells(Z, subshell);
		dist.resize(subshell.size());

		for (unsigned i = 0; i < subshell.size(); i++) {
			dist[i] = subshell[i] + 0.5;
		}
		randSubshell[Z - 1].setRandomValues(dist);
		double energy;

		for (unsigned k = 0; k <= noOfEnergySteps; k++) {

			energy = minEnergy * exp(k * log(maxEnergy / minEnergy) / noOfEnergySteps);

			for (unsigned m = 0; m < subshell.size(); m++) {

				dist[m] = data.getTotalCrossSection(subshell[m], energy);
			}
			randSubshell[Z - 1].addProbabilities(dist, energy);

		}
		randSubshell[Z - 1].prepare();
	}

}
void cRandMC::prepareAugerChar(unsigned Z) {

	vector<unsigned> subshell;
	vector<cAtomicData::tAugerEm> augerEm;
	vector<cAtomicData::tCharEm> charEm;

	cAtomicData data;
	data.getPhotoSubshells(Z, subshell);

	for (unsigned k = 0; k < subshell.size(); k++) {

		if (data.getBindingEnergy(Z, subshell[k]) > minEnergy) {

			data.getAugerEmissionTable(Z, subshell[k], augerEm);
			data.getCharEmissionTable(Z, subshell[k], charEm);

			randAugerChar[Z -1][subshell[k] -1].reset();
			randAugerChar[Z-1][subshell[k]-1].setRandomBase(*pRand);

			for (unsigned j = 0; j < augerEm.size(); j++) {

				randAugerChar[Z - 1][subshell[k] - 1].addDistPoint(augerEm[j].vacancy1 + (augerEm[j].vacancy2 << 6) + 0.5, augerEm[j].prob);
			}

			for (unsigned m = 0; m < charEm.size(); m++) {
				randAugerChar[Z - 1][subshell[k] - 1].addDistPoint(charEm[m].vacancy + 0.5, charEm[m].prob);
			}
			randAugerChar[Z - 1][subshell[k] - 1].prepare();
		}

	}
}
void cRandMC::prepareComptonPolar(unsigned Z) {

	double energy;
	double theta;
	double prob;
	cAtomicData data;
	randComPolAngle[Z-1].reset();
	randComPolAngle[Z-1].setRandomBase(*pRand);

	for (unsigned k = 0; k <= noOfEnergySteps; k++) {

		energy = minEnergy * exp(k * log(maxEnergy / minEnergy) / noOfEnergySteps);

		randComPolAngle[Z - 1].add(0, 0, energy);

		for (unsigned m = 1; m < 179; m++) {

			theta = m * M_PI / 180;

			prob = sin(theta) * data.getComptonDiff(Z, energy, theta);

			randComPolAngle[Z-1].add(theta, prob, energy);
		}
	}

	randComPolAngle[Z-1].add(M_PI, 0, energy);
	randComPolAngle[Z - 1].prepare();
}
void cRandMC::prepareRayleighPolar(unsigned Z) {

	double energy;
	double theta;
	double prob;
	cAtomicData data;
	randRayPolAngle[Z-1].reset();
	randRayPolAngle[Z-1].setRandomBase(*pRand);

	for (unsigned k = 0; k <= noOfEnergySteps; k++) {

		energy = minEnergy * exp(k * log(maxEnergy / minEnergy) / noOfEnergySteps);
		randRayPolAngle[Z - 1].add(0, 0, energy);

		for (unsigned m = 1; m < 179; m++) {

			theta = m * M_PI / 180;
			prob = sin(theta) * data.getComptonDiff(Z, energy, theta);

			randRayPolAngle[Z - 1].add(theta, prob, energy);
		}
		randRayPolAngle[Z - 1].add(M_PI, 0, energy);
	}
	randRayPolAngle[Z - 1].prepare();
}