Skip to content

Commit

Permalink
Merge pull request #152 from aglowacki/master
Browse files Browse the repository at this point in the history 
Added ability to disable Ka lines from fitting in override file
  • Loading branch information
Arthur Glowacki authored Dec 19, 2022
2 parents 17b0d3a + 60f96fc commit a23b977
Show file tree
Hide file tree
Showing 11 changed files with 173 additions and 78 deletions.
1 change: 0 additions & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -335,7 +335,6 @@ ELSE()
ENDIF()

#--------------- start pyxrf-maps lib -----------------

IF (BUILD_WITH_PYBIND11)
add_subdirectory(src/support/pybind11)
# .PYD file extension on Windows
Expand Down
1 change: 1 addition & 0 deletions reference/Scaler_to_PV_map.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ BeamLines:
OCR2: 2iddXMAP:dxp2:OutputCountRate
OCR3: 2iddXMAP:dxp3:OutputCountRate
OCR4: 2iddXMAP:dxp4:OutputCountRate
XBIC: 2idd:scaler1_cts1.D
TimeNormalizedScalers:
Timing: [2idd:3820:mca1.VAL, 50000000.0]
Scalers:
Expand Down
66 changes: 52 additions & 14 deletions src/data_struct/fit_element_map.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ POSSIBILITY OF SUCH DAMAGE.

#include "fit_element_map.h"
#include <iostream>
//#include "xraylib.h"

namespace data_struct
{
Expand Down Expand Up @@ -91,6 +92,11 @@ Fit_Element_Map<T_real>::Fit_Element_Map(std::string name, Element_Info<T_real>*
{
num_ratios = 12;
}
else if (_shell_type == "M")
{
num_ratios = 4;
}

for(size_t i=0; i<num_ratios; i++)
{
_energy_ratio_custom_multipliers.push_back(1.0);
Expand Down Expand Up @@ -127,19 +133,19 @@ void Fit_Element_Map<T_real>::init_energy_ratio_for_detector_element(const Eleme

if (_shell_type == "K") // K line
{
if(_pileup_element_info != nullptr)
if(_pileup_element_info != nullptr) // pileup's
{
_center = _element_info->xrf["ka1"] + _pileup_element_info->xrf["ka1"];
if (false == disable_Ka)
{
generate_energy_ratio(_element_info->xrf["ka1"] + _pileup_element_info->xrf["ka1"], 1.0, Element_Param_Type::Ka1_Line, detector_element);
generate_energy_ratio(_element_info->xrf["ka1"] + _pileup_element_info->xrf["kb1"], (_pileup_element_info->xrf_abs_yield["kb1"] / _pileup_element_info->xrf_abs_yield["ka1"]), Element_Param_Type::Kb1_Line, detector_element);
generate_energy_ratio(_element_info->xrf["ka1"] + _pileup_element_info->xrf["ka1"], _energy_ratio_custom_multipliers[0], Element_Param_Type::Ka1_Line, detector_element);
generate_energy_ratio(_element_info->xrf["ka1"] + _pileup_element_info->xrf["kb1"], (_pileup_element_info->xrf_abs_yield["kb1"] / _pileup_element_info->xrf_abs_yield["ka1"]) * _energy_ratio_custom_multipliers[1], Element_Param_Type::Kb1_Line, detector_element);
}
generate_energy_ratio(_element_info->xrf["kb1"] + _pileup_element_info->xrf["kb1"], (_element_info->xrf_abs_yield["kb1"] / _element_info->xrf_abs_yield["ka1"]) * (_pileup_element_info->xrf_abs_yield["kb1"] / _pileup_element_info->xrf_abs_yield["ka1"]), Element_Param_Type::Kb1_Line, detector_element);
generate_energy_ratio(_element_info->xrf["kb1"] + _pileup_element_info->xrf["kb1"], (_element_info->xrf_abs_yield["kb1"] / _element_info->xrf_abs_yield["ka1"]) * (_pileup_element_info->xrf_abs_yield["kb1"] / _pileup_element_info->xrf_abs_yield["ka1"]) * _energy_ratio_custom_multipliers[2], Element_Param_Type::Kb1_Line, detector_element);

if(_element_info != _pileup_element_info)
{
generate_energy_ratio(_element_info->xrf["kb1"] + _pileup_element_info->xrf["ka1"], (_element_info->xrf_abs_yield["kb1"] / _element_info->xrf_abs_yield["ka1"]), Element_Param_Type::Kb1_Line, detector_element);
generate_energy_ratio(_element_info->xrf["kb1"] + _pileup_element_info->xrf["ka1"], (_element_info->xrf_abs_yield["kb1"] / _element_info->xrf_abs_yield["ka1"]) * _energy_ratio_custom_multipliers[3], Element_Param_Type::Kb1_Line, detector_element);
}

}
Expand All @@ -148,7 +154,7 @@ void Fit_Element_Map<T_real>::init_energy_ratio_for_detector_element(const Eleme
_center = _element_info->xrf["ka1"];
if (false == disable_Ka)
{
generate_energy_ratio(_element_info->xrf["ka1"], 1.0, Element_Param_Type::Ka1_Line, detector_element);
generate_energy_ratio(_element_info->xrf["ka1"], _energy_ratio_custom_multipliers[0], Element_Param_Type::Ka1_Line, detector_element);
generate_energy_ratio(_element_info->xrf["ka2"], (_element_info->xrf_abs_yield["ka2"] / _element_info->xrf_abs_yield["ka1"]) * _energy_ratio_custom_multipliers[1], Element_Param_Type::Ka2_Line, detector_element);
}
generate_energy_ratio(_element_info->xrf["kb1"], (_element_info->xrf_abs_yield["kb1"] / _element_info->xrf_abs_yield["ka1"]) * _energy_ratio_custom_multipliers[2], Element_Param_Type::Kb1_Line, detector_element);
Expand Down Expand Up @@ -176,13 +182,45 @@ void Fit_Element_Map<T_real>::init_energy_ratio_for_detector_element(const Eleme
generate_energy_ratio(_element_info->xrf["ll"], (_element_info->xrf_abs_yield["ll"] / _element_info->xrf_abs_yield["la1"]) * _energy_ratio_custom_multipliers[10], Element_Param_Type::Ll_Line, detector_element);
generate_energy_ratio(_element_info->xrf["ln"], (_element_info->xrf_abs_yield["ln"] / _element_info->xrf_abs_yield["la1"]) * _energy_ratio_custom_multipliers[11], Element_Param_Type::Ln_Line, detector_element);
}
else if(_shell_type == "M")
else if (_shell_type == "M")
{
//TODO: finish adding M info
generate_energy_ratio(_element_info->xrf["ma1"], 1.0, Element_Param_Type::Ma1_Line, detector_element);
//generate_energy_ratio(_element_info->xrf["ma2"], 1.0, Element_Param_Type::Ma2_Line, detector_element);
//generate_energy_ratio(_element_info->xrf["mb"], 1.0, Element_Param_Type::Mb_Line, detector_element);
//generate_energy_ratio(_element_info->xrf["mg"], 1.0, Element_Param_Type::Mg_Line, detector_element);
// M5 - N7
T_real ratio = 1.0;
//ratio = RadRate(_element_info->number, M5N7_LINE, nullptr);
generate_energy_ratio(_element_info->xrf["ma1"], ratio, Element_Param_Type::Ma1_Line, detector_element);
// M5 - N6
if (_element_info->xrf_abs_yield["ma2"] == 0)
{
ratio = _energy_ratio_custom_multipliers[1];
}
else
{
ratio = (_element_info->xrf_abs_yield["ma2"] / _element_info->xrf_abs_yield["ma1"]) * _energy_ratio_custom_multipliers[1];
}
//ratio = RadRate(_element_info->number, M5N6_LINE, nullptr) * _energy_ratio_custom_multipliers[1]; //xraylib
generate_energy_ratio(_element_info->xrf["ma2"], ratio, Element_Param_Type::Ma2_Line, detector_element);
// M4 - N6
if (_element_info->xrf_abs_yield["mb"] == 0)
{
ratio = _energy_ratio_custom_multipliers[2];
}
else
{
ratio = (_element_info->xrf_abs_yield["mb"] / _element_info->xrf_abs_yield["ma1"]) * _energy_ratio_custom_multipliers[2];
}
//ratio = RadRate(_element_info->number, M4N6_LINE, nullptr) * _energy_ratio_custom_multipliers[2]; // xraylib
generate_energy_ratio(_element_info->xrf["mb"], ratio, Element_Param_Type::Mb_Line, detector_element);
// M3 - N5
if (_element_info->xrf_abs_yield["mg"] == 0)
{
ratio = _energy_ratio_custom_multipliers[3];
}
else
{
ratio = (_element_info->xrf_abs_yield["mg"] / _element_info->xrf_abs_yield["ma1"]) * _energy_ratio_custom_multipliers[3];
}
//ratio = RadRate(_element_info->number, M3N5_LINE, nullptr) * _energy_ratio_custom_multipliers[3]; //xraylib
generate_energy_ratio(_element_info->xrf["mg"], ratio, Element_Param_Type::Mg_Line, detector_element);
}

_width = int( std::sqrt( std::pow(ENERGY_RES_OFFSET, 2) + std::pow( (_center * ENERGY_RES_SQRT), 2) ) );
Expand Down Expand Up @@ -264,7 +302,7 @@ void Fit_Element_Map<T_real>::generate_energy_ratio(T_real energy, T_real ratio,
template<typename T_real>
void Fit_Element_Map<T_real>::set_custom_multiply_ratio(unsigned int idx, T_real multi)
{
if (idx > 0 && idx < _energy_ratio_custom_multipliers.size()) // index 0 has to be 1.0
if (idx < _energy_ratio_custom_multipliers.size())
{
_energy_ratio_custom_multipliers[idx] = multi;
}
Expand All @@ -275,7 +313,7 @@ void Fit_Element_Map<T_real>::set_custom_multiply_ratio(unsigned int idx, T_real
template<typename T_real>
void Fit_Element_Map<T_real>::multiply_custom_multiply_ratio(unsigned int idx, T_real multi)
{
if (idx > 0 && idx < _energy_ratio_custom_multipliers.size()) // index 0 has to be 1.0
if ( idx < _energy_ratio_custom_multipliers.size())
{
_energy_ratio_custom_multipliers[idx] *= multi;
}
Expand Down
60 changes: 26 additions & 34 deletions src/data_struct/params_override.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,28 @@ class DLL_EXPORT Params_Override
elements_to_fit.clear();
}

void _parse_ratios(vector<string> &branching_ratio_shell, int amt)
{
for (const auto& itr : branching_ratio_shell)
{
istringstream f(itr);
string s;
string el_name;
getline(f, s, ':');
getline(f, el_name, ',');
el_name.erase(std::remove_if(el_name.begin(), el_name.end(), ::isspace), el_name.end());


for (unsigned int i = 0; i < amt; i++)
{
T_real factor = 1.0;
std::getline(f, s, ',');
factor = parse_input_real<T_real>(s);
branching_ratios[el_name][i] = factor;
}
}
}

void parse_and_gen_branching_ratios()
{
branching_ratios.clear();
Expand Down Expand Up @@ -165,43 +187,12 @@ class DLL_EXPORT Params_Override
branching_ratios[el_name][10] = factor;

}
for (const auto& itr : branching_ratio_L)
{
istringstream f(itr);
string s;
string el_name;
getline(f, s, ':');
getline(f, el_name, ',');
el_name.erase(std::remove_if(el_name.begin(), el_name.end(), ::isspace), el_name.end());


for (unsigned int i = 0; i < 12; i++)
{
T_real factor = 1.0;
std::getline(f, s, ',');
factor = parse_input_real<T_real>(s);
branching_ratios[el_name][i] = factor;
}
}
for (const auto& itr : branching_ratio_K)
{
istringstream f(itr);
string s;
string el_name;
getline(f, s, ':');
getline(f, el_name, ',');
el_name.erase(std::remove_if(el_name.begin(), el_name.end(), ::isspace), el_name.end());

_parse_ratios(branching_ratio_K, 4);

for (unsigned int i = 0; i < 4; i++)
{
T_real factor = 1.0;
std::getline(f, s, ',');
factor = parse_input_real<T_real>(s);
branching_ratios[el_name][i] = factor;
}
}
_parse_ratios(branching_ratio_L, 12);

_parse_ratios(branching_ratio_M, 4);
}

map<int, T_real> get_custom_factor(string el_name)
Expand Down Expand Up @@ -237,6 +228,7 @@ class DLL_EXPORT Params_Override
vector<string> branching_family_L;
vector<string> branching_ratio_L;
vector<string> branching_ratio_K;
vector<string> branching_ratio_M;

unordered_map< string, map<int, T_real> > branching_ratios;

Expand Down
1 change: 1 addition & 0 deletions src/fitting/models/gaussian_model.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -677,6 +677,7 @@ const ArrayTr<T_real> Gaussian_Model<T_real>::compton_peak(const Fit_Parameters<
template<typename T_real>
const ArrayTr<T_real> Gaussian_Model<T_real>::escape_peak(const Fit_Parameters<T_real>* const fitp, const ArrayTr<T_real>& ev, T_real gain) const
{
// check detector element for how much to move escape peaks.
ArrayTr<T_real>counts(ev.size());
counts.setZero();
/*
Expand Down
33 changes: 32 additions & 1 deletion src/io/file/aps/aps_fit_params_import.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -368,7 +368,7 @@ DLL_EXPORT bool load_parameters_override(std::string path, Params_Override<T_rea
params_override->fit_params[tag_name].value = fvalue;
}
}
else if (tag == "BRANCHING_FAMILY_ADJUSTMENT_L" || tag == "BRANCHING_RATIO_ADJUSTMENT_L" || tag == "BRANCHING_RATIO_ADJUSTMENT_K")
else if (tag == "BRANCHING_FAMILY_ADJUSTMENT_L" || tag == "BRANCHING_RATIO_ADJUSTMENT_L" || tag == "BRANCHING_RATIO_ADJUSTMENT_K" || tag == "BRANCHING_RATIO_ADJUSTMENT_M")
{
unsigned int cnt = 0;

Expand All @@ -387,6 +387,11 @@ DLL_EXPORT bool load_parameters_override(std::string path, Params_Override<T_rea
params_override->branching_ratio_L.push_back(line);
cnt = 12;
}
else if (tag == "BRANCHING_RATIO_ADJUSTMENT_M")
{
params_override->branching_ratio_M.push_back(line);
cnt = 4;
}

std::string element_symb;
std::string str_value;
Expand Down Expand Up @@ -891,6 +896,10 @@ DLL_EXPORT bool save_parameters_override(std::string path, Params_Override<T_rea
{
out_stream << "BRANCHING_RATIO_ADJUSTMENT_L: " << itr.first;
}
else if (itr.second->shell_type_as_string() == "M")
{
out_stream << "BRANCHING_RATIO_ADJUSTMENT_M: " << itr.first;
}
for (float r : itr.second->energy_ratio_multipliers())
{
out_stream << "," << r;
Expand Down Expand Up @@ -962,6 +971,28 @@ DLL_EXPORT bool save_parameters_override(std::string path, Params_Override<T_rea
out_stream << itr << "\n";
}
}
out_stream << " this allows manual adjustment of the branhcing ratios between the different M lines, such as Ma1, Ma2, Mb, Mg\n";
out_stream << " Please note, these are all RELATIVE RELATIVE modifications, i.e., a 1 will correspond to exactly the literature value, etc.\n";
out_stream << " all will be normalized to the Ma1 line, and the values need to be in the following order:\n";
out_stream << " Ma1, Ma2, Mb, Mg\n";
out_stream << " please note, the first value (Ma1) MUST BE A 1. !!!\n";
for (const auto& itr : params_override->branching_ratio_M)
{
bool found = false;
for (std::string element_name : branching_ratios_updated)
{
if (itr.find(element_name) != std::string::npos)
{
found = true;
break;
}
}
if (false == found)
{
out_stream << itr << "\n";
}
}

out_stream << " the parameter adds the escape peaks (offset) to the fit if larger than 0. You should not enable Si and Ge at the same time, ie, one of these two values should be zero\n";
out_stream << "SI_ESCAPE_FACTOR: " << params_override->si_escape_factor << "\n";
out_stream << "GE_ESCAPE_FACTOR: " << params_override->ge_escape_factor << "\n";
Expand Down
9 changes: 6 additions & 3 deletions src/io/file/csv_io.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -215,10 +215,13 @@ namespace csv

for (const auto& s_itr : *standards)
{
if (s_itr.second.element_counts.at(routine).count(name) > 0)
if (s_itr.second.element_counts.count(routine) > 0)
{
counts = s_itr.second.element_counts.at(routine).at(name);
break;
if (s_itr.second.element_counts.at(routine).count(name) > 0)
{
counts = s_itr.second.element_counts.at(routine).at(name);
break;
}
}
}

Expand Down
6 changes: 5 additions & 1 deletion src/io/file/file_scan.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,7 @@ namespace io
_bnp_netcdf_files = find_all_dataset_files(dataset_dir + "flyXRF" + DIR_END_CHAR, "_001.nc");
_hdf_files = find_all_dataset_files(dataset_dir + "flyXRF.h5" + DIR_END_CHAR, "_0.h5");
_hdf_xspress_files = find_all_dataset_files(dataset_dir + "flyXspress" + DIR_END_CHAR, "_0.h5");
_hdf_xspress_files = find_all_dataset_files(dataset_dir + "flyXRF" + DIR_END_CHAR, "_0.hdf5");
//_hdf_confocal_files = find_all_dataset_files(dataset_dir , ".hdf5");
_hdf_emd_files = find_all_dataset_files(dataset_dir, ".emd");
}
Expand Down Expand Up @@ -277,7 +278,10 @@ namespace io
{
std::string full_path = dataset_directory + DIR_END_CHAR + "mda" + DIR_END_CHAR + itr;
long fsize = file::mda_get_multiplied_dims(full_path);
f_list.push_back(file_name_size(itr, fsize));
if (fsize > -1)
{
f_list.push_back(file_name_size(itr, fsize));
}
}
else
{
Expand Down
8 changes: 6 additions & 2 deletions src/io/file/hdf5_io.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1265,7 +1265,7 @@ class DLL_EXPORT HDF5_IO

buffer = new T_real[dims_in[0] * dims_in[2]]; // cols x spectra_size
count_row[0] = dims_in[0];
count_row[1] = 1;
count_row[1] = detector_num; // this should be detector number.
count_row[2] = dims_in[2];

size_t greater_cols = std::max(spec_row->size(), (size_t)dims_in[0]);
Expand Down Expand Up @@ -5531,7 +5531,11 @@ class DLL_EXPORT HDF5_IO
}

//save quantification_standard counts
std::unordered_map<std::string, T_real> element_counts = quant_itr.second.element_counts.at(fit_itr.first);
std::unordered_map<std::string, T_real> element_counts;
if (quant_itr.second.element_counts.count(fit_itr.first) > 0)
{
element_counts = quant_itr.second.element_counts.at(fit_itr.first);
}
count[0] = element_counts.size();

_create_memory_space(1, count, memoryspace_id);
Expand Down
22 changes: 21 additions & 1 deletion src/io/file/mca_io.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,6 +227,15 @@ DLL_EXPORT bool load_integrated_spectra(std::string path, data_struct::Spectra<T
if (f_idx > -1)
{
string pv_name = value.substr(0, f_idx);
string new_name;
string beamline;
bool is_time_norm;
// Translate pv_name with scaler ref file
if (data_struct::Scaler_Lookup::inst()->search_pv(pv_name, new_name, is_time_norm, beamline))
{
pv_name = new_name;
}

value = value.substr(f_idx + 1);
if (value[0] == '"')
{
Expand All @@ -237,7 +246,17 @@ DLL_EXPORT bool load_integrated_spectra(std::string path, data_struct::Spectra<T
try
{
pv_map[pv_name] = parse_input_real<T_real>(pv_value);

// update icr ocr in spectra
// if starts with ICR
if (pv_name.find("ICR") == 0)
{
spectra->input_counts(pv_map[pv_name]);
}
// if starts with OCR
if (pv_name.find("OCR") == 0)
{
spectra->output_counts(pv_map[pv_name]);
}
}
catch (std::exception& e)
{
Expand All @@ -259,6 +278,7 @@ DLL_EXPORT bool load_integrated_spectra(std::string path, data_struct::Spectra<T
float fvalue = parse_input_real<T_real>(value);
(*spectra)[i] = fvalue;
}
spectra->recalc_elapsed_livetime();
}
}
}
Expand Down
Loading

0 comments on commit a23b977

Please sign in to comment.