Commit e4f0a039 authored by jannick.wolters's avatar jannick.wolters
Browse files

added correct scaling to scattering cross sections

parent 39ca461d
......@@ -3,7 +3,7 @@ project( RTSN VERSION 0.1.0 LANGUAGES CXX )
set( CMAKE_CXX_STANDARD 17 )
set( CMAKE_CXX_STANDARD_REQUIRED ON )
set( CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -march=native -Wno-dev" )
set( CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3 -march=native -Wno-dev" )
set( CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -march=native -pg -no-pie" )
set( CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -Wall" )
......
......@@ -7,7 +7,7 @@ class Interpolation
{
public:
enum BOUNDARY { first_deriv, second_deriv };
enum TYPE { linear, cubic };
enum TYPE { linear, loglinear, cubic };
private:
Vector _x, _y;
......
......@@ -14,9 +14,8 @@ class Physics
private:
enum Element { H = 0, O = 1 };
std::vector<VectorVector> _xsScatteringH2O;
std::vector<VectorVector> _xsTotalH2O;
std::vector<VectorVector> _xsTransportH2O;
std::vector<VectorVector> _xsScatteringH2O; // only holds angular distribution, no cross section terms
std::vector<VectorVector> _xsTotalH2O; // equal to scatteringXS integrated over angle
VectorVector _stpowH2O;
const Vector _H20MassFractions{ 0.11189400, 0.88810600 };
......@@ -59,15 +58,6 @@ class Physics
*/
Vector GetStoppingPower( Vector energies );
/**
* @brief GetTransportXS gives back vector of vectors of stopping powers for materials defined by density and energies in vector energy
* @param energies is vector with energies
* @param density is vector with patient densities (at different spatial cells)
*/
VectorVector GetTransportXS( Vector energies, Vector density );
Vector GetTransportXSE( Vector energies );
/**
* @brief Physics constructor
*/
......
......@@ -95,6 +95,16 @@ void Interpolation::set_points( const Vector& x, const Vector& y, TYPE type ) {
_c[i] = ( _y[i + 1] - _y[i] ) / ( _x[i + 1] - _x[i] );
}
}
else if( type == loglinear ) {
_a.resize( n );
_b.resize( n );
_c.resize( n );
for( int i = 0; i < n - 1; i++ ) {
_a[i] = 0.0;
_b[i] = 0.0;
_c[i] = ( std::log( _y[i + 1] ) - std::log( _y[i] ) ) / ( _x[i + 1] - _x[i] );
}
}
else {
ErrorMessages::Error( "Invalid interpolation type!", CURRENT_FUNCTION );
}
......@@ -125,5 +135,8 @@ double Interpolation::operator()( double x ) const {
else {
interpol = ( ( _a[idx] * h + _b[idx] ) * h + _c[idx] ) * h + _y[idx];
}
return interpol;
if( _type == loglinear )
return std::exp( interpol );
else
return interpol;
}
......@@ -3,16 +3,12 @@
Physics::Physics( std::string fileName_H, std::string fileName_O, std::string fileName_stppower ) {
_xsScatteringH2O.resize( 2 );
_xsTotalH2O.resize( 2 );
_xsTransportH2O.resize( 2 );
LoadDatabase( fileName_H, fileName_O, fileName_stppower );
}
Physics::~Physics() {}
void Physics::LoadDatabase( std::string fileName_H, std::string fileName_O, std::string fileName_stppower ) {
VectorVector transport_XS_H;
VectorVector transport_XS_O;
VectorVector scattering_XS_H;
VectorVector scattering_XS_O;
......@@ -30,11 +26,9 @@ void Physics::LoadDatabase( std::string fileName_H, std::string fileName_O, std:
for( unsigned i = 0; i < headers_H.size(); ++i ) {
auto header = headers_H[i];
auto data = data_H[i];
if( header[1][0] == 7 && header[1][1] == 0 ) {
transport_XS_H = data;
}
// Integrated elastic scattering XS
else if( header[1][0] == 10 && header[1][1] == 0 ) {
if( header[1][0] == 10 && header[1][1] == 0 ) {
total_XS_H = data;
}
// Angular distribution large angle elastic scattering XS
......@@ -50,11 +44,8 @@ void Physics::LoadDatabase( std::string fileName_H, std::string fileName_O, std:
for( unsigned i = 0; i < headers_O.size(); ++i ) {
auto header = headers_O[i];
auto data = data_O[i];
if( header[1][0] == 7 && header[1][1] == 0 ) {
transport_XS_O = data;
}
// Integrated elastic scattering XS
else if( header[1][0] == 10 && header[1][1] == 0 ) {
if( header[1][0] == 10 && header[1][1] == 0 ) {
total_XS_O = data;
}
// Angular distribution large angle elastic scattering XS
......@@ -66,9 +57,6 @@ void Physics::LoadDatabase( std::string fileName_H, std::string fileName_O, std:
_xsScatteringH2O[H] = scattering_XS_H;
_xsScatteringH2O[O] = scattering_XS_O;
_xsTransportH2O[H] = transport_XS_H;
_xsTransportH2O[O] = transport_XS_O;
_xsTotalH2O[H] = total_XS_H;
_xsTotalH2O[O] = total_XS_O;
......@@ -84,13 +72,15 @@ std::vector<Matrix> Physics::GetScatteringXS( const Vector& energies, const Matr
angleVec[i * angle.columns() + j] = angle( i, j );
}
}
VectorVector outVec = GetScatteringXS( energies, angleVec );
Vector totalXS = GetTotalXSE( energies );
// rearrange output to matrix format
for( unsigned n = 0; n < energies.size(); ++n ) {
for( unsigned i = 0; i < angle.rows(); ++i ) {
for( unsigned j = 0; j < angle.columns(); ++j ) {
out[n]( i, j ) = outVec[n][i * angle.columns() + j];
out[n]( i, j ) = outVec[n][i * angle.columns() + j] * totalXS[n];
}
}
}
......@@ -224,36 +214,6 @@ Vector Physics::GetStoppingPower( Vector energies ) {
return stopping_power;
}
VectorVector Physics::GetTransportXS( Vector energies, Vector density ) {
VectorVector transport_XS( energies.size() );
Interpolation xsH( _xsTransportH2O[H][0], _xsTransportH2O[H][1], Interpolation::linear );
Interpolation xsO( _xsTransportH2O[O][0], _xsTransportH2O[O][1], Interpolation::linear );
for( unsigned i = 0; i < energies.size(); i++ ) {
for( unsigned i = 0; i < energies.size(); i++ ) {
transport_XS[i] = ( _H20MassFractions[H] * xsH( energies[i] ) * 1e-24 + _H20MassFractions[O] * xsO( energies[i] ) * 1e-24 ) * density;
}
}
return transport_XS;
}
Vector Physics::GetTransportXSE( Vector energies ) {
Vector transport_XS( energies.size() );
Interpolation xsH( _xsTransportH2O[H][0], _xsTransportH2O[H][1], Interpolation::linear );
Interpolation xsO( _xsTransportH2O[O][0], _xsTransportH2O[O][1], Interpolation::linear );
for( unsigned i = 0; i < energies.size(); i++ ) {
for( unsigned i = 0; i < energies.size(); i++ ) {
transport_XS[i] = ( _H20MassFractions[H] * xsH( energies[i] ) * 1e-24 + _H20MassFractions[O] * xsO( energies[i] ) * 1e-24 );
}
}
return transport_XS;
}
std::tuple<std::vector<VectorVector>, std::vector<VectorVector>> Physics::ReadENDL( std::string filename ) {
std::vector<VectorVector> _headers;
std::vector<VectorVector> _data;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment