<repo_name>gbzan/sim-xrf<file_sep>/python/snr/__init__.py
__author__ = "Yue"
__all__ = ["snr"]<file_sep>/Include/global.hpp
/* global.hpp */
#ifndef GLOABAL_HPP
#define GLOABAL_HPP

extern vector<double> ev_xrf, y_xrf;

#endif<file_sep>/python/__init__.py
__author__ = 'Yue'
__all__ = ["pyapi"]
path = __path__[0]
libpath = path + '/../Lib'

# Load .pyconfig file
configf = open(path + '/../.pyconfig')
config = {}
for line in configf.readlines():
    key, value = line.split('=')
    try:
        config[key.strip()] = [int(x) for x in value.split('#')[0].split(',')]
    except ValueError:
        config[key.strip()] = [float(x) for x in value.split('#')[0].split(',')]

configf.close()

plt_kwargs = {key: value for key, value in config.items()}
del (plt_kwargs['nout'])<file_sep>/src/dose.cpp
//dose.cpp
#include <cmath>

#include "dose.hpp"
#include "constants.hpp"

Dose::Dose(const Sample & sp, const Illumination & il) : n_photons(_n_photons), beam_cross_section(_beam_cross_section), dose_vec(_dose_vec), total_dose(_total_dose)
{
    _n_photons = il.n_photons;
	_beam_cross_section = il.beam_cross_section;
    _total_dose = 0;

	double _n_photons_calc, _beam_cross_section_calc;
    if ( n_photons > 0 && beam_cross_section > 0)
    {
		_n_photons_calc = n_photons;
		_beam_cross_section_calc = beam_cross_section;
    }
    else
    {
    	_n_photons_calc = 1;
		_beam_cross_section_calc = 1;
    }

    double joules_absorbed, matrix_mass, n_photons_in=_n_photons_calc, n_photons_out;
    for (auto ml : sp.layer_vec)
    {
        if (ml.thickness > 0)
        {
        	n_photons_out = n_photons_in*std::exp(-ml.mac_tot(il.ev0)*ml.thickness);
			joules_absorbed = (n_photons_in - n_photons_out)*il.ev0*eV_in_Joules;
			matrix_mass = ml.thickness*_beam_cross_section_calc*ml.density*1000/std::sin(il.psi);
			_dose_vec.push_back(joules_absorbed/matrix_mass);
			_total_dose += dose_vec.back();
			n_photons_in = n_photons_out;
        }
        else
        	_dose_vec.push_back(-1);
    }
}


Dose::~Dose()
{

}

Dose & Dose::operator=(const Dose & d)
{
	if (this != &d)
	{
		_n_photons = d.n_photons;
		_beam_cross_section = d.beam_cross_section;
		_dose_vec = d.dose_vec;
		_total_dose = d.total_dose;
	}
	return *this;
}

void Dose::show() const
{
	out(std::cout);
}

void Dose::out(std::ostream & ost) const
{
	ost << "# ======================================== #" << std::endl;
	ost << "# # # # # Radiation dosage: # # # # #" << std::endl;
	if ( n_photons > 0 && beam_cross_section > 0)
    {
		ost << "# For " << n_photons << " photon(s) incident, "
			<< "beam cross-sectional area = " << beam_cross_section << " cm^2 "
			<< "(d = " << std::sqrt(beam_cross_section/Pi)*2e7 << " nm):" << std::endl;
	}
	else
		ost << "# Dose per incident photon number density (1 photon/cm^2)" << std::endl;
	for (int i=0; i<dose_vec.size(); i++)
	    ost << "# \tLayer #" << i+1 << ":\t" << dose_vec[i] << " Gy" << std::endl;
	ost << "# \tTotal:\t" << total_dose << " Gy" << std::endl;
	ost << "# # # # # End of radiation dosage # # # # #" << std::endl;
	ost << "# ======================================== #" << std::endl;
}
<file_sep>/Include/geometry.hpp
//geometry.hpp
#ifndef GEOMETRY_HPP
#define GEOMETRY_HPP

#include <vector>
#include <ostream>
#include "constants.hpp"
#include "sample.hpp"

void spherical2cartesian(const double & r, const double & theta, const double & beta, double & x, double & y, double & z); //Convert spherical coordinates (r,theta_rad,beta_rad) to Cartesian coordinates (x,y,z).
double arctan2(const double & x);
void rotationx(const double & rad, double & x, double & y, double & z);
void rotationy(const double & rad, double & x, double & y, double & z);

class Illumination
{
private:
	double _ev0;
	double _psi;
	double _alpha;
	double _n_photons;
	double _beam_cross_section;
public:
	Illumination(double ev0_ = 1e4, double psi_ = Pi/4, double alpha_ = 0, double n_photons_=-1, double beam_cross_section_=-1);
	~Illumination();
	Illumination & operator=(const Illumination & il);
	const double & ev0;
	const double & psi;
	const double & alpha;
	const double & n_photons;
	const double & beam_cross_section;
	double psi_prime(const double & theta, const double & beta) const;
	void show() const;
	void out(std::ostream & ost) const;
};

class solid_angle
{
private:
	double _angle_range[4];
	std::vector<double> _theta;
	std::vector<double> _beta;
	double _theta_inc;
	double _beta_inc;
	double _subtend;
	void update();
public:
	solid_angle();
	solid_angle(const double *ar, double th_inc = Pi/180, double be_inc = Pi/180);
	~solid_angle();
	solid_angle & operator=(const solid_angle & sa);
	// const std::vector<double> & get_theta(double th_inc = -1);
	// const std::vector<double> & get_beta(double be_inc = -1);
	const double* const angle_range = _angle_range;
	const std::vector<double> & theta;
	const std::vector<double> & beta;
	const double & theta_inc;
	const double & beta_inc;
	const double & subtend;
	double domega(const double & theta) const;
	void show() const;
	void out(std::ostream & ost) const;
};

double atten_mono(const double & ev0, 
		const double & ev, 
		const double & psi, 
		const double & psiprime, 
		const Monolayer & c);
		
double atten_refl(const double & ev0, 
		const double & ev, 
		const double & psi, 
		const double & psiprime, 
		const Monolayer & ml);	
		
double atten_trans_in(const double & ev0, 
		const double & psi, 
		const Monolayer & ml);	
		
double atten_trans_out(const double & ev, 
		const double & psiprime, 
		const Monolayer & ml);

#endif<file_sep>/src/cs.cpp
#include <cmath>

#include "xraylib.h"

#include "constants.hpp"
#include "math.hpp"

inline double ev2nm(double ev) //Convert photon energy (in ev) to wavelength (in nm).
{
	return 1239.84187/ev;
}

double theta2x(double theta_rad, double ev) //theta in rad -> x in cm^-1
{
	return sin(theta_rad/2.)/(1.e-7*ev2nm(ev));
}


/*differential cross section for single electron*/

double thomson_pol(double theta_rad, double beta_rad)
{
	return sq(r_e)*(sq(cos(beta_rad)*cos(theta_rad))+sq(sin(beta_rad)));
}

double thomson_unpol(double theta_rad)
{
	return sq(r_e)*(1+sq(cos(theta_rad)))/2.;
}

double dcs_rayleigh_pol(double theta_rad, double beta_rad, double ev, int Z)
{
	return sq(FF_Rayl(Z, 1e-8*theta2x(theta_rad, ev)))*thomson_pol(theta_rad, beta_rad);
}
	
	
double dcs_rayleigh_unpol(double theta_rad, double ev, int Z)
{
	return sq(FF_Rayl(Z, 1e-8*theta2x(theta_rad, ev)))*thomson_unpol(theta_rad);
}	

	
/* Compton scattering cross sections */
	
/* compton dcs of single electron */
//dcs

double ev_scattered(double theta_rad, double ev) //return the energy in ev of the Compton scattered photons.
{
	return ev/(1+ev/m_e*(1-cos(theta_rad)));
}

double klein_nishina_unpol(double theta_rad, double ev) //Differential Klein-Nishina cross section in cm^2 for unpolarized radiation.

{
	double ev1 = ev_scattered(theta_rad, ev);
	return sq(r_e*(ev1/ev))/2.*(ev1/ev+ev/ev1-sq(sin(theta_rad)));
}

double klein_nishina_pol(double theta_rad, double beta_rad, double ev) //Differential Klein-Nishina cross section in cm^2 for polarized radiation.

{
	double ev1 = ev_scattered(theta_rad, ev);
	return sq(r_e*(ev1/ev))/2.*(ev1/ev+ev/ev1-2*sq(sin(theta_rad)*cos(beta_rad)));
}

//total cs
double klein_nishina_total_col(double ev) // Total Klein-Nishina collision cross section in cm^2.
{
	double a = ev/m_e;
	return 2*Pi*sq(r_e)*((1+a)/(a*a*a)*(2*a*(1+a)/(1+2*a)-log(1+2*a))+log(1+2*a)/(2*a)-(1+3*a)/sq(1+2*a));
}


double klein_nishina_total_sca(double ev) //Total Klein-Nishina scattering cross section in cm^2.
{
	double a = ev/m_e;
	return Pi*sq(r_e)*(log(1+2*a)/(a*a*a)+2*(1+a)*(2*sq(a)-2*a-1)/sq(a*(1+2*a))+8*sq(a)/(3*sq(1+2*a)*(1+2*a)));
}
//compton dcs of atom
double dcs_compton_pol(double theta_rad, double beta_rad, double ev, int Z) //Differential Compton cross section in cm^2 for polarized radiation.
{
	return klein_nishina_pol(theta_rad, beta_rad, ev)*SF_Compt(Z, theta2x(theta_rad,ev)*1e-8);
}
	

double dcs_compton_unpol(double theta_rad, double ev, int Z) //Differential Compton cross section in cm^2 for unpolarized radiation.
{
	return klein_nishina_unpol(theta_rad, ev)*SF_Compt(Z, theta2x(theta_rad,ev)*1e-8);
}<file_sep>/Include/xrf.hpp
/* xrf.hpp */
#ifndef XRF_HPP
#define XRF_HPP
#include <vector>

using namespace std;

#define K_LINES -29
#define	L_LINES -113
#define M_LINES -219
#define N_LINES -324
#define O_LINES -374
#define P_LINES -383

vector<int> mac_xrf(double ev0, int Z, vector<double> & ev_vec, vector<double> & y_vec, double weight = 1.0, int line_end = P_LINES);

#endif


<file_sep>/Include/dose.hpp
//dose.hpp
#ifndef DOSE_HPP
#define DOSE_HPP

#include <vector>
#include <iostream>

#include "geometry.hpp"
#include "sample.hpp"

class Dose
{
private:
	double _n_photons;
	double _beam_cross_section;
	double _total_dose;
	std::vector<double> _dose_vec;
public:
	Dose(const Sample & sp, const Illumination & il);
	~Dose();
	Dose & operator=(const Dose & d);
	const double & n_photons;
	const double & beam_cross_section;
	const double & total_dose;
	const std::vector<double> & dose_vec;
	void show() const;
	void out(std::ostream & ost) const;
};


#endif<file_sep>/simpy.py
# simpy.py
import python.pyapi as pyapi
from python import plt_kwargs


if __name__ == "__main__":
    import sys
    import time

    # Input and output file names
    try:
        input_file = sys.argv[1]
        try:
            output_file = sys.argv[2]
        except IndexError:
            output_file = "./output.txt"
    except IndexError:
        input_file = "./input.txt"
        output_file = "./output.txt"

    start = time.time()
    spec = pyapi.calc(input_file=input_file, output_file=output_file)
    print 'Calculation finished.'
    print "Time used:", time.time() - start, 'sec.'
    spec.show(**plt_kwargs)


    # import matplotlib.pyplot as plt
    # plt.plot(spec.y_sep.T)
    # plt.plot(spec.y_vec)
    # plt.yscale('log')
    # plt.ylim(1e-14,1e-4)
    # plt.legend(spec.labels+['total'])
    # plt.show()<file_sep>/src/input.cpp
/*input.cpp*/
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>

#include "input.hpp"
#include "constants.hpp"
#include "math.hpp"

inline std::string trim_comment(const std::string & s, const std::string & delimiter="#")
{
    if (s.empty())
        return s;
    else
        return s.substr(0, s.find(delimiter));
}

inline std::string trim_right(
  const std::string & s,
  const std::string & delimiters = " \f\n\r\t\v" )
{
	if (s.empty())
		return s;
	return s.substr(0, s.find_last_not_of(delimiters) + 1);
}

inline std::string trim_left(
  const std::string & s,
  const std::string & delimiters = " \f\n\r\t\v" )
{
	if (s.empty())
		return s;
	return s.substr(s.find_first_not_of(delimiters));
}

inline std::string trim(
  const std::string & s,
  const std::string & delimiters = " \f\n\r\t\v" )
{
	return trim_left(trim_right(trim_comment(s), delimiters), delimiters);
}

inline std::vector<int> parse_vec_int(const std::string & s)
{
	std::vector<int> v;
	std::istringstream is(s);
	std::string t;
	while(std::getline(is, t, ','))
	// {
		// std::cout << trim(t) << std::endl;
		v.push_back(std::stoi(trim(t)));
	// }
	return v;
}

inline std::vector<double> parse_vec_double(const std::string & s)
{
	std::vector<double> v;
	std::istringstream is(s);
	std::string t;
	while(std::getline(is, t, ','))
		v.push_back(std::stod(trim(t)));
	return v;
}

inline bool parse_parameter(const std::string & line, std::string & name, std::string & value, const std::string & sep = "=")
{
	std::size_t pos = line.find(sep);
	if (pos == std::string::npos)
		return false;
	else
	{
		name = trim(line.substr(0, pos));
		value = trim(line.substr(pos+1));
		return true;
	}
}

inline double deg_to_rad(const double & deg)
{
	return deg*Pi/180;
}


void readfile(std::string filename, Sample & sp, Illumination & il, solid_angle & omega, Detector & det)
{
	std::ifstream fin;
	std::string s;
	fin.open(filename.c_str());

	std::size_t pos;
	while (getline(fin, s))
	{
		if (trim(s).empty())
			continue;
		else if (s.find("Sample:") != std::string::npos)
		{
			std::cout << "Reading in sample data..." << std::endl;
			// read sample
			while (getline(fin, s) && !trim(s).empty())
			{
				int layer = std::stoi(trim(s.substr(s.find("layer")+5, s.find(":")-6)));
				std::cout << "Reading layer #" << layer << "..." << std::endl;
				//read layer
				std::vector<int> Z_vec;
				std::vector<double> p_vec;
				double thickness, density;
				while (getline(fin, s) && !trim(s).empty())
				{
					std::string name, value;
					if (parse_parameter(s, name, value))
					{
						if (name == "Z")
						{
							Z_vec = parse_vec_int(value);					
						}
						else if (name == "p")
						{
							p_vec = parse_vec_double(value);
						}
						else if (name == "thickness")
						{
							thickness = std::stod(value);
						}
						else if (name == "density")
						{
							density = std::stod(value);
						}
					}
				}
				sp.add_layer(Monolayer(Z_vec, p_vec, density, thickness, layer));
			}
			
		}
		else if (s.find("Detector:") != std::string::npos)
		{
			std::cout << "Reading in detector data..." << std::endl;
			// read detector
			double ev_offset, ev_gain, n_channels, 
				noise, fano, gamma, fs, ft, 
				thickness, density;
			std::string material;
			while (getline(fin, s) && !trim(s).empty())
			{
				std::string name, value;
				if(parse_parameter(s, name, value))
				{
					if (name == "ev_offset")
					{
						ev_offset = std::stod(value);
					}
					else if (name == "ev_gain")
					{
						ev_gain = std::stod(value);
					}
					else if (name == "n_channels")
					{
						n_channels = std::stod(value);
					}
					else if (name == "noise")
					{
						noise = std::stod(value);
					}
					else if (name == "fano")
					{
						fano = std::stod(value);
					}
					else if (name == "gamma")
					{
						gamma = std::stod(value);
					}
					else if (name == "fs")
					{
						fs = std::stod(value);
					}
					else if (name == "ft")
					{
						ft = std::stod(value);
					}
					else if (name == "material")
					{
						material = value;
					}
					else if (name == "thickness")
					{
						thickness = std::stod(value);
					}
					else if (name == "density")
					{
						density = std::stod(value);
					}
				}
				
				det = Detector(Channel(ev_offset, ev_gain, n_channels), 
					Response(noise, fano, gamma, fs, ft, ev_gain), 
					Window(material, thickness, density));
			}
		}
		else if (s.find("Illumination:") != std::string::npos)
		{
			std::cout << "Reading in illumination data..." << std::endl;
			// read illumination
			double ev0, psi, alpha, n_photons=-1, beam_cross_section=-1;
			while (getline(fin, s) && !trim(s).empty())
			{
				std::string name, value;
				if(parse_parameter(s, name, value))
				{
					if (name == "ev0")
					{
						ev0 = std::stod(value);
					}
					else if (name == "psi")
					{
						psi = deg_to_rad(std::stod(value));
					}
					else if (name == "alpha")
					{
						alpha = deg_to_rad(std::stod(value));
					}
					else if (name == "n_photons")
					{
						n_photons = std::stod(value);
					}
					else if (name == "beam_cross_section")
					{
						beam_cross_section = std::stod(value);
					}
					else if (name == "beam_diameter")
					{
						beam_cross_section = sq(std::stod(value)/2.)*Pi;
					}
				}
				
			}
			il = Illumination(ev0, psi, alpha, n_photons, beam_cross_section);
		}
		
		else if (s.find("Solid angle:") != std::string::npos)
		{
			std::cout << "Reading in solid angle data..." << std::endl;
			// read solid angle
			std::vector<double> ar;
			double theta_inc, beta_inc;
			while (getline(fin, s) && !trim(s).empty())
			{
				std::string name, value;
				if(parse_parameter(s, name, value))
				{
					if (name == "angle_range")
					{
						ar = parse_vec_double(value);
						for (std::vector<double>::iterator i = ar.begin(); i < ar.end(); i++)
							*i = deg_to_rad(*i);
					}
					else if (name == "theta_inc")
					{
						theta_inc = deg_to_rad(std::stod(value));
					}
					else if (name == "beta_inc")
					{
						beta_inc = deg_to_rad(std::stod(value));
					}
				}
			}
			omega = solid_angle(ar.data(), theta_inc, beta_inc);
		}
	}
	std::cout << "Done!" << std::endl;
}


<file_sep>/python/snr/pysnip.py
import numpy as np
from scipy.optimize import curve_fit


def FWHM(x, noise=100, fano=0.114):
    sigma = np.sqrt((noise / 2.3548) ** 2 + 3.58 * fano * x)
    return 2.3548 * sigma


def fit_FWHM(x, F):
    def _FWHM(x, noise, fano):
        return (noise / 2.3548) ** 2 + 3.58 * fano * x

    popt, pcov = curve_fit(_FWHM, x, (F / 2.3548) ** 2, p0=[100, 0.114])
    return popt


def energy_to_channel(energy, offset=2.97, gain=12.26952):
    return 1. * (energy - offset) / gain


# # # Low statistics digital filter

def lsdf(E, y, FWHM=FWHM,
         f=1.5,
         A=75,
         M=10,
         r=1.3):
    def _reduce(x, length_start):
        for i in range(length_start):
            length = length_start - i
            if x < length:
                raise IndexError
            L = y[x - length:x].sum()
            R = y[x + 1:x + length + 1].sum()
            S = y[x] + L + R
            slope = (R + 1.) / (L + 1.)
            if S < M or S < A * np.sqrt(y[x]) or (1. / r <= slope <= r):
                return S / (2. * length + 1)
        print 'Not found for x = %d!' % x
        return y[x]

    y_out = y.copy()

    for x in range(len(E)):
        try:
            len_0 = int(energy_to_channel(f * FWHM(E[x]), E[0], E[1] - E[0]))
            y_out[x] = _reduce(x, len_0)
        except IndexError:
            pass
    return y_out


# # # Peak-clipping

def snip(E, y, FWHM=FWHM, offset=0., gain=10., **kwargs):
    det = kwargs.get('detector')
    loops = kwargs.get('loops', 24)
    end_loops = kwargs.get('end_loops', 8)
    reduce_factor = kwargs.get('reduce_factor', np.sqrt(2))
    factor = kwargs.get('factor', 2)

    if det is not None:
        FWHM = det.response.FWHM
        offset = det.channel.offset
        gain = det.channel.gain

    def G(y):
        return np.log(np.log(y + 1) + 1)

    def w(x, factor=2):
        return energy_to_channel(factor * FWHM(E[x]), offset=offset, gain=gain)

    def G_inv(z):
        return np.exp(np.exp(z) - 1) - 1

    z_out = G(y)
    for i in range(loops):
        if i >= loops - end_loops:
            factor /= 1. * reduce_factor
        z = z_out.copy()
        for x in range(len(E)):
            try:
                _w = w(x, factor=factor)
                if _w > x:
                    raise IndexError
                z_bar = (z[x + _w] + z[x - _w]) / 2.
                z_out[x] = min(z[x], z_bar)
            except IndexError:
                pass
    return G_inv(z_out)


<file_sep>/src/detector.cpp
/*detector.cpp*/
#include <iostream>
#include<cmath>

#include "xraylib.h"
#include "detector.hpp"
#include "math.hpp"
#include "constants.hpp"

#define MP 1

Channel::Channel(double ev_offset_, double ev_gain_, int n_channels_) 
	: ev_offset(_ev_offset), ev_gain(_ev_gain), n_channels(_n_channels)
{
	_ev_offset = ev_offset_;
	_ev_gain = ev_gain_;
	_n_channels = n_channels_;
}

Channel & Channel::operator=(const Channel & ch)
{
	if (this != & ch)
	{
		_ev_offset = ch.ev_offset;
		_ev_gain = ch.ev_gain;
		_n_channels = ch.n_channels;
	}
	return *this;
}

int Channel::ev_to_channel(double ev) const
{
	return int((ev - ev_offset)/ev_gain + 0.5);
}

void Channel::bin(const std::vector<double> & ev_raw, const std::vector<double> & y_raw, std::vector<double> & y_binned) const
{
	if (y_binned.size() != n_channels)
	{
		std::cout << "Error: Bin size does not match the number of channels!" << std::endl;
		return;
	}
	#pragma omp parallel for if(MP)
	for (int i=0; i<ev_raw.size(); i++)
		y_binned[ev_to_channel(ev_raw[i])] += y_raw[i];
}

void Channel::bin(const double & ev_raw, const double & y_raw, std::vector<double> & y_binned) const
{
	if (y_binned.size() != n_channels)
	{
		std::cout << "Error: Bin size does not match the number of channels!" << std::endl;
		return;
	}
	y_binned[ev_to_channel(ev_raw)] += y_raw;
}

void Channel::bin(const std::vector<double> & ev_raw, const std::vector<double> & y_raw, std::vector<double> & y_binned, std::vector<double> & y_separate, const std::vector<int> & row) const
{
	if (y_binned.size() != n_channels)
	{
		std::cout << "Error: Bin size does not match the number of channels!" << std::endl;
		return;
	}
	int n_old = y_separate.size();
	y_separate.resize(n_old+n_channels);
	#pragma omp parallel for collapse(2) if(MP)
	for (int i = 0; i < row.size()-1; i++)
		for (int j = row[i]; j < ((row[i+1] > 0) ? row[i+1] : ev_raw.size()); j++)
		{
			int chn = ev_to_channel(ev_raw[j]);
			y_binned[chn] += y_raw[j];
			// for separate output
			y_separate[n_old + i*n_channels + chn] = y_raw[j];
		}
}

void Channel::bin(const double & ev_raw, const double & y_raw, std::vector<double> & y_binned, std::vector<double> & y_separate) const
{
	if (y_binned.size() != n_channels)
	{
		std::cout << "Error: Bin size does not match the number of channels!" << std::endl;
		return;
	}
	int n_old = y_separate.size();
	y_separate.resize(n_old+n_channels);
	int chn = ev_to_channel(ev_raw);
	y_binned[chn] += y_raw;
	y_separate[n_old+chn] = y_raw;
}

void Channel::show() const
{
    out(std::cout);
}

void Channel::out(std::ostream & ost) const
{
	ost << "# Detector channel parameters:" << std::endl;
	ost << "# \tev_offset = " << ev_offset << std::endl;
	ost << "# \tev_gain = " << ev_gain << std::endl;
	ost << "# \tn_channels = " << n_channels << std::endl;
    ost << std::endl;
}


Response::Response(double noise_,
		double fano_,
		double gamma_,
		double fs_,
		double ft_,
		double ev_gain_)
		: noise(_noise), fano(_fano), gamma(_gamma), fs(_fs), ft(_ft), ev_gain(_ev_gain)
{
	_noise = noise_;
	_fano = fano_;
	_gamma = gamma_;
	_fs = fs_;
	_ft = ft_;
	_ev_gain = ev_gain_;
}

Response & Response::operator=(const Response & r)
{
	if (this != & r)
	{
		_fano = r.fano;
		_gamma = r.gamma;
		_fs = r.fs;
		_ft = r.ft;
		_ev_gain = r.ev_gain;
	}
	return *this;
}

void Response::set_gain(const double ev_gain)
{
	_ev_gain = ev_gain;
}

double Response::FWHM(double ev)
{
	return 2.3548*std::sqrt(sq(noise/2.3548)+3.58*fano*ev);
}

void Response::show() const
{
    out(std::cout);
}

void Response::out(std::ostream & ost) const
{
	ost << "# Detector response parameters:" << std::endl;
	ost << "# \tfano = " << fano << std::endl;
	ost << "# \tgamma = " << gamma << std::endl;
	ost << "# \tfs = " << fs << std::endl;
	ost << "# \tft = " << ft << std::endl;
    ost << std::endl;
}



Window::Window(std::string material_ ,
	double thickness_,
	double density_)
	: material(_material), thickness(_thickness), density(_density)
{
	_material = material_;
	_thickness = thickness_;
	_density = density_;
}

Window & Window::operator=(const Window & w)
{
	if (this != & w)
	{
		_material = w.material;
		_thickness = w.thickness;
		_density = w.density;
	}	
	return *this;
}

double Window::transmission(double ev) const
{
	double _mac = CS_Total_CP(material.c_str(), ev/1000.);
	return std::exp(-_mac*density*thickness);
}

void Window::show() const
{
    out(std::cout);
}

void Window::out(std::ostream & ost) const
{
	ost << "# Detector window parameters:" << std::endl;
	ost << "# \tmaterial = " << material << std::endl;
	ost << "# \tthickness = " << thickness << std::endl;
	ost << "# \tdensity = " << density << std::endl;
    ost << std::endl;
}


Detector::Detector(Channel channel_,
			Response response_,
			Window window_)
			: channel(_channel), response(_response), window(_window)
{
	_channel = channel_;
	_response = response_;
	_window = window_;
	_response.set_gain(channel.ev_gain);
}

Detector & Detector::operator=(const Detector & d)
{
	if (this != & d)
	{
		_channel = d.channel;
		_response = d.response;
		_window = d.window;
	}
	return *this;
}

void Detector::genspec(const std::vector<double> & ev_raw, const std::vector<double> & y_raw, std::vector<double> & y_binned, bool det_response, bool det_window) const
{
	if (y_binned.size() != channel.n_channels)
	{
		std::cout << "Error: Bin size does not match the number of channels!" << std::endl;
		return;
	}
	//apply detector response
	if (det_response)
	{
		std::cout << "Generating spectrum with detector response..." << std::endl;
		double sigma, y_temp;
		for (std::vector<double>::const_iterator ev = ev_raw.begin(), y = y_raw.begin();
			ev < ev_raw.end(); ev++, y++)
		{
			sigma = std::sqrt(sq(response.noise/2.3548)+3.58*response.fano*(*ev));
			y_temp = *y;
			//apply detector filter window
			if (det_window)
				y_temp *= window.transmission(*ev);
			#pragma omp parallel for if(MP)
			for (int i = 0; i < y_binned.size(); i++)
			{
				double ev_ch = channel.ev_offset + i*channel.ev_gain;
				//Gaussian
				y_binned[i] += y_temp*channel.ev_gain/(sigma*std::sqrt(2*Pi))*std::exp(-sq((*ev)-ev_ch)/(2*sq(sigma)));
				//Step function
				y_binned[i] += response.fs*y_temp*channel.ev_gain/(2.*(*ev))*std::erfc((ev_ch-(*ev))/(std::sqrt(2)*sigma));
				//Tailing function
				y_binned[i] += response.ft*y_temp*channel.ev_gain/(2.*response.gamma*sigma*std::exp(-1./(2*sq(response.gamma))))*std::exp((ev_ch-(*ev))/(response.gamma*sigma))*std::erfc((ev_ch-(*ev))/(std::sqrt(2)*sigma)+1/(std::sqrt(2)*response.gamma));
			}
		}
	}
	else
	{
		std::cout << "Generating spectrum without detector response..." << std::endl;
		if (det_window)//apply detector filter window
			#pragma omp parallel for if(MP)
			for (int i = 0;	i < ev_raw.size(); i++)
				y_binned[channel.ev_to_channel(ev_raw[i])] += (y_raw[i])*window.transmission(ev_raw[i]);
		else
			channel.bin(ev_raw, y_raw, y_binned);
	}
}

void Detector::genspec(const double & ev_raw, const double & y_raw, std::vector<double> & y_binned, bool det_response, bool det_window) const
{
	if (y_binned.size() != channel.n_channels)
	{
		std::cout << "Error: Bin size does not match the number of channels!" << std::endl;
		return;
	}
	//apply detector response
	if (det_response)
	{
		std::cout << "Generating spectrum with detector response..." << std::endl;
		double sigma, y_temp;
		sigma = std::sqrt(sq(response.noise/2.3548)+3.58*response.fano*ev_raw);
		y_temp = y_raw;
		//apply detector filter window
		if (det_window)
			y_temp *= window.transmission(ev_raw);
		for (int i = 0; i < y_binned.size(); i++)
		{
			double ev_ch = channel.ev_offset + i*channel.ev_gain;
			//Gaussian
			y_binned[i] += y_temp*channel.ev_gain/(sigma*std::sqrt(2*Pi))*std::exp(-sq(ev_raw-ev_ch)/(2*sq(sigma)));
			//Step function
			y_binned[i] += response.fs*y_temp*channel.ev_gain/(2.*ev_raw)*std::erfc((ev_ch-ev_raw)/(std::sqrt(2)*sigma));
			//Tailing function
			y_binned[i] += response.ft*y_temp*channel.ev_gain/(2.*response.gamma*sigma*std::exp(-1./(2*sq(response.gamma))))*std::exp((ev_ch-ev_raw)/(response.gamma*sigma))*std::erfc((ev_ch-ev_raw)/(std::sqrt(2)*sigma)+1/(std::sqrt(2)*response.gamma));
		}
	}
	else
	{
		std::cout << "Generating spectrum without detector response..." << std::endl;
		if (det_window)//apply detector filter window
			y_binned[channel.ev_to_channel(ev_raw)] += y_raw*window.transmission(ev_raw);
		else
			channel.bin(ev_raw, y_raw, y_binned);
	}
}

void Detector::genspec(const std::vector<double> & ev_raw, const std::vector<double> & y_raw, std::vector<double> & y_binned, std::vector<double> & y_separate, const std::vector<int> & row, bool det_response, bool det_window) const
{
	if (y_binned.size() != channel.n_channels)
	{
		std::cout << "Error: Bin size does not match the number of channels!" << std::endl;
		return;
	}
	//apply detector response
	if (det_response)
	{
		std::cout << "Generating spectrum with detector response..." << std::endl;

		std::vector<double>::iterator ys0; // ys0 pointing to the first channel
		int n_old = y_separate.size();
		y_separate.resize(n_old+(row.size()-1)*channel.n_channels);
		ys0 = y_separate.begin() + n_old;

		int i = 0;
		double sigma, y_temp;
		for (std::vector<double>::const_iterator ev = ev_raw.begin(), y = y_raw.begin();
			ev < ev_raw.end(); ev++, y++)
		{
			// for separate output
			if (ev - ev_raw.begin() == row[i+1])
				i++;
			
			sigma = std::sqrt(sq(response.noise/2.3548)+3.58*response.fano*(*ev));
			y_temp = *y;
			//apply detector filter window
			if (det_window)
				y_temp *= window.transmission(*ev);
			#pragma omp parallel for if(MP)
			for (int j = 0; j < y_binned.size(); j++)
			{
				double y_temp1 = 0;
				double ev_ch = channel.ev_offset + j*channel.ev_gain;
				
				//Gaussian
				y_temp1 += y_temp*channel.ev_gain/(sigma*std::sqrt(2*Pi))*std::exp(-sq((*ev)-ev_ch)/(2*sq(sigma)));
				//Step function
				y_temp1 += response.fs*y_temp*channel.ev_gain/(2.*(*ev))*std::erfc((ev_ch-(*ev))/(std::sqrt(2)*sigma));
				//Tailing function
				y_temp1 += response.ft*y_temp*channel.ev_gain/(2.*response.gamma*sigma*std::exp(-1./(2*sq(response.gamma))))*std::exp((ev_ch-(*ev))/(response.gamma*sigma))*std::erfc((ev_ch-(*ev))/(std::sqrt(2)*sigma)+1/(std::sqrt(2)*response.gamma));
				
				// Total output
				y_binned[j] += y_temp1;
				//Separate output
				*(ys0 + i*channel.n_channels + j) += y_temp1;
			}
		}
	}
	else
	{
		std::cout << "Generating spectrum without detector response..." << std::endl;
		std::vector<double>::iterator ys0; // ys0 pointing to the first channel
		int n_old = y_separate.size();
		y_separate.resize(n_old+(row.size()-1)*channel.n_channels);
		ys0 = y_separate.begin() + n_old;
		int i = 0;
		if (det_window)//apply detector filter window
		{
			#pragma omp parallel for collapse(2) if(MP)
			for (int i = 0; i < row.size()-1; i++)
				for (int j = row[i]; j < ((row[i+1] > 0) ? row[i+1] : ev_raw.size()); j++)
				{
					int chn = channel.ev_to_channel(ev_raw[j]);
					double y_temp = (y_raw[j])*window.transmission(ev_raw[j]);
					y_binned[chn] += y_temp;
					// for separate output
					*(ys0 + i*channel.n_channels + chn) = y_temp;
				}
		}
		else
			channel.bin(ev_raw, y_raw, y_binned, y_separate, row);
	}
}

void Detector::genspec(const double & ev_raw, const double & y_raw, std::vector<double> & y_binned, std::vector<double> & y_separate, bool det_response, bool det_window) const
{
	if (y_binned.size() != channel.n_channels)
	{
		std::cout << "Error: Bin size does not match the number of channels!" << std::endl;
		return;
	}
	//apply detector response
	if (det_response)
	{
		std::cout << "Generating spectrum with detector response..." << std::endl;
		std::vector<double>::iterator ys0; // ys0 pointing to the first channel
		int n_old = y_separate.size();
		y_separate.resize(n_old+channel.n_channels);
		ys0 = y_separate.begin() + n_old;
		double sigma, y_temp;
		sigma = std::sqrt(sq(response.noise/2.3548)+3.58*response.fano*ev_raw);
		y_temp = y_raw;
		//apply detector filter window
		if (det_window)
			y_temp *= window.transmission(ev_raw);
		#pragma omp parallel for if(MP)
		for (int i = 0; i < y_binned.size(); i++)
		{
			double y_temp1 = 0;
			double ev_ch = channel.ev_offset + i*channel.ev_gain;
			//Gaussian
			y_temp1 += y_temp*channel.ev_gain/(sigma*std::sqrt(2*Pi))*std::exp(-sq(ev_raw-ev_ch)/(2*sq(sigma)));
			//Step function
			y_temp1 += response.fs*y_temp*channel.ev_gain/(2.*ev_raw)*std::erfc((ev_ch-ev_raw)/(std::sqrt(2)*sigma));
			//Tailing function
			y_temp1 += response.ft*y_temp*channel.ev_gain/(2.*response.gamma*sigma*std::exp(-1./(2*sq(response.gamma))))*std::exp((ev_ch-ev_raw)/(response.gamma*sigma))*std::erfc((ev_ch-ev_raw)/(std::sqrt(2)*sigma)+1/(std::sqrt(2)*response.gamma));
			
			// Total output
			y_binned[i] += y_temp1;
			//Separate output
			*(ys0+i) = y_temp1;
		}
	}
	else
	{
		std::cout << "Generating spectrum without detector response..." << std::endl;
		std::vector<double>::iterator ys0; // ys0 pointing to the first channel
		int n_old = y_separate.size();
		y_separate.resize(n_old+channel.n_channels);
		ys0 = y_separate.begin() + n_old;
		
		if (det_window)//apply detector filter window
		{
			int chn = channel.ev_to_channel(ev_raw);
			double y_temp = y_raw*window.transmission(ev_raw);
			y_binned[chn] += y_temp;
			*(ys0+chn) = y_temp;
		}
		else
			channel.bin(ev_raw, y_raw, y_binned, y_separate);
	}
}

void Detector::show() const
{
    out(std::cout);
}

void Detector::out(std::ostream & ost) const
{
    ost << "# ======================================== #" << std::endl;
	ost << "# # # # # Detector parameters # # # # #" << std::endl;
	channel.out(ost);
	response.out(ost);
	window.out(ost);
	ost << "# # # # # End of detector parameters # # # # #" << std::endl;
    ost << "# ======================================== #" << std::endl;
}<file_sep>/python/snr/snr.py
import numpy as np
import snip


def FWHM(x, noise=100, fano=0.114):
    sigma = np.sqrt((noise / 2.3548) ** 2 + 3.58 * fano * x)
    return 2.3548 * sigma


def pnb(ev, y, bg, peak_center, peak_width, bg_range=None, normalize=True):
    """

    :param ev: Energy array in eV.
    :param y: Total spectrum array. First dimension must be energy/channel number.
    :param bg: Background spectrum array. First dimension must be energy/channel number.
    :param peak_center: peak center in eV.
    :param peak_width: peak width in eV.
    :param bg_range: Range in which background is sampled. List with length of multiples of two.
    :param normalize: Boolean. If true, background is normalized with respect to the width of bg compared to width of y.
    :return: P,B in P/B.
    """

    def ev2ch(x):
        return np.abs(ev - x).argmin()

    p = np.array([ev2ch(peak_center - peak_width / 2.), ev2ch(peak_center + peak_width / 2.)])
    if bg_range is None:
        b = p
    else:
        b = np.array([ev2ch(x) for x in bg_range])

    Nt = np.sum(y[p[0]:p[1]])
    B = 0
    nchB = 0
    for i in range(0, len(b), 2):
        B += np.sum(bg[b[i]:b[i + 1]])
        nchB = b[i + 1] - b[i]
    if normalize:
        B *= 1. * (p[1] - p[0]) / nchB
    return Nt - B, B


def pnb_fixedwidth(ev, y, peak_center, peak_width, bg_range):
    return pnb(ev, y, y, peak_center, peak_width, bg_range)


def pnb_snip(ev, y, peak_center, peak_width=100, FWHM=snip.FWHM, offset=0., gain=10., **kwargs):
    bg = snip.snip(ev, y, FWHM, offset, gain, **kwargs)
    return pnb(ev, y, bg, peak_center, peak_width)


def pnb_ascalc(ev, y, total, peak_center, peak_width):
    return pnb(ev, total, total - y, peak_center, peak_width, normalize=False)


def snr(P, B):
    return P/np.sqrt(P + 2*B)<file_sep>/src/cwrapper.cc
#include <iostream>
#include <fstream>
#include <string>

#include "spectrum.hpp"
#include "input.hpp"
#include "dose.hpp"
extern "C"
{
	// __declspec(dllexport) 
	void sim(char * input_file, 
	char * output_file, 
	// Total
	double * y_vec,
	// Separate
	double * y_sep,
	// XRF
	int * Z_vec,
	int * row,
	int * lines,
	double * xrf_ev,
	double * xrf_y,
	// Compton
	double * comp_ev,
	double * comp_y,
	// Rayleigh
	double * ray_y,
	// Detector
	double * det_,
	int * n_channels,
	char * win_mat,
	// Illumination
	double * il_,
	// Solid angle
	double * sa,
	double * dose,
	int * nout);
}


void sim(char * input_file, 
	char * output_file, 
	// Total
	double * y_vec,
	// Separate
	double * y_sep,
	// XRF
	int * Z_vec,
	int * row,
	int * lines,
	double * xrf_ev,
	double * xrf_y,
	// Compton
	double * comp_ev,
	double * comp_y,
	// Rayleigh
	double * ray_y,
	// Detector
	double * det_,
	int * n_channels,
	char * win_mat,
	// Illumination
	double * il_,
	// Solid angle
	double * sa,
	// Radiation dose
	double * dose,
	int * nout)
{
	Sample sp;
	Illumination il;
	solid_angle omega;
	Detector det;
	
	// Read input file
	readfile(input_file, sp, il, omega, det);
	
	// Calculate spectrum
	Spectrum spec(sp, il, omega, det, true, true, true);
	Dose ds(sp, il);
	
	// Save results to file
	static std::ofstream fout;
	fout.open(output_file);
	//Outputs
	spec.out(fout);
	fout << std::endl;
	ds.out(fout);
	//Inputs
	il.out(fout);
	det.out(fout);
	omega.out(fout);
	sp.out(fout);
	fout.close();

   	// Return results
	
	// Total
	if (spec.y_vec.size() > nout[0])
		std::cerr << "Error: Output for y_vec out of range! ("
		<< nout[0] << " given, needs "
		<< spec.y_vec.size() << ".)" << std::endl;
	else
		for (auto i : spec.y_vec)
			*(y_vec++) = i;
		// *y_vec = -1;

	// Separate
	if (spec.y_sep.size() > (nout[1]+2)*nout[0])
		std::cerr << "Error: Output for y_sep out of range! ("
		<< (nout[1]+2)*nout[0] << " given, needs "
		<< spec.y_sep.size() << ".)" << std::endl;
	else
		for (auto i : spec.y_sep)
			*(y_sep++) = i;
		// *y_sep = -1;


	// XRF
	if (spec.xrf.Z_vec.size() > nout[1])
		std::cerr << "Error: Output for xrf.Z_vec out of range! ("
		<< nout[1] << " given, needs "
		<< spec.xrf.Z_vec.size() << ".)" << std::endl;
	else
	{
		for (auto i : spec.xrf.Z_vec)
			*(Z_vec++) = i;
		// *Z_vec = -1;

		for (auto i : spec.xrf.row)
			*(row++) = i;
		// *row = -1;
	}

	if (spec.xrf.lines.size() > nout[2])
		std::cerr << "Error: Output for xrf.lines out of range! ("
		<< nout[2] << " given, needs "
		<< spec.xrf.lines.size() << ".)" << std::endl;
	else
	{
		for (auto i : spec.xrf.lines)
			*(lines++) = i;
		// *lines = -1;

		for (auto i : spec.xrf.ev_vec)
			*(xrf_ev++) = i;
		// *xrf_ev = -1;

		for (auto i : spec.xrf.y_vec)
			*(xrf_y++) = i;
		// *xrf_y = -1;
	}

	// Compton
	if (spec.comp.ev_vec.size() > nout[3])
		std::cerr << "Error: Output for comp.ev_vec out of range! ("
		<< nout[3] << " given, needs "
		<< spec.comp.ev_vec.size() << ".)" << std::endl;
	else
	{
		for (auto i : spec.comp.ev_vec)
			*(comp_ev++) = i;
		// *comp_ev = -1;

		for (auto i : spec.comp.y_vec)
			*(comp_y++) = i;
		// *comp_y = -1;
	}

	// nout
	nout[0] = spec.y_vec.size();
	nout[1] = spec.xrf.Z_vec.size();
	nout[2] = spec.xrf.lines.size();
	nout[3] = spec.comp.ev_vec.size();
	nout[4] = ds.dose_vec.size();


	// Rayleigh
	*ray_y = spec.ray.y;

	// Detector
	*(det_++) = det.channel.ev_offset;
	*(det_++) = det.channel.ev_gain;
	*n_channels = det.channel.n_channels;
	*(det_++) = det.response.noise;
	*(det_++) = det.response.fano;
	*(det_++) = det.response.gamma;
	*(det_++) = det.response.fs;
	*(det_++) = det.response.ft;
	*(det_++) = det.window.thickness;
	*det_ = det.window.density;
//	for (int i = 0; i < det.window.material.length(); i++)
//		*(win_mat++) = int(det.window.material.at(i));
//	*(win_mat) = int('\n');

    for (auto i: det.window.material)
        *(win_mat++) = i;

	// Illumination
	*(il_++) = il.ev0;
	*(il_++) = il.psi;
	*(il_++) = il.alpha;
	*(il_++) = il.n_photons;
	*(il_) = il.beam_cross_section;

	// Solid angle
	for (int i = 0; i < 4; i++)
		*(sa++) = omega.angle_range[i];
	*(sa++) = omega.theta_inc;
	*(sa) = omega.beta_inc;

	// Radiation dosage
	for (auto i : ds.dose_vec)
		*(dose++) = i;

	return;
}

<file_sep>/Include/input.hpp
/*input.hpp*/
#ifndef INPUT_HPP
#define INPUT_HPP

#include "detector.hpp"
#include "geometry.hpp"
#include "sample.hpp"

void readfile(std::string filename, Sample & sp, Illumination & il, solid_angle & omega, Detector & det);

#endif<file_sep>/python/tools/xraylib.py
__author__ = 'Yue'
import ctypes
from .. import libpath

try:
    lib = ctypes.cdll.LoadLibrary(libpath + '/../xraylib/Lib/linux/libxrl.so')
except OSError:
    lib = ctypes.cdll.LoadLibrary(libpath + '/../xraylib/Lib/windows/libxrl-7.dll')

lib.CS_Total_CP.restype = ctypes.c_double
lib.CS_Total_CP.argtypes = [ctypes.c_char_p, ctypes.c_double]


def CS_Total_CP(CP, kev):
    assert isinstance(CP, str)
    return lib.CS_Total_CP(CP, kev)


def ElementDensity(*args):
    return lib.ElementDensity(*args)


def SymbolToAtomicNumber(*args):
    return lib.SymbolToAtomicNumber(*args)<file_sep>/python/classes/__init__.py
__author__ = 'Yue'
__all__ = ["detector", "geometry", "spectrum"]<file_sep>/Include/math.hpp
/*math.hpp*/
#ifndef	MATH_HPP
#define MATH_HPP

namespace
{
	inline double sq(double x)
	{
		return x*x;
	}
}
#endif<file_sep>/variables.mk
# # # variables used in makefile

# Makefile Paths
mkfile_path := $(shell dirname $(abspath $(lastword $(MAKEFILE_LIST))))
INCLUDE_PATH = $(abspath $(mkfile_path)/Include)
LIB_PATH = $(abspath $(mkfile_path)/Lib)

# XrayLib paths
XRL_PATH = "D:\Program Files\xraylib 64-bit"

# C++11 compiler
g++ = g++

# Others
ifeq ($(OS),Windows_NT)
	#Windows stuff
	# RM = del
	DLIB = dll
	EXE = exe
	XRL_INCLUDE = $(XRL_PATH)/Include
	XRL_LIB = $(XRL_PATH)/Lib
else
	#Linux stuff
	# RM = rm
	DLIB = so
	EXE = out
	XRL_INCLUDE = $(XRL_PATH)/include
	XRL_LIB = $(XRL_PATH)/src/.lib
endif


<file_sep>/src/sample.cpp
//sample.cpp
#include "xraylib.h"
#include "cs.hpp"
#include "constants.hpp"
#include "sample.hpp"


//Compound
Compound::Compound() : Z_vec(_Z_vec), p_vec(_p_vec), molecular_weight(_molecular_weight)
{

}

Compound::Compound(const std::vector<int> & Z_vec_, const std::vector<double> & p_vec_) : Z_vec(_Z_vec), p_vec(_p_vec), molecular_weight(_molecular_weight)
{
	_Z_vec = Z_vec_;
	_p_vec = p_vec_;
	set_mw();
}

Compound::~Compound()
{

}

Compound & Compound::operator=(const Compound & c)
{
	if (this == &c)
		return *this;
	_Z_vec = c.Z_vec;
	_p_vec = c.p_vec;
	set_mw();
	return *this;
}

void Compound::set_mw()
{
	_molecular_weight = 0;
	for (int i=0; i<Z_vec.size(); i++)
		_molecular_weight += p_vec[i]*AtomicWeight(Z_vec[i]);
}

double Compound::mac_tot(const double & ev) const
{
	double mac = 0;
	for (int i=0; i<Z_vec.size(); i++)
		mac += CS_Total(Z_vec[i], ev/1e3)*p_vec[i]*AtomicWeight(Z_vec[i])/molecular_weight;
	return mac;
}

double Compound::dmac_rayleigh_pol(const double & ev, const double & theta, const double & beta) const
{
	double dmac = 0;
	for (int i=0; i<Z_vec.size(); i++)
		dmac += dcs_rayleigh_pol(theta, beta, ev, Z_vec[i])*p_vec[i]*N_A/molecular_weight;
	return dmac;
}

double Compound::dmac_compton_pol(const double & ev, const double & theta, const double & beta) const
{
	double dmac = 0;
	for (int i=0; i<Z_vec.size(); i++)
		dmac += dcs_compton_pol(theta, beta, ev, Z_vec[i])*p_vec[i]*N_A/molecular_weight;
	return dmac;
}

void Compound::show() const
{
	out(std::cout);
}

void Compound::out(std::ostream & ost) const
{
	ost << "# \tZ = ";
	for (auto Z : Z_vec)
		ost << Z << ", ";
	ost << std::endl;
	
	ost << "# \tp = ";
	for (auto p : p_vec)
		ost << p << ", ";
	ost << std::endl;
	
	ost << "# \tMolecular weight = ";
	ost << molecular_weight;
	ost << std::endl;
}

//Monolayer
Monolayer::Monolayer() : Compound::Compound(), density(_density), thickness(_thickness), layer(_layer)
{

}

// Copy constructor: used when passing through function calls
Monolayer::Monolayer(const Monolayer & ml) : Compound::Compound(), density(_density), thickness(_thickness), layer(_layer)
{
	operator=(ml);
}

Monolayer::Monolayer(const std::vector<int> & _Z_vec, const std::vector<double> & _p_vec, const double & density_, const double & thickness_, const double & layer_) : Compound::Compound(_Z_vec, _p_vec), density(_density), thickness(_thickness), layer(_layer)
{
	_density = density_;
	_thickness = thickness_;
	_layer = layer_;
}

Monolayer::~Monolayer()
{
	
}

// overloaded assignment operator
Monolayer & Monolayer::operator=(const Monolayer & ml)
{
	if (this == &ml)
		return *this;
	Compound::operator=(ml);
	_density = ml.density;
	_thickness = ml.thickness;
	_layer = ml.layer;
	return *this;
}


void Monolayer::set_layer(int i)
{
	_layer = i;
}

void Monolayer::show() const
{
	out(std::cout);
}

void Monolayer::out(std::ostream & ost) const
{
	ost << "# Parameters for layer #" << layer << ":" << std::endl;
	Compound::out(ost);
	ost << "# \tDensity = " << density << std::endl;
	ost << "# \tThickness = " << thickness << std::endl;
}

//Sample
Sample::Sample() : layer_vec(_layer_vec), nlayers(_nlayers)
{
	update();
}

Sample::~Sample()
{

}

Sample & Sample::operator=(const Sample & s)
{
	if (this == & s)
		return *this;
	_layer_vec = s.layer_vec;
	_nlayers = s.nlayers;
	return *this;
}

void Sample::update()
{
	for (int i = 0; i < _layer_vec.size(); i++)
		_layer_vec[i].set_layer(i);	
	_nlayers = _layer_vec.size();
}

void Sample::add_layer(const Monolayer & monolayer_)
{
	_layer_vec.push_back(monolayer_);
	if (monolayer_.layer == 0)
		_layer_vec.back().set_layer(nlayers);
	_nlayers++;
}

void Sample::show() const
{
	out(std::cout);
}

void Sample::out(std::ostream & ost) const
{
	ost << "# ======================================== #" << std::endl;
	ost << "# # # # # Sample parameters: # # # # #" << std::endl;
	ost << "# Total number of layers = " << nlayers << std::endl;
	for (auto m : layer_vec)
		m.out(ost);
	ost << "# # # # # End of sample parameters # # # # #" << std::endl;
	ost << "# ======================================== #" << std::endl;
}<file_sep>/src/geometry.cpp
//geometry.cpp
#include <iostream>
#include <cmath>

#include "geometry.hpp"

using std::vector;


//Illumination
Illumination::Illumination(double ev0_, double psi_, double alpha_, double n_photons_, double beam_cross_section_) : ev0(_ev0), psi(_psi), alpha(_alpha), n_photons(_n_photons), beam_cross_section(_beam_cross_section)
{
	_ev0 = ev0_;
	_psi = psi_;
	_alpha = alpha_;
	_n_photons = n_photons_;
	_beam_cross_section = beam_cross_section_;
}
Illumination::~Illumination()
{

}

Illumination & Illumination::operator=(const Illumination & il)
{
	if (this != &il)
	{
		_ev0 = il.ev0;
		_psi = il.psi;
		_alpha = il.alpha;
		_n_photons = il.n_photons;
		_beam_cross_section = il.beam_cross_section;
	}
	return *this;
}

void Illumination::show() const
{
	out(std::cout);
}

void Illumination::out(std::ostream & ost) const
{
	ost << "# ======================================== #" << std::endl;
	ost << "# # # # # Illumination: # # # # #" << std::endl;
	ost << "# ev0 = " << ev0 << std::endl;
	ost << "# psi = " << psi << std::endl;
	ost << "# alpha = " << alpha << std::endl;
	ost << "# Number of incident photons = " << n_photons << std::endl;
	ost << "# Beam cross sectional area = " << beam_cross_section << " cm^2" << std::endl;
	ost << "# # # # # End of illumination # # # # #" << std::endl;
	ost << "# ======================================== #" << std::endl;
	
}


void spherical2cartesian(const double & r, const double & theta, const double & beta, double & x, double & y, double & z) //Convert spherical coordinates (r,theta,beta) to Cartesian coordinates (x,y,z).
{
	using std::sin;
	using std::cos;
	x = r*sin(theta)*cos(beta);
	y = r*sin(theta)*sin(beta);
	z = r*cos(theta);
}

double arctan2(const double & x) //like arctan, but the returned value is in [0,np.pi]
{	
	using std::atan;
	double rad = atan(x);
	if (rad<0)
		rad += Pi;
	return rad;
}

void rotationy(const double & rad, double & x, double & y, double & z)
//rotate the axes about original y axis, angle = rad.
{
	if (rad)
	{
		using std::cos;
		using std::sin;
		double x1;
		x1 = cos(rad)*x - sin(rad)*z;
		z = sin(rad)*x + cos(rad)*z;
		x = x1;
	}
}

void rotationx(const double & rad, double & x, double & y, double & z)
//rotate the axes about original x axis, angle = rad.
{
	if (rad)
	{
		using std::cos;
		using std::sin;
		double y1;
		y1 = cos(rad)*y + sin(rad)*z;
		z = -sin(rad)*y + cos(rad)*z;
		y = y1;
	}
}
double Illumination::psi_prime(const double & theta, const double & beta) const
{
	using std::sqrt;
	double theta0 = Pi/2.+psi;
	double x, y, z;
	spherical2cartesian(1, theta, beta, x, y, z);
	rotationy(theta0, x, y, z);
	rotationx(alpha, x, y, z);
	double theta_prime = arctan2(sqrt(x*x+y*y)/(1.*z));
	return Pi/2.-theta_prime;
}
	
// solid_angle
solid_angle::solid_angle():theta(_theta), beta(_beta), subtend(_subtend), theta_inc(_theta_inc), beta_inc(_beta_inc)//, angle_range(_angle_range)
{
	_angle_range[0] = 0;
	_angle_range[1] = Pi;
	_angle_range[2] = 0;
	_angle_range[3] = Pi/2;
	_theta_inc = Pi/180;
	_beta_inc = Pi/180;
	update();
}
solid_angle::solid_angle(const double *ar, double th_inc, double be_inc):theta(_theta), beta(_beta), subtend(_subtend), theta_inc(_theta_inc), beta_inc(_beta_inc)//, angle_range(_angle_range)
{
	for (int i=0; i<4; i++)
		_angle_range[i] = ar[i];
	_theta_inc = th_inc;
	_beta_inc = be_inc;
	update();
}
solid_angle::~solid_angle()
{

}

solid_angle & solid_angle::operator=(const solid_angle & sa)
{
	if (this == &sa)
		return *this;
	for (int i=0; i<4; i++)
		_angle_range[i] = sa.angle_range[i];
	_theta_inc = sa._theta_inc;
	_beta_inc = sa._beta_inc;
	update();
	return *this;
}

double solid_angle::domega(const double & theta) const
{
	return std::sin(theta)*_theta_inc*_beta_inc;
}

void solid_angle::update()
{
	_theta.clear();
	for (double x = angle_range[0]; x <= angle_range[1]; x += _theta_inc)
		_theta.push_back(x);
	_beta.clear();
	for (double x = angle_range[2]; x <= angle_range[3]; x += _beta_inc)
		_beta.push_back(x);
	_subtend = (std::cos(angle_range[0])-std::cos(angle_range[1]))*(angle_range[3]-angle_range[2]);
}

// const vector<double> & solid_angle::get_theta(double th_inc)
// {
	// if ((th_inc > 0) && (th_inc != theta_inc))
	// {
		// theta_inc = th_inc;
		// _theta.clear();
		// for (double x = angle_range[0]; x <= angle_range[1]; x += theta_inc)
			// _theta.push_back(x);
	// }
	// else if (theta.empty())
		// for (double x = angle_range[2]; x <= angle_range[3]; x += theta_inc)
			// _theta.push_back(x);
	// return theta;
// }



// const vector<double> & solid_angle::get_beta(double be_inc)
// {
	// if ((be_inc > 0) && (be_inc != beta_inc))
	// {
		// beta_inc = be_inc;
		// _beta.clear();
		// for (double x = angle_range[2]; x <= angle_range[3]; x += beta_inc)
			// _beta.push_back(x);
	// }
	// else if (beta.empty())
		// for (double x = angle_range[2]; x <= angle_range[3]; x += beta_inc)
			// _beta.push_back(x);
	// return beta;
// } 

void solid_angle::show() const
{
	out(std::cout);
}

void solid_angle::out(std::ostream & ost) const
{
	ost << "# ======================================== #" << std::endl;
	ost << "# # # # # Solid angle definition # # # # #" << std::endl;
	ost << "# angle range: ";
	for (int i = 0; i < 4; i++)
		ost << angle_range[i] << ' ';
	ost << std::endl;
	ost << "# theta_inc = " << _theta_inc << std::endl;
	ost << "# beta_inc = " << _beta_inc << std::endl;
	ost << "# theta: ";
	for (vector<double>::const_iterator i = theta.begin(); i < theta.end(); i++)
		ost << *i << ' ';
	ost << std::endl;
	ost << "# beta: ";
	for (vector<double>::const_iterator i = beta.begin(); i < beta.end(); i++)
		ost << *i << ' ';
	ost << std::endl;
	ost << "# Total subtended solid angle: " << subtend << std::endl;
	ost << "# # # # # End of solid angle definition # # # # #" << std::endl;
	ost << "# ======================================== #" << std::endl;
}

double atten_mono(const double & ev0, 
		const double & ev, 
		const double & psi, 
		const double & psiprime, 
		const Monolayer & ml)
{
	if ((psiprime < 1e-6) && (psiprime > -1e-6)) // Check if psiprime == 0
		return 0;
	double t = ml.thickness;
	double rho = ml.density;
	if (t == 0 || rho == 0)
		return 1;
	double mac0 = ml.mac_tot(ev0);
	double mac1 = ml.mac_tot(ev);
	double sp0 = std::sin(psi);
	double sp1 = std::sin(psiprime);
	double temp = (mac0/sp0+mac1/sp1)*rho;
	if (mac0+mac1*sp0/sp1 < 1e-50 && mac0+mac1*sp0/sp1 > -1e-50) // To avoid division by zero
		return 0;
	if (psiprime > 0)
		return (1-std::exp(-(mac0/sp0+mac1/sp1)*rho*t))/(mac0+mac1*sp0/sp1);
	else
		return (std::exp(mac1*rho*t/sp1)-std::exp(-mac0*rho*t/sp0))/(mac0+mac1*sp0/sp1);
}

double atten_refl(const double & ev0, 
		const double & ev, 
		const double & psi, 
		const double & psiprime, 
		const Monolayer & ml)
{
	double t = ml.thickness;
	double rho = ml.density;
	if (t == 0 || rho == 0)
		return 1;
	else if (psi == 0 || psiprime == 0)
	    return 0;
	return std::exp(-(ml.mac_tot(ev0)/std::sin(psi)+ml.mac_tot(ev)/std::sin(psiprime))*ml.density*ml.thickness);
}		
		
double atten_trans_in(const double & ev0, 
		const double & psi, 
		const Monolayer & ml)
{
	if (ml.density == 0 || ml.thickness == 0)
		return 1;
	else if (psi == 0)
	    return 0;
	return std::exp(-ml.mac_tot(ev0)/std::sin(psi)*ml.density*ml.thickness);
}		
		
double atten_trans_out(const double & ev, 
		const double & psiprime, 
		const Monolayer & ml)
{
	if (ml.density == 0 || ml.thickness == 0)
		return 1;
	else if (psiprime == 0)
	    return 0;
	return std::exp(ml.mac_tot(ev)/std::sin(psiprime)*ml.density*ml.thickness);
}<file_sep>/read.py
import numpy as np
import matplotlib.pyplot as plt


def read(fname):
    # Read channel values
    a = np.genfromtxt(fname, invalid_raise=False, skiprows=0)

    # Read detector configurations and line labels
    labels = ['Total']
    f = open(fname)
    for line in f.readlines():
        if 'ev_offset' in line:
            ev_offset = float(line.split('=')[1])
        elif 'ev_gain' in line:
            ev_gain = float(line.split('=')[1])
        elif 'n_channels' in line:
            n_channels = float(line.split('=')[1])
        elif ('Z =' in line) and ('(' in line):
            labels.append(line.split('(')[-1].strip().strip(')'))
    f.close()
    labels += ['Rayleigh', 'Compton']

    ev = ev_offset + ev_gain * np.arange(n_channels)

    return ev, a, labels

def plot(ev, a, labels, xlim=None, ylim=None, show=True):
    ax = plt.subplot(111)
    plt.plot(ev / 1e3, a.T)
    plt.yscale('log')
    plt.ylim(ylim)
    plt.xlim(xlim)
    plt.xlabel('E (KeV)')
    plt.ylabel(r'$I(E)/I_0(E_0)$')
    # Shrink current axis by 20%
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    plt.legend(labels, loc='upper right', ncol=1, bbox_to_anchor=(1.35, 1), borderaxespad=0.)
    if show:
        plt.show()

def dose(fname):
    f = open(fname)
    dose = []
    for line in f.readlines():
        if ('Layer' in line) and ('Gy' in line):
            dose.append(float(line.split(':')[-1].split()[0]))
    f.close()
    return np.array(dose)



if __name__ == '__main__':
    import sys
    from python import plt_kwargs

    try:
        fname = sys.argv[1]
    except IndexError:
        fname = 'out.txt'

    # Read data file
    ev, a, labels = read(fname)

    # Plot spectrum
    plt.figure()
    plt.title('Spectrum read from %s' % fname)
    plot(ev, a, labels, xlim=plt_kwargs['xlim'], ylim=plt_kwargs['ylim'])<file_sep>/src/xrf.cpp
/* xrf.cpp */
#include <iostream>
#include <vector>
#include "xraylib.h"

#define MP 1 // OpenMP

using namespace std;

#define K_LINES -29
#define	L_LINES -113
#define M_LINES -219
#define N_LINES -324
#define O_LINES -374
#define P_LINES -383

vector<int> mac_xrf(double ev0, int Z, vector<double> & ev_vec, vector<double> & y_vec, double weight = 1.0, int line_end = P_LINES)
{
	vector<int> lines;
	int line_start;
	if (ev0 >= 1e3*EdgeEnergy(Z, K_SHELL))
		line_start = -1;
	else if (line_end <= L_LINES && ev0 >= 1e3*EdgeEnergy(Z, L3_SHELL))
		line_start = K_LINES-1;
	else if (line_end <= M_LINES && ev0 >= 1e3*EdgeEnergy(Z, M5_SHELL))
		line_start = L_LINES-1;
	else if (line_end <= N_LINES && ev0 >= 1e3*EdgeEnergy(Z, N7_SHELL))
		line_start = M_LINES-1;
	else if (line_end <= O_LINES && ev0 >= 1e3*EdgeEnergy(Z, O7_SHELL))
		line_start = N_LINES-1;	
	else if (line_end <= P_LINES && ev0 >= 1e3*EdgeEnergy(Z, P5_SHELL))
		line_start = N_LINES-1;
	else 
	{
		cout << "XRF lines not found for Z = " << Z << " @ " << ev0 << " eV!" << endl;
		return lines;
	}
	
	#pragma omp parallel for ordered if(MP)
	for (int line = line_start; line >= line_end; line--)
	{
		double ev = LineEnergy(Z, line)*1e3;
		if (ev > 1e-6)
		{
			double y = CS_FluorLine_Kissel_Cascade(Z, line, ev0/1e3);
			if (y > 1e-50)
			{
				#pragma omp ordered
				{
					ev_vec.push_back(ev);
					y_vec.push_back(y*weight);
					lines.push_back(line);
				}
			}
		}
	}
	return lines;
}




<file_sep>/python/tools/savespec.py
#savespec
import h5py
import os

def dim_scale(f,dset):
	ndims = len(dset.dims)
	dset.dims[ndims-2].label = 'theta'
	dset.dims[ndims-1].label = 'beta'
	dset.dims.create_scale(f['omega/theta'])
	dset.dims.create_scale(f['omega/theta_rad'])
	dset.dims.create_scale(f['omega/beta'])
	dset.dims.create_scale(f['omega/beta_rad'])
	dset.dims[ndims-2].attach_scale(f['omega/theta'])
	dset.dims[ndims-2].attach_scale(f['omega/theta_rad'])
	dset.dims[ndims-1].attach_scale(f['omega/beta'])
	dset.dims[ndims-1].attach_scale(f['omega/beta_rad'])

def savespec(rawspec, filename = 'specdata.h5', dir = os.getcwd()):
	if not os.path.exists(dir):
		os.mkdir(dir)
	rawspec_arr = flatten([rawspec])
	try:
		f = h5py.File(dir+filename, 'w')
		
		f['omega/theta'] = rawspec_arr[0].omega.theta_arr
		f['omega/theta_rad'] = rawspec_arr[0].omega.theta_rad_arr
		f['omega/beta'] = rawspec_arr[0].omega.beta_arr
		f['omega/beta_rad'] = rawspec_arr[0].omega.beta_rad_arr
		f['omega/range'] = rawspec_arr[0].omega.range
		f['omega/n_theta'] = rawspec_arr[0].omega.n_theta
		f['omega/n_beta'] = rawspec_arr[0].omega.n_beta
		
		for rawspec in rawspec_arr:
			
			group = f.require_group(rawspec.stype+'/'+rawspec.name)
			# try:
				# group = f[rawspec.stype+'/'+rawspec.name]
			# except KeyError:
				# group = f.create_group(rawspec.stype+'/'+rawspec.name)
			dset_list = rawspec.__dict__.keys()
			dset_list.remove('stype')
			dset_list.remove('name')
			dset_list.remove('omega')
			for dset in dset_list:
				try:
					group[dset] = rawspec.__getattribute__(dset)
					dim_scale(f,group[dset])
				except TypeError:
					pass
	finally:
		f.close()

def readspec(filename = 'specdata.h5', dir = os.getcwd()):
	try:
		f = h5py.File(dir+filename, 'r')
		out = []
		omega = solid_angle(angle_range=f['omega/range'].value,n_theta = f['omega/n_theta'].value,n_beta = f['omega/n_beta'].value)
		stype_list = f.keys()
		stype_list.remove('omega')
		for stype in stype_list[::-1]:
			for name in f[stype]:
				group = f[stype][name]
				out.append(spectrum(ev_mat = group['ev_mat'], intensity_mat = group['intensity_mat'], omega = omega, name = name, stype = stype))
		return out
	finally:
		f.close()<file_sep>/src/makefile
# makefile

include ../variables.mk

SRC_DIR = "."
OBJ_DIR = "../obj"
SRC_FILES = $(wildcard *.cpp)
OBJ_FILES = $(SRC_FILES:.cpp=.o)


all: $(OBJ_FILES) cwrapper.o $(LIB_PATH)
	ar crf $(LIB_PATH)/libsim.a $(OBJ_FILES)
	$(g++) $(OBJ_FILES) cwrapper.o -L$(XRL_LIB) -lxrl -shared -fPIC -o $(LIB_PATH)/libsim.$(DLIB) -fopenmp
	# $(g++) cwrapper.o $(LIB_PATH)/libsim.a -L$(XRL_LIB) -lxrl -shared -fPIC -o $(LIB_PATH)/libsim.$(DLIB) -fopenmp # works for mingw-w64 on windows, not for linux.

%.o: %.cpp
	@echo "============="
	@echo "Compiling $<"
	@$(g++) -c $< -fPIC -o $@ -I$(XRL_INCLUDE) -I$(INCLUDE_PATH) -std=c++11 -fopenmp

cwrapper.o: cwrapper.cc
	@echo "============="
	@echo "Compiling $<"
	@$(g++) -c $< -fPIC -o $@ -I$(XRL_INCLUDE) -I$(INCLUDE_PATH) -std=c++11 -fopenmp

$(LIB_PATH):
	- mkdir $(LIB_PATH)

clean:
	$(RM) -rvf *.o $(LIB_PATH)/*.$(DLIB) $(LIB_PATH)/*.a all<file_sep>/python/classes/detector.py
# detector.py
import numpy as np

try:
    import xraylib as xrl
except ImportError:
    from ..tools import xraylib as xrl


def energy_to_channel(energy, offset=0., gain=10.):
    """
    Convert energy to channel number.
    :param energy:
        The energy in eV.
    :param offset:
        Energy offset in eV. Or the energy of channel 0.
    :param gain:
        Energy gain in eV. The increment of energy in eV from channel to channel.
    :return:
        Channel number.
    """
    return (energy - offset) // gain


class Channel(object):
    """
    Class for detector channels.
    Attributes:
        ev_offset:  Energy offset in eV.
        ev_gain:    Channel gain in eV.
        n_channels: Total number of channels.
    """

    def __init__(self, ev_offset=0, ev_gain=10, n_channels=2048):
        self.ev_offset = ev_offset
        self.ev_gain = ev_gain
        self.n_channels = n_channels
        self.ev_arr = np.arange(ev_offset, ev_offset + ev_gain * n_channels, ev_gain)

    def channel(self, ev):
        """
        Convert energy in eV to channel number.
        :param ev:
            The energy in eV.
        :return:
            Channel number
        """
        return energy_to_channel(ev, self.ev_offset, self.ev_gain)


class Response(object):
    """
    Class for detector response.
    Attributes:
        noise:
        fano:
        fs:
        ft:
        ev_gain:
    """

    def __init__(self,
                 noise=100,
                 fano=0.114,
                 gamma=2.5,
                 fs=0.03,
                 ft=0.02,
                 ev_gain=None):
        self.noise = noise
        self.fano = fano
        self.gamma = gamma
        self.fs = fs
        self.ft = ft
        self.ev_gain = ev_gain

    def FWHM(self, ev, **kwargs):
        """
        FWHM function (of energy in eV).
        :rtype : float
        :param ev: Energy in eV.
        :param kwargs:
        :return:
            FWHM in eV.
        """
        noise = kwargs.get('noise', self.noise)
        fano = kwargs.get('fano', self.fano)
        sigma = np.sqrt((noise / 2.3548) ** 2 + 3.58 * fano * ev)
        return 2.3548 * sigma


class Window(object):
    """
    Class for detector window.
    :Attributes:
        material:
            String. Material of the filtering window on the detector.
        thickness:
            Thickness of the filtering window material. In cm.
        density:
            Mass density of the filtering window material. In g/cm^3.
    """

    def __init__(self,
                 material='Be',
                 thickness=24e-4,
                 density=None):
        self.material = material
        self.thickness = thickness
        if density is None:
            density = xrl.ElementDensity(xrl.SymbolToAtomicNumber(material))
            if not density:
                density = 1
        self.density = density

    def transmission(self, ev):
        """
        The transmission function of the window.
        :param ev:
            The energy in eV of the incident beam.
        :return:
            The ratio of the intensity of the transmitted beam through the window vs that of the incident.
        """
        _mac = xrl.CS_Total_CP(self.material, ev / 1000.)
        # _mac = compound(CP=self.material).mac_total(ev)
        return np.exp(-_mac * self.density * self.thickness)


class Detector(object):
    """
    Class for detector.
    :Attributes:
        channel: A detector channel object.
        response:   A detector response object.
        window: A detector window object.
    """

    def __init__(self, channel=Channel(), response=Response(), window=Window()):
        self.channel = channel
        self.response = response
        response.ev_gain = channel.ev_gain
        self.window = window<file_sep>/python/tools/elementlookup.py
from . import __path__ as datapath

datapath = datapath[0]

def lookup(x):
    """

    :param x:
    :return:
    """
    f = open(datapath + '/periodictable.txt', 'r')  # open periodic table file
    data = f.readlines()
    f.close()

    # dictionaries
    d = {}
    ivd1 = {}
    ivd2 = {}

    # create dictionaries from data in the file
    for line in data:
        Z, symbol, name = (s.strip() for s in line.split('-'))
        Z = int(Z)
        d[Z] = symbol, name
        ivd1[symbol] = Z, name
        ivd2[name] = Z, symbol

    # lookup in the dictionaries
    try:
        x = x.title()
        try:
            return ivd1[x]
        except KeyError:
            return ivd2[x]
    except AttributeError:
        return d[x]


def symbol2number(Z):
    if type(Z) == str:
        try:
            Z = int(Z)
        except ValueError:
            return lookup(Z)[0]
    return Z


def number2symbol(Z):
    if type(Z) == str:
        try:
            Z = int(Z)
        except ValueError:
            return Z
    return lookup(Z)[0]
<file_sep>/makefile
# makefile

include variables.mk

all: main.o
	$(MAKE) -C ./src all
	$(g++) main.o -L./Lib -lsim -L$(XRL_LIB) -lxrl -o main.$(EXE) -fopenmp #-static
	echo export LD_LIBRARY_PATH="$(LIB_PATH)" > init.sh
	echo export OMP_NUM_THREADS=16 >> init.sh
	echo export PYTHONPATH="$(mkfile_path):\$$PYTHONPATH" >> init.sh

main.o: main.cpp
	$(g++) -c main.cpp -I$(INCLUDE_PATH) -I$(XRL_INCLUDE) -std=c++11 -fopenmp

clean:
	$(MAKE) -C ./src clean
	$(RM) -rvf *.o main.exe $(LIB_PATH)/*.a $(LIB_PATH)/*.$(DLIB) $(LIB_PATH) all
<file_sep>/main.cpp
#include <iostream>
#include <fstream>
#include <string>

#include "spectrum.hpp"
#include "dose.hpp"
#include "input.hpp"

int main(int argc, char ** argv)
{
	std::string input_file = "input.txt", output_file = "out.txt";
	if (argc > 2)
	{
		output_file = argv[2];
		input_file = argv[1];
	}
	else if (argc > 1)
		input_file = argv[1];

	Sample sp;
	Illumination il;
	solid_angle omega;
	Detector det;
	
	// Read input file
	std::cout << "Reading file \"" << input_file << "\"..." << std::endl;
	readfile(input_file, sp, il, omega, det);
	
	// Show configurations
	sp.show();
	omega.show();
	det.show();

	// Calculate spectrum
	std::cout << "Calculating spectrum..." << std::endl;
	Spectrum spec(sp, il, omega, det, true, true, true);

	// Calculate dose
	Dose ds(sp, il);

	// Show results on screen
	spec.show();

	// Save results to file
	std::cout << "Writing output into \"" << output_file << "\"..." << std::endl;
	static std::ofstream fout;
	fout.open(output_file);
	// Outputs
	spec.out(fout);
	fout << std::endl;
	ds.out(fout);
	//Inputs
	il.out(fout);
	det.out(fout);
	omega.out(fout);
	sp.out(fout);
	std::cout << "Done!" << std::endl;

	return 0;
}<file_sep>/Include/sample.hpp
//sample.hpp
#ifndef SAMPLE_HPP
#define SAMPLE_HPP

#include <vector>
#include <iostream>

class Compound
{
private:
	std::vector<int> _Z_vec;
	std::vector<double> _p_vec;
	double _molecular_weight;
	void set_mw();
public:
	Compound();
	Compound(const std::vector<int> & Z_vec_, const std::vector<double> & p_vec_);
	~Compound();
	Compound & operator=(const Compound & c);
	const std::vector<int> & Z_vec;
	const std::vector<double> & p_vec;
	const double & molecular_weight;
	double mac_tot(const double & ev) const;
	double dmac_rayleigh_pol(const double & ev, const double & theta, const double & beta) const;
	double dmac_compton_pol(const double & ev, const double & theta, const double & beta) const;
	void show() const;
	void out(std::ostream & ost) const;
};

class Monolayer : public Compound
{
private:
	double _density;
	double _thickness;
	double _layer;
public:
	Monolayer();
	Monolayer(const Monolayer & ml);
	Monolayer(const std::vector<int> & _Z_vec, const std::vector<double> & _p_vec, const double & density_, const double & thickness_, const double & layer_ = 0);
	~Monolayer();
	Monolayer & operator=(const Monolayer & ml);
	const double & density;
	const double & thickness;
	const double & layer;
	void set_layer(int i);
	void show() const;
	void out(std::ostream & ost) const;
};

class Sample
{
private:
	std::vector<Monolayer> _layer_vec;
	int _nlayers;
	void update();
public:
	Sample();//const std::vector<Monolayer> & layer_vec_);
	~Sample();
	Sample & operator=(const Sample & s);
	const std::vector<Monolayer> & layer_vec;
	const int & nlayers;
	void add_layer(const Monolayer & monolayer_);
	void show() const;
	void out(std::ostream & ost) const;
};

#endif<file_sep>/python/classes/spectrum.py
# spectrum.py
import matplotlib.pyplot as plt

from .detector import *
from .geometry import *

# from matplotlib.colors import LogNorm
# from matplotlib.ticker import LogFormatterMathtext

# # # #fonts# # #
# import matplotlib
# from matplotlib import rc

# matplotlib.rcParams['pdf.fonttype'] = 'truetype'
# fontProperties = {'family':'serif','serif':['Arial'],
# 'weight' : 'normal', 'size' : '12'}
# rc('font',**fontProperties)
# # # #


class Xrf(object):
    def __init__(self, y_mat=None,
                 y_vec=None, ev_vec=None, lines=None, Z_vec=None, row=None):
        """


        :param y_mat:
        :param y_vec:
        :param ev_vec:
        :param lines:
        :param Z_vec:
        :param row:
        """
        self.y_mat = y_mat
        self.y_vec = y_vec
        self.ev_vec = ev_vec
        self.lines = lines
        self.Z_vec = Z_vec
        self.row = row


class Rayleigh(object):
    def __init__(self, y_mat=None,
                 y=None, ev0=None):
        self.y_mat = y_mat
        self.y = y
        self.ev0 = ev0


class Compton(object):
    def __init__(self, y_mat=None,
                 y_vec=None, ev_vec=None):
        self.y_mat = y_mat
        self.y_vec = y_vec
        self.ev_vec = ev_vec


class Spectrum():
    def __init__(self,
                 y_vec=None,
                 y_sep=None,
                 labels=[],
                 xrf=Xrf(),
                 ray=Rayleigh(),
                 comp=Compton(),
                 det=Detector(),
                 il=Illumination(),
                 omega=SolidAngle()):
        self.y_vec = y_vec
        self.y_sep = y_sep
        self.labels = labels
        self.xrf = xrf
        self.rayleigh = ray
        self.compton = comp
        self.detector = det
        self.illumination = il
        self.omega = omega

    def show(self, **kwargs):
        xlim = kwargs.get('xlim')
        show_total = kwargs.get('show_total', True)
        logy = kwargs.get('logy', True)
        show = kwargs.get('show', True)
        yshift = kwargs.get('yshift', 0)
        ax = plt.subplot(111)
        plt.plot(self.detector.channel.ev_arr / 1.e3, self.y_sep.T + yshift)
        labels = self.labels
        if show_total:
            plt.plot(self.detector.channel.ev_arr / 1.e3, self.y_vec + yshift, linestyle='--')
            # plt.fill_between(x.ev_arr[845:880]/1.e3,1e-18,total[845:880],color = 'black',alpha = 0.15)
            labels += ['Total']
        ymax = self.y_vec.max()
        ylim = kwargs.get('ylim', [1e-18, 5 * ymax])


        # Shrink current axis by 20%
        box = ax.get_position()
        ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
        # plt.legend(loc=0,ncol=3)
        plt.legend(self.labels, ncol=1, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
        if logy:
            plt.yscale('log')
        plt.xlim(xlim)
        plt.ylim(ylim)
        plt.xlabel('E (KeV)')
        plt.ylabel(r'$I(E)/I_0(E_0)$')
        if show:
            plt.show()

    def save_npz(self, filename='spec.npz'):
        np.savez(filename, ev_arr=self.detector.channel.ev_arr, y_vec=self.y_vec, y_sep=self.y_sep)

    def save_hdf(self, filename='spec.hdf5'):
        pass
<file_sep>/python/pyapi.py
# pyapi.py
import ctypes
import os

from tools.elementlookup import number2symbol  # Must be imported before loading the DLL.

from . import libpath


try:
    lib = ctypes.cdll.LoadLibrary(libpath + '/libsim.so')
except OSError:
    lib = ctypes.cdll.LoadLibrary(libpath + '/libsim.dll')

from .classes.spectrum import *
from . import config

lib.sim.restype = None
lib.sim.argtypes = [ctypes.c_char_p,
                    ctypes.c_char_p,
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_int),
                    ctypes.POINTER(ctypes.c_int),
                    ctypes.POINTER(ctypes.c_int),
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_int),
                    ctypes.c_char_p,
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_double),
                    ctypes.POINTER(ctypes.c_int)]


def sim(input_file,
        output_file,
        y_vec,
        y_sep,
        Z_vec,
        row,
        lines,
        xrf_ev,
        xrf_y,
        comp_ev,
        comp_y,
        ray_y,
        det,
        n_channels,
        win_mat,
        il,
        sa,
        dose,
        nout):
    if not os.path.isfile(input_file):
        raise IOError("File %s does not exit!" % input_file)

    lib.sim(ctypes.c_char_p(input_file),
            ctypes.c_char_p(output_file),
            y_vec.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            y_sep.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            Z_vec.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
            row.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
            lines.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
            xrf_ev.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            xrf_y.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            comp_ev.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            comp_y.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            ray_y.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            det.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            n_channels.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
            win_mat,
            il.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            sa.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            dose.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            nout.ctypes.data_as(ctypes.POINTER(ctypes.c_int)))


def calc(input_file="input.txt",
         output_file="output.txt",
         nout=config['nout'], return_dose=False):  # N of channels, N of Z, N of lines, N of thetas, N of layers
    y_vec = np.zeros(nout[0])
    y_sep = np.zeros(nout[0] * (nout[1] + 2))
    Z_vec = np.zeros(nout[1], dtype=np.dtype(ctypes.c_int))
    row = np.zeros(nout[1] + 2, dtype=np.dtype(ctypes.c_int))
    lines = np.zeros(nout[2], dtype=np.dtype(ctypes.c_int))
    xrf_ev = np.zeros(nout[2])
    xrf_y = np.zeros(nout[2])
    comp_ev = np.zeros(nout[3])
    comp_y = np.zeros(nout[3])
    ray_y = np.zeros(1)
    det = np.zeros(9)
    n_channels = np.zeros(1, dtype=np.dtype(ctypes.c_int))
    win_mat = ctypes.create_string_buffer(20)
    il = np.zeros(5)
    sa = np.zeros(6)
    dose = np.zeros(nout[4])
    nout = np.array(nout, dtype=np.dtype(ctypes.c_int))

    _nout = nout.copy()

    sim(input_file,
        output_file,
        y_vec,
        y_sep,
        Z_vec,
        row,
        lines,
        xrf_ev,
        xrf_y,
        comp_ev,
        comp_y,
        ray_y,
        det,
        n_channels,
        win_mat,
        il,
        sa,
        dose,
        nout)

    # check out-of-range errors
    if (nout > _nout).any():
        raise IOError("Output out of range! \nGiven: " + str(_nout) + "\nNeeds: " + str(nout))

    xrf = Xrf(None, xrf_y[:nout[2]], xrf_ev[:nout[2]], lines[:nout[2]], Z_vec[:nout[1]], row[:nout[1] + 1])
    comp = Compton(None, comp_y[:nout[3]], comp_ev[:nout[3]])
    ray = Rayleigh(None, ray_y, il[0])
    _det = Detector(Channel(*det[:2], n_channels=n_channels),
                    Response(*det[2:-2], ev_gain=det[1]),
                    Window(win_mat, *det[-2:]))
    _il = Illumination(*il)
    omega = SolidAngle(sa[:4], *sa[-2:])

    # labels = Z_vec[:nout[1]].tolist()
    labels = [number2symbol(Z) for Z in Z_vec[:nout[1]]]
    labels.append('Rayleigh')
    labels.append('Compton')

    if return_dose:
        return Spectrum(y_vec[:n_channels], y_sep[:n_channels * (nout[1] + 2)].reshape(-1, n_channels), labels, xrf,
                        ray, comp, _det, _il, omega), dose[:nout[5]]
    else:
        return Spectrum(y_vec[:n_channels], y_sep[:n_channels * (nout[1] + 2)].reshape(-1, n_channels), labels, xrf,
                        ray, comp, _det, _il, omega)


        # print y_vec
        # print Z_vec
        # print row
        # print lines
        # print xrf_ev
        # print xrf_y
        # print comp_ev
        # print comp_y
        # print ray_y
        # print det
        # print n_channels
        # print win_mat
        # print il
        # print sa
        # print nout

        # spec = calc()

        # print spec.detector.window.material
        # print spec.labels
        # print spec.y_sep[0]<file_sep>/python/classes/geometry.py
# geometry.py
import numpy as np


class SolidAngle(object):
    def __init__(self, angle_range=[0, np.pi, -np.pi / 2, np.pi / 2], theta_inc=np.pi / 180, beta_inc=np.pi / 180,
                 unit="rad"):
        """

        :param angle_range: 
        :param theta_inc: 
        :param beta_inc: 
        :param unit: 
        """
        if unit == "deg":
            angle_range = np.radians(angle_range)
            theta_inc = np.radians(theta_inc)
            beta_inc = np.radians(beta_inc)
        self.range = angle_range
        self.theta_inc = theta_inc
        self.beta_inc = beta_inc
        self.theta_arr = np.arange(angle_range[0], angle_range[1], self.theta_inc)
        self.beta_arr = np.arange(angle_range[2], angle_range[3], self.beta_inc)
        self.subtend_calc = np.sum(np.sin(self.theta_arr) * self.theta_inc * self.beta_inc) * self.beta_arr.size
        self.subtend = (np.cos(self.theta_arr[0]) - np.cos(self.theta_arr[-1])) * (self.beta_arr[-1] - self.beta_arr[0])

    def shift(self, new_center):
        pass

    def slice(self, angle_range=[0, np.pi, -np.pi / 2, np.pi / 2]):
        """


        :rtype : 2D array.
        :param angle_range: 
        :return: 
        """
        theta_slice = np.logical_and(self.theta_arr >= angle_range[0], self.theta_arr <= angle_range[1])
        beta_slice = np.logical_and(self.beta_arr >= angle_range[2], self.beta_arr <= angle_range[3])
        omega_slice = np.multiply(*np.meshgrid(theta_slice, beta_slice))
        return omega_slice

    def domega(self, theta_rad):
        return (np.cos(theta_rad - self.theta_inc / 2.) - np.cos(theta_rad + self.theta_inc / 2.)) * self.beta_inc


class Illumination(object):
    def __init__(self, ev0=1e4, psi=0, alpha=0, n_photons=1, beam_cross_section=1e-12*np.pi):
        self.ev0 = ev0
        self.psi = psi
        self.alpha = alpha
        self.n_photons = n_photons
        self.beam_cross_section = beam_cross_section
<file_sep>/Include/constants.hpp
//constants.hpp
#ifndef CONSTANTS_HPP
#define CONSTANTS_HPP
namespace
{
	const double r_e = 2.8179403267e-13; //classical radius of electron in cm
	const double m_e = 0.510998e6; //electron rest mass in eV/c^2
	const double Pi  = 3.1415926535897932384626433832795;
	const double N_A = 6.0221413e+23; //Avogadro constant
	const double eV_in_Joules = 1.60217657e-19; // 1 eV in Joules
}
#endif<file_sep>/README.md
sim-xrf
=======
This is a python/C++ program for simulating X-Ray Fluorescence (XRF) spectrum acquired at a synchrotron radiation beamline.


1. Installation
=======
The source files need to be compiled with C++11 (tested with g++ version 4.7 and mingw-w64 x86_64-4.9.1-posix-seh-rt_v3-rev1). 

Xraylib with python bindings is also needed at compiling.

Xraylib is available at https://github.com/tschoonj/xraylib.

1.1. Compiling Xraylib
-------
Instructions on compiling XrayLib can be found at https://github.com/tschoonj/xraylib/wiki.

After compilation of XrayLib, change the "XRL_PATH" in "variables.mk" to the actual installation directory of XrayLib.

1.2. Compiling sim-xrf
-------
In the top directory do "make". If your C++ compiler is not called with "g++", modify the variable "g++" in "variables.mk" according to your actual compiler.

1.3. Initialization
-------
Initialize the program by running ". init.sh" before running the program.


2. Running sim-xrf
=======
Sim-xrf can be run in C++ mode and Python mode. In both modes, an input file defining the parameters needed in the calculation is required.
The output of the C++ mode is a text file; the Python mode can in addition plot the spectrum or read the output files and plot them. 

2.1. Input file
-------
The input file is a text file defining the parameters needed in the calculation. See the example "input.txt" file for instructions.

2.2. Running in C++ mode
-------
On linux:

    ./main.out [input file name] [output file name]

On windows:

    main.exe [input file name] [output file name]

2.3. Running in Python mode
-------
### 2.3.1. Direct run in command line

####  Calculation + plotting:


    python simpy.py [input file name] [output file name]
    
#### Reading output file + plotting:


    python read.py [input file name] [output file name]
    
### 2.3.2. Import as a module:

#### Calculate spectrum


    from python.pyapi import calc
    
    nout = 3000, 30, 500, 500  # N of channels, N of Z, N of lines, N of thetas
    xlim = 0, 11
    ylim = 1e-14, 1e-4
    
    spec = calc(input_file="input.txt", output_file="output.txt")
    spec.show(xlim=xlim, ylim=ylim)

#### Read/plot output files


    from read import read, plot
        
    fname = "output.txt"
    xlim = 0, 11
    ylim = 1e-14, 1e-4
    
    # Read data file
    ev, a, labels = read(fname)
    
    # Plot spectrum
    plt.figure()
    plt.title('Spectrum read from %s' % fname)
    plot(ev, a, labels, xlim=xlim, ylim=ylim)<file_sep>/__init__.py
__author__ = 'Yue'
<file_sep>/Include/cs.hpp
#ifndef CS_HPP
#define CS_HPP

double sq(double x);

double ev2nm(double ev); //Convert photon energy (in ev) to wavelength (in nm).

double theta2x(double theta_rad, double ev); //theta in rad -> x in cm^-1


/*differential cross section for single electron*/

double thomson_pol(double theta_rad, double beta_rad);

double thomson_unpol(double theta_rad);

double dcs_rayleigh_pol(double theta_rad, double beta_rad, double ev, int Z);
	
double dcs_rayleigh_unpol(double theta_rad, double ev, int Z);

	
/* Compton scattering cross sections */
	
/* compton dcs of single electron */
//dcs

double ev_scattered(double theta_rad, double ev); //return the energy in ev of the Compton scattered photons.

double klein_nishina_unpol(double theta_rad, double ev); //Differential Klein-Nishina cross section in cm^2 for unpolarized radiation.

double klein_nishina_pol(double theta_rad, double beta_rad, double ev); //Differential Klein-Nishina cross section in cm^2 for polarized radiation.

//total cs
double klein_nishina_total_col(double ev); // Total Klein-Nishina collision cross section in cm^2.

double klein_nishina_total_sca(double ev); //Total Klein-Nishina scattering cross section in cm^2.


//compton dcs of atom
double dcs_compton_pol(double theta_rad, double beta_rad, double ev, int Z); //Differential Compton cross section in cm^2 for polarized radiation.

double dcs_compton_unpol(double theta_rad, double ev, int Z); //Differential Compton cross section in cm^2 for unpolarized radiation.

#endif<file_sep>/Include/spectrum.hpp
//spectrum.hpp
#ifndef SPECTRUM_HPP
#define SPECTRUM_HPP

#include <vector>

#include "sample.hpp"
#include "geometry.hpp"
#include "detector.hpp"

class Xrf
{
private:
	std::vector<double> _y_mat;
	solid_angle _omega;
	std::vector<double> _ev_vec;
	std::vector<double> _y_vec;
	std::vector<int> _lines;
	std::vector<int> _Z_vec;
	std::vector<int> _row;
public:
	Xrf(); //default constructor
	Xrf(double ev0, const Compound & c, double n_photons=1);
	Xrf(double ev0, const Compound & c, const solid_angle omega, double n_photons=1);
	Xrf(const Sample & s, const Illumination & il, const solid_angle & omega_);
	// ~Xrf();
	Xrf & operator=(const Xrf & x);
	const std::vector<double> & y_mat;
	const solid_angle & omega;
	const std::vector<double> & ev_vec;
	const std::vector<double> & y_vec;
	const std::vector<int> & lines;
	const std::vector<int> & Z_vec;
	const std::vector<int> & row;
	// Xrf & operator+(const Xrf & x) const;
	void add(const Xrf & x, bool mat_only = false);
	void sum();
	void show() const;
	void out(std::ostream & ost) const;
};

class Rayleigh
{
private:
	std::vector<double> _y_mat;
	solid_angle _omega;
	double _y;
	double _ev0;
public:
	Rayleigh(); //default constructor
	Rayleigh(double ev0_, const Compound & c, const solid_angle & omega_, double n_photons=1, bool calc_sum = true);
	Rayleigh(const Sample & s, const Illumination & il, const solid_angle & omega_);
	// ~Rayleigh();
	Rayleigh & operator=(const Rayleigh & r);
	const std::vector<double> & y_mat;
	const solid_angle & omega;
	const double & y;
	const double & ev0;
	// Rayleigh & operator+(const Rayleigh & r) const;
	void add(const Rayleigh & x, bool mat_only = false);
	void sum();
	void show() const;
	void out(std::ostream & ost) const;
};

class Compton
{
private:
	double _ev0;
	std::vector<double> _y_mat;
	solid_angle _omega;
	std::vector<double> _ev_vec;
	std::vector<double> _y_vec;
public:
	Compton(); //default constructor
	Compton(double ev0_, const Compound & c, const solid_angle & omega_, double n_photons=1, bool calc_ev = true, bool calc_sum = true);
	Compton(const Sample & s, const Illumination & il, const solid_angle & omega_);
	// ~Compton();
	Compton & operator=(const Compton & c);
	const std::vector<double> & y_mat;
	const solid_angle & omega;
	const std::vector<double> & ev_vec;
	const std::vector<double> & y_vec;
	// Compton & operator+(const Compton & c) const;
	void add(const Compton & x, bool mat_only = false);
	const double & ev0;
	void sum();
	void ev();
	void sum_ev();
	void show() const;
	void out(std::ostream & ost) const;
};

class Spectrum
{
private:
	std::vector<double> _y_vec;
	std::vector<double> _y_sep;
	// std::vector<int> _Z_vec;
	// std::vector<int> _lines;
	// std::vector<int> _row;
	Xrf _xrf;
	Rayleigh _ray;
	Compton _comp;
	solid_angle _omega;
	Sample _sample;
	Illumination _illumination;
	Detector _detector;
public:
	Spectrum();
	Spectrum(const Sample & s, const Illumination & il, const solid_angle & omega, const Detector & det, bool det_response = true, bool det_window = true, bool separate = false);
	Spectrum(double ev0, const Compound & c, const solid_angle & omega, const Detector & det, bool det_response = true, bool det_window = true, bool separate = false);
	// ~Spectrum();
	Spectrum & operator=(const Spectrum & s);	
	
	const std::vector<double> & y_vec;
	const std::vector<double> & y_sep;
	const std::vector<int> & Z_vec;
	const std::vector<int> & lines;
	const std::vector<int> & row;
	const Xrf & xrf;
	const Rayleigh & ray;
	const Compton & comp;
	const solid_angle & omega;
	
	const Sample & sample;
	const Illumination & illumination;
	const Detector & detector;
	
	void set_geom(const Sample & s, const Illumination & il, const solid_angle & omega);
	void set_det(const Detector & det);
	
	void calc(bool self_abs = true);
	
	void genspec(bool det_response = true);
	// void genspec_raw(std::vector<double> & ev_vec_, std::vector<double> & y_vec_);	
	void show() const;
	void out(std::ostream & ost) const;
};

	
	
#endif<file_sep>/Include/detector.hpp
/*detector.hpp*/
#ifndef DETECTOR_HPP
#define DETECTOR_HPP

#include<vector>
#include<string>


class Channel
{
private:
	double _ev_offset;
	double _ev_gain;
	int _n_channels;
public:
	Channel(double ev_offset_ = 0, double ev_gain_ = 10, int n_channels_ = 2048);
	// ~Channel();
	Channel & operator=(const Channel & ch);
	const double & ev_offset;
	const double & ev_gain;
	const int & n_channels;
	int ev_to_channel(double ev) const;
	void bin(const std::vector<double> & ev_raw, const std::vector<double> & y_raw, std::vector<double> & y_binned) const;
	void bin(const double & ev_raw, const double & y_raw, std::vector<double> & y_binned) const;
	void bin(const std::vector<double> & ev_raw, const std::vector<double> & y_raw, std::vector<double> & y_binned, std::vector<double> & y_separate, const std::vector<int> & row) const;
	void bin(const double & ev_raw, const double & y_raw, std::vector<double> & y_binned, std::vector<double> & y_separate) const;
	void show() const;
    void out(std::ostream & ost) const;
};

class Response
{
private:
	double _noise;
	double _fano;
	double _gamma;
	double _fs;
	double _ft;
	double _ev_gain;
public:
	Response(double noise_ = 100,
		double fano_ = 0.114,
		double gamma_ = 2.5,
		double fs_ = 0.03,
		double ft_ = 0.02,
		double ev_gain_ = 0);
	// ~Response();
	Response & operator=(const Response & r);
	const double & noise;
	const double & fano;
	const double & gamma;
	const double & fs;
	const double & ft;
	const double & ev_gain;
	
	void set_gain(const double ev_gain);
	double FWHM(double ev);

	void show() const;
    void out(std::ostream & ost) const;
};

		
class Window
{
private:
	std::string _material;
	double _thickness;
	double _density;
public:
	Window(std::string material_ = "Be",
	double thickness_ = 24e-4,
	double density_ = 1.848);
	// ~Window();
	Window & operator=(const Window & w);
	const std::string & material;
	const double & thickness;
	const double & density;
	double transmission(double ev) const;
	void show() const;
    void out(std::ostream & ost) const;
};
		
class Detector
{
private:
	Channel _channel;
	Response _response;
	Window _window;
public:
	Detector(Channel channel_ = Channel(), Response response_ = Response(), Window window_ = Window());
	// ~Detector();
	Detector & operator=(const Detector & d);
	const Channel & channel;
	const Response & response;
	const Window & window;
	void genspec(const std::vector<double> & ev_raw, const std::vector<double> & y_raw, std::vector<double> & y_binned, bool det_response = true, bool det_window = true) const;
	void genspec(const double & ev_raw, const double & y_raw, std::vector<double> & y_binned, bool det_response = true, bool det_window = true) const;
	void genspec(const std::vector<double> & ev_raw, const std::vector<double> & y_raw, std::vector<double> & y_binned, std::vector<double> & y_separate, const std::vector<int> & row, bool det_response = true, bool det_window = true) const;
	void genspec(const double & ev_raw, const double & y_raw, std::vector<double> & y_binned, std::vector<double> & y_separate, bool det_response = true, bool det_window = true) const;
	void show() const;
    void out(std::ostream & ost) const;
};

#endif<file_sep>/python/tools/__init__.py
__author__ = 'Yue'
__all__ = ["elementlookup"]<file_sep>/python/snr/snip.py
from .pysnip import *<file_sep>/src/spectrum.cpp
/*spectrum.cpp*/

#include "spectrum.hpp"

#include <iostream>
#include <numeric> //std::accumulate
#include <algorithm> //std::find

#include "xraylib.h"
#include "xrf.hpp"
#include "cs.hpp"

#define MP 1

Xrf::Xrf() : ev_vec(_ev_vec), y_vec(_y_vec), lines(_lines), Z_vec(_Z_vec), row(_row), y_mat(_y_mat), omega(_omega) {}

Xrf::Xrf(double ev0, const Compound & c, double n_photons)
	: ev_vec(_ev_vec), y_vec(_y_vec), lines(_lines), Z_vec(_Z_vec), row(_row), y_mat(_y_mat), omega(_omega)
{
	std::vector<int> lines_temp;
	double weight;
	int Z;
	if (n_photons < 0)
		n_photons = 1;
	_row.push_back(0);
	for (int i = 0; i < c.Z_vec.size(); i++)
	{
		Z = c.Z_vec[i];
		weight = c.p_vec[i]*AtomicWeight(Z)/c.molecular_weight*n_photons;
		lines_temp = mac_xrf(ev0, Z, _ev_vec, _y_vec, weight);
		if (lines_temp.size() > 0)
		{
			_Z_vec.push_back(Z);
			_lines.insert(_lines.end(),lines_temp.begin(),lines_temp.end());
			_row.push_back(lines.size());
		}
	}
}

//c++98:
// Xrf::Xrf(double ev0, const Compound & c, const solid_angle omega, n_photons)
	// : ev_vec(_ev_vec), y_vec(_y_vec), lines(_lines), Z_vec(_Z_vec), row(_row), y_mat(_y_mat), omega(_omega)
// {
	// *this = Xrf(ev0, c, n_photons);
	// for (int i = 0; i < y_vec.size(); i++)
		// _y_vec[i] *= omega.subtend/(4*Pi);
// }

// c++11:
Xrf::Xrf(double ev0, const Compound & c, const solid_angle omega, double n_photons)
	: Xrf(ev0, c, n_photons)
{
	#pragma omp parallel for if(MP)
	for (int i = 0; i < y_vec.size(); i++)
		_y_vec[i] *= omega.subtend/(4*Pi);
}

Xrf::Xrf(const Sample & s, const Illumination & il, const solid_angle & omega_)
	: ev_vec(_ev_vec), y_vec(_y_vec), lines(_lines), Z_vec(_Z_vec), row(_row), y_mat(_y_mat), omega(_omega)
{
	_omega = omega_;
	_row.push_back(0);
	int n_pixels = omega.theta.size()*omega.beta.size();
	for (std::vector<Monolayer>::const_iterator layer = s.layer_vec.begin(); layer < s.layer_vec.end(); layer++)
	{
		std::cout << "Calculating XRF for layer " << (*layer).layer << std::endl;
		Xrf *temp_p = new Xrf(il.ev0, *layer, il.n_photons), & temp = *temp_p;
		std::cout << "Calculating XRF attenuation for layer " << (*layer).layer << std::endl;
		temp._y_mat.resize(temp.lines.size()*n_pixels,0);
		#pragma omp parallel for collapse(3) if(MP)
		for (int line=0; line < temp.lines.size(); line++)
		{
			for (int theta = 0; theta < omega.theta.size(); theta++)
			{
				for (int beta = 0; beta < omega.beta.size(); beta++)
				{
					double psiprime = il.psi_prime(omega.theta[theta], omega.beta[beta]);
					int i = line*n_pixels+theta*omega.beta.size()+beta; // the matrix pixel identifier in y_mat
					//attenuation within the layer
					temp._y_mat[i] = temp.y_vec[line]*atten_mono(il.ev0, temp.ev_vec[line], il.psi, psiprime, *layer)*omega.domega(omega.theta[theta])/(4*Pi);
					if (psiprime > 0) //reflection geometry
						//attenuation from layers upstream
						for (std::vector<Monolayer>::const_iterator layer_up = s.layer_vec.begin(); layer_up < layer; layer_up++)
							temp._y_mat[i] *= atten_refl(il.ev0, temp.ev_vec[line], il.psi, psiprime, *layer_up);
					else if (psiprime < 0) //transmission geometry
					{
						//attenuation from layers upstream
						for (std::vector<Monolayer>::const_iterator layer_up = s.layer_vec.begin(); layer_up < layer; layer_up++)
							temp._y_mat[i] *= atten_trans_in(il.ev0, il.psi, *layer_up);
						//attenuation from layers downstream
						for (std::vector<Monolayer>::const_iterator layer_down = layer+1; layer_down < s.layer_vec.end(); layer_down++)
							temp._y_mat[i] *= atten_trans_out(temp.ev_vec[line], psiprime, *layer_down);
					}
				}
			}
		}
		add(temp, true);
		delete temp_p;
	}
	sum();
}

Xrf & Xrf::operator=(const Xrf & x)
{
	if (this != & x)
	{
		_y_mat = x.y_mat;
		_omega = x.omega;
		_ev_vec = x.ev_vec;
		_y_vec = x.y_vec;
		_Z_vec = x.Z_vec;
		_lines = x.lines;
		_row = x.row;
	}
	return *this;
}

void Xrf::add(const Xrf & x, bool mat_only)
{
	int Z, n_pixels;
	std::vector<int>::const_iterator Z_iter;
	if (!x.y_mat.empty())
		n_pixels = omega.theta.size()*omega.beta.size();
	for (int i = 0; i < x.Z_vec.size(); i++)
	{
		Z = x.Z_vec[i];
		Z_iter = std::find(Z_vec.begin(), Z_vec.end(), Z);
		if (Z_iter != Z_vec.end())
		{
			int p = Z_iter - Z_vec.begin();
			if (!mat_only)
			{
				std::vector<double>::const_iterator i2 = x.y_vec.begin()+x.row[i];
				std::vector<double>::iterator i1 = _y_vec.begin()+row[p];
				#pragma omp parallel for if(MP)
				for (int j = 0;	j < row[p+1]-row[p]; j++)
					*(i1+j) += *(i2+j);
			}
			if (n_pixels >= 1)
			{
				std::vector<double>::const_iterator i2 = x.y_mat.begin()+x.row[i]*n_pixels;
				std::vector<double>::iterator i1 = _y_mat.begin()+row[p]*n_pixels;
				#pragma omp parallel for if(MP)
				for (int j = 0;	j < (row[p+1] - row[p])*n_pixels; j++)
					*(i1+j) += *(i2+j);
			}
		}
		else
		{
			#pragma omp parallel sections if(MP)
			{
				{
					_Z_vec.push_back(Z);
					_lines.insert(_lines.end(),x.lines.begin()+x.row[i],x.lines.begin()+x.row[i+1]);
					_row.push_back(lines.size());
				}
				
				#pragma omp section
				{
					_ev_vec.insert(_ev_vec.end(),x.ev_vec.begin()+x.row[i],x.ev_vec.begin()+x.row[i+1]);
					_y_vec.insert(_y_vec.end(),x.y_vec.begin()+x.row[i],x.y_vec.begin()+x.row[i+1]);
				}
				
				#pragma omp section
				{
					if (n_pixels)
						_y_mat.insert(_y_mat.end(),x.y_mat.begin()+x.row[i]*n_pixels,x.y_mat.begin()+x.row[i+1]*n_pixels);
				}
			}
		}
	}
}

void Xrf::sum()
{
	if (!y_mat.empty())
	{
		_y_vec.resize(ev_vec.size(), 0.0);
		int n_pixels = omega.theta.size()*omega.beta.size();
		#pragma omp parallel for if(MP)
		for (int i = 0; i < ev_vec.size(); i++)
			_y_vec[i] = std::accumulate(y_mat.begin()+i*n_pixels,y_mat.begin()+(i+1)*n_pixels,0.0);
	}
}

void Xrf::show() const
{
	out(std::cout);
}

void Xrf::out(std::ostream & ost) const
{
	std::vector<int>::const_iterator Z, r, line;
	std::vector<double>::const_iterator ev, y;
	r = row.begin();
	y = y_vec.begin();
	ev = ev_vec.begin();
	line = lines.begin();
	
	ost << "# ======================================== #" << std::endl;
	ost << "# # # # # XRF lines: # # # # #" << std::endl;
	ost << "# ======================================== #" << std::endl;
	for (Z = Z_vec.begin(); Z < Z_vec.end(); Z++)
	{
		ost << "# ======================================== #" << std::endl;
		ost << "# Z = " << *Z <<  " (" << AtomicNumberToSymbol(*Z) << ")" << std::endl;
		ost << "# Line#\tE0(eV)\tIntensity\n";
		ost << "# ======================================== #" << std::endl;
		for (int i=(*r); i<(*(r+1)); i++)
		{
			ost << *(line++) << "\t";
			ost << *(ev++) << "\t";
			ost << *(y++) << std::endl;
		}
		ost << "# ======================================== #" << std::endl;
		r++;
	}
	ost << "# # # # # End of XRF # # # # #" << std::endl;
}

Rayleigh::Rayleigh()
	: y_mat(_y_mat), omega(_omega), y(_y), ev0(_ev0)
{}

Rayleigh::Rayleigh(double ev0_, const Compound & c, const solid_angle & omega_, double n_photons, bool calc_sum)
	: y_mat(_y_mat), omega(_omega), y(_y), ev0(_ev0)
{
	_ev0 = ev0_;
	_omega = omega_;
	_y_mat.resize(omega.theta.size()*omega.beta.size(),0.0);
	if (n_photons < 0)
		n_photons = 1;
	#pragma omp parallel for collapse(2) if(MP)
	for (int j = 0; j < omega.theta.size(); j++)
		for (int k = 0; k < omega.beta.size(); k++)
			_y_mat[j*omega.beta.size()+k] = c.dmac_rayleigh_pol(ev0, omega.theta[j], omega.beta[k])*omega.domega(omega.theta[j])*n_photons;
	if (calc_sum)
		sum();
}

Rayleigh::Rayleigh(const Sample & s, const Illumination & il, const solid_angle & omega_)
	: y_mat(_y_mat), omega(_omega), y(_y), ev0(_ev0)
{
	_ev0 = il.ev0;
	_omega = omega_;
	for (std::vector<Monolayer>::const_iterator layer = s.layer_vec.begin(); layer < s.layer_vec.end(); layer++)
	{
		std::cout << "Calculating Rayleigh scattering for layer " << (*layer).layer << std::endl;
		Rayleigh *temp_p = new Rayleigh(il.ev0, *layer, omega, il.n_photons, false), & temp = *temp_p;
		std::cout << "Calculating Rayleigh scattering attenuation for layer " << (*layer).layer << std::endl;
		#pragma omp parallel for collapse(2) if(MP)
		for (int theta = 0; theta < omega.theta.size(); theta++)
		{
			for (int beta = 0; beta < omega.beta.size(); beta++)
			{
				double psiprime = il.psi_prime(omega.theta[theta], omega.beta[beta]);
				int i = theta*omega.beta.size() + beta;
				//attenuation within the layer
				temp._y_mat[i] *= atten_mono(il.ev0, il.ev0, il.psi, psiprime, *layer);
				if (psiprime > 0) //reflection geometry
					//attenuation from layers upstream
					for (std::vector<Monolayer>::const_iterator layer_up = s.layer_vec.begin(); layer_up < layer; layer_up++)
						temp._y_mat[i] *= atten_refl(il.ev0, il.ev0, il.psi, psiprime, *layer_up);
				else if (psiprime < 0) //transmission geometry
				{
					//attenuation from layers upstream
					for (std::vector<Monolayer>::const_iterator layer_up = s.layer_vec.begin(); layer_up < layer; layer_up++)
						temp._y_mat[i] *= atten_trans_in(il.ev0, il.psi, *layer_up);
					//attenuation from layers downstream
					for (std::vector<Monolayer>::const_iterator layer_down = layer+1; layer_down < s.layer_vec.end(); layer_down++)
						temp._y_mat[i] *= atten_trans_out(il.ev0, psiprime, *layer_down);
				}
			}
		}
		add(temp, true);
		delete temp_p;
	}
	sum();
}

Rayleigh & Rayleigh::operator=(const Rayleigh & r)
{
	if (this != & r)
	{
		_y_mat = r.y_mat;
		_omega = r.omega;
		_y = r.y;
		_ev0 = r.ev0;
	}
	return *this;
}

void Rayleigh::add(const Rayleigh & x, bool mat_only)
{
	if (!x.y_mat.empty())
	{
		if (y_mat.empty())
			_y_mat = x.y_mat;
		else
			#pragma omp parallel for if(MP)
			for (int i = 0; i < y_mat.size(); i++)
				_y_mat[i] += x.y_mat[i];
	}
	if (!mat_only)
		_y += x.y;
}

void Rayleigh::sum()
{
	if (!y_mat.empty())
		_y = std::accumulate(y_mat.begin(), y_mat.end(), 0.0);
}

void Rayleigh::show() const
{
	out(std::cout);
}

void Rayleigh::out(std::ostream & ost) const
{
	ost << "# ======================================== #" << std::endl;
	ost << "# # # # # Rayleigh scattering # # # # #" << std::endl;
	ost << "# Peak @ " << ev0 << std::endl;
	ost << "# Intensity: " << y << std::endl;
	ost << "# # # # # End of Rayleigh scattering # # # # #" << std::endl;
	ost << "# ======================================== #" << std::endl;
}


Compton::Compton()
	: y_mat(_y_mat), omega(_omega), y_vec(_y_vec), ev_vec(_ev_vec), ev0(_ev0)
{}

Compton::Compton(double ev0_, const Compound & c, const solid_angle & omega_, double n_photons, bool calc_ev, bool calc_sum)
	: y_mat(_y_mat), omega(_omega), y_vec(_y_vec), ev_vec(_ev_vec), ev0(_ev0)
{
	_ev0 = ev0_;
	_omega = omega_;
	_y_mat.resize(omega.theta.size()*omega.beta.size(),0.0);
	double theta, beta;

	if (n_photons < 0)
		n_photons = 1;
	
	// #pragma omp parallel for collapse(2) private(theta, beta) if(MP)
	// for (int j = 0; j < omega.theta.size(); j++)
		// for (int k = 0; k < omega.beta.size(); k++)
		// {
			// theta = omega.theta[j];
			// if (theta < 1e-6)
				// continue;
			// beta = omega.beta[k];
			// _y_mat[j*omega.beta.size()+k] = c.dmac_compton_pol(ev0, theta, beta)*omega.domega(theta);
		// }
	
	for (int j = 0; j < omega.theta.size(); j++)
	{
		theta = omega.theta[j];
		if (theta < 1e-6)
			continue;
		#pragma omp parallel for if(MP)
		for (int k = 0; k < omega.beta.size(); k++)
		{
			beta = omega.beta[k];
			_y_mat[j*omega.beta.size()+k] = c.dmac_compton_pol(ev0, theta, beta)*omega.domega(theta)*n_photons;
		}
	}
	
	if (calc_ev)
	{
		_ev_vec.resize(omega.theta.size(),0.0);
		#pragma omp parallel for if(MP)
		for (int i = 0; i < omega.theta.size(); i++)
			_ev_vec[i] = ev_scattered(omega.theta[i], ev0);
	}
	if (calc_sum)
	{
		_y_vec.resize(omega.theta.size(),0.0);
		#pragma omp parallel for if(MP)
		for (int i = 0; i < omega.theta.size(); i++)
			_y_vec[i] = std::accumulate(y_mat.begin()+i*omega.beta.size(),y_mat.begin()+(i+1)*omega.beta.size(),0.0);
	}
}

Compton::Compton(const Sample & s, const Illumination & il, const solid_angle & omega_)
	: y_mat(_y_mat), omega(_omega), y_vec(_y_vec), ev_vec(_ev_vec), ev0(_ev0)
{
	_ev0 = il.ev0;
	_omega = omega_;
	ev();
	for (std::vector<Monolayer>::const_iterator layer = s.layer_vec.begin(); layer < s.layer_vec.end(); layer++)
	{
		std::cout << "Calculating Compton scattering for layer " << (*layer).layer << std::endl;
		Compton *temp_p = new Compton(il.ev0, *layer, omega, il.n_photons, false, false), & temp = *temp_p;
		std::cout << "Calculating Compton scattering attenuation for layer " << (*layer).layer << std::endl;
		#pragma omp parallel for collapse(2) if(MP)
		// for (std::vector<double>::const_iterator theta = omega.theta.begin(), ev = ev_vec.begin(); theta < omega.theta.end(); theta++, ev++)
		for (int theta = 0; theta < omega.theta.size(); theta++)
			for (int beta = 0; beta < omega.beta.size(); beta++)
			{
				double psiprime = il.psi_prime(omega.theta[theta], omega.beta[beta]);
				int i = theta*omega.beta.size() + beta;
				//attenuation within the layer
				temp._y_mat[i] *= atten_mono(il.ev0, ev_vec[theta], il.psi, psiprime, *layer);
				if (psiprime > 0) //reflection geometry
					//attenuation from layers upstream
					for (std::vector<Monolayer>::const_iterator layer_up = s.layer_vec.begin(); layer_up < layer; layer_up++)
						temp._y_mat[i] *= atten_refl(il.ev0, ev_vec[theta], il.psi, psiprime, *layer_up);
				else if (psiprime < 0) //transmission geometry
				{
					//attenuation from layers upstream
					for (std::vector<Monolayer>::const_iterator layer_up = s.layer_vec.begin(); layer_up < layer; layer_up++)
						temp._y_mat[i] *= atten_trans_in(il.ev0, il.psi, *layer_up);
					//attenuation from layers downstream
					for (std::vector<Monolayer>::const_iterator layer_down = layer+1; layer_down < s.layer_vec.end(); layer_down++)
						temp._y_mat[i] *= atten_trans_out(ev_vec[theta], psiprime, *layer_down);
				}
			}
		add(temp, true);
		delete temp_p;
	}
	sum();
}


Compton & Compton::operator=(const Compton & c)
{
	if (this != & c)
	{
		_y_mat = c.y_mat;
		_omega = c.omega;
		_y_vec = c.y_vec;
		_ev_vec = c.ev_vec;
		_ev0 = ev0;
	}
	return *this;
}

void Compton::add(const Compton & x, bool mat_only)
{
	if (!x.y_mat.empty())
	{
		if (y_mat.empty())
			_y_mat = x.y_mat;
		else
			#pragma omp parallel for if(MP)
			for (int i = 0; i < y_mat.size(); i++)
				_y_mat[i] += x.y_mat[i];
	}
	if (!mat_only)
		#pragma omp parallel for if(MP)
		for (int j = 0; j < x.y_vec.size(); j++)
			_y_vec[j] += x.y_vec[j];
}

void Compton::sum()
{
	if (!y_mat.empty())
	{
		_y_vec.resize(omega.theta.size());
		#pragma omp parallel for if(MP)
		for (int j = 0; j < omega.theta.size(); j++)
			_y_vec[j] = std::accumulate(y_mat.begin()+j*omega.beta.size(),y_mat.begin()+(j+1)*omega.beta.size(),0.0);
	}
}

void Compton::sum_ev()
{
	if (!y_mat.empty())
	{
		_y_vec.resize(omega.theta.size());
		_ev_vec.resize(omega.theta.size());
		double theta;
		#pragma omp parallel for private(theta) if(MP)
		for (int j = 0; j < omega.theta.size(); j++)
		{
			theta = omega.theta[j];
			_y_vec[j] = std::accumulate(y_mat.begin()+j*omega.beta.size(),y_mat.begin()+(j+1)*omega.beta.size(),0.0);
			_ev_vec[j] = ev_scattered(theta, ev0);
		}
	}
}

void Compton::ev()
{
	_ev_vec.resize(omega.theta.size());
	double theta;
	#pragma omp parallel for private(theta) if(MP)
	for (int j = 0; j < omega.theta.size(); j++)
	{
		theta = omega.theta[j];
		_ev_vec[j] = ev_scattered(theta, ev0);
	}
}

void Compton::show() const
{
	out(std::cout);
}

void Compton::out(std::ostream & ost) const
{
	ost << "# ======================================== #" << std::endl;
	ost << "# # # # # Compton scattering # # # # #" << std::endl;
	ost << "# Energy range: " << ev_vec.back() << " - " << ev_vec.front() << std::endl;
	ost << "# Total intensity: " << std::accumulate(y_vec.begin(), y_vec.end(), 0.0) << std::endl;
	ost << "# # # # # End of Compton scattering # # # # #" << std::endl;
	ost << "# ======================================== #" << std::endl;
}

Spectrum::Spectrum() 
	: y_vec(_y_vec), xrf(_xrf), ray(_ray), comp(_comp), omega(_omega), sample(_sample), illumination(_illumination), detector(_detector), Z_vec(xrf.Z_vec), lines(xrf.lines), row(xrf.row), y_sep(_y_sep) {}

Spectrum::Spectrum(const Sample & s, const Illumination & il, const solid_angle & omega_, const Detector & det, bool det_response, bool det_window, bool separate) 
	: y_vec(_y_vec), xrf(_xrf), ray(_ray), comp(_comp), omega(_omega), sample(_sample), illumination(_illumination), detector(_detector), Z_vec(xrf.Z_vec), lines(xrf.lines), row(xrf.row), y_sep(_y_sep)
{
	_omega = omega_;
	_illumination = il;
	_sample = s;
	_detector = det;
	
		_xrf = Xrf(sample, illumination, omega);
	#pragma omp parallel sections if(MP)
	{
		{_ray = Rayleigh(sample, illumination, omega);}
		#pragma omp section
		{_comp = Compton(sample, illumination, omega);}
	}
	
	_y_vec.clear();
	_y_vec.resize(detector.channel.n_channels);
	if (!separate)
	{
		detector.genspec(xrf.ev_vec, xrf.y_vec, _y_vec);
		detector.genspec(illumination.ev0, ray.y, _y_vec);
		detector.genspec(comp.ev_vec, comp.y_vec, _y_vec);
	}
	else
	{
		std::vector<int> row_temp{0, -1};
		detector.genspec(xrf.ev_vec, xrf.y_vec, _y_vec, _y_sep, xrf.row);
		detector.genspec(illumination.ev0, ray.y, _y_vec, _y_sep);
		detector.genspec(comp.ev_vec, comp.y_vec, _y_vec, _y_sep, row_temp);
	}
}

Spectrum::Spectrum(double ev0, const Compound & c, const solid_angle & omega_, const Detector & det, bool det_response, bool det_window, bool separate) 
	: y_vec(_y_vec), xrf(_xrf), ray(_ray), comp(_comp), omega(_omega), sample(_sample), illumination(_illumination), detector(_detector), Z_vec(xrf.Z_vec), lines(xrf.lines), row(xrf.row), y_sep(_y_sep)
{
	_omega = omega_;
	_detector = det;
	
	#pragma omp parallel sections if(MP)
	{
		{_xrf = Xrf(sample, illumination, omega);}
		#pragma omp section
		{_ray = Rayleigh(sample, illumination, omega);
		_comp = Compton(sample, illumination, omega);}
	}
	
	_y_vec.clear();
	_y_vec.resize(detector.channel.n_channels);
	std::cout << "Generating detector binned spectra..." << std::endl;
	if (!separate)
	{
		detector.genspec(xrf.ev_vec, xrf.y_vec, _y_vec);
		detector.genspec(illumination.ev0, ray.y, _y_vec);
		detector.genspec(comp.ev_vec, comp.y_vec, _y_vec);
	}
	else
	{
		std::vector<int> row_temp{0, -1};
		detector.genspec(xrf.ev_vec, xrf.y_vec, _y_vec, _y_sep, xrf.row);
		detector.genspec(illumination.ev0, ray.y, _y_vec, _y_sep);
		detector.genspec(comp.ev_vec, comp.y_vec, _y_vec, _y_sep, row_temp);
	}
}

void Spectrum::show() const
{
	out(std::cout);
}

void Spectrum::out(std::ostream & ost) const
{
	if (ost != std::cout)
	{
		ost << "# Total spectrum for each channel " << std::endl;
		for (auto i : y_vec)
			ost << i << "\t";
		ost << std::endl;
		ost << std::endl;
		if (!y_sep.empty())
		{
			ost << "# Separate spectra for each channel " << std::endl;
			int j = 0;
			for (auto i : y_sep)
			{
				ost << i << "\t";
				if (++j == y_vec.size())
				{
					ost << std::endl;
					ost << std::endl;
					j = 0;
				}
			}
			ost << std::endl;
		}
	}
	
	ost << "# ======================================== #" << std::endl;
	ost << "# # # # # Spectrum: # # # # #" << std::endl;
	ost << "# ======================================== #" << std::endl;
	xrf.out(ost);
	ray.out(ost);
	comp.out(ost);
	ost << "# ======================================== #" << std::endl;
	ost << "# # # # # End of spectrum # # # # #" << std::endl;
	ost << "# ======================================== #" << std::endl;
}
