Skip to content

Commit

Permalink
Merge pull request #198 from aglowacki/master
Browse files Browse the repository at this point in the history 
Fixed typo in compton peak calc
  • Loading branch information
@aglowacki
aglowacki authored Aug 16, 2024
2 parents a66a5c1 + 4fc0789 commit dad709b
Show file tree
Hide file tree
Showing 4 changed files with 44 additions and 25 deletions.
31 changes: 21 additions & 10 deletions src/fitting/models/gaussian_model.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ Fit_Parameters<T_real> Gaussian_Model<T_real>::_generate_default_fit_parameters(
fit_params.add_parameter(Fit_Param<T_real>(STR_COMPTON_FWHM_CORR, (T_real)1.0, (T_real)4.0, (T_real)1.0, (T_real)0.1, E_Bound_Type::LIMITED_LO_HI));
fit_params.add_parameter(Fit_Param<T_real>(STR_COMPTON_AMPLITUDE, MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, (T_real)5.0, STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI));
fit_params.add_parameter(Fit_Param<T_real>(STR_COMPTON_F_STEP, (T_real)0.0, (T_real)1.0, (T_real)0.0, (T_real)0.1, E_Bound_Type::FIXED));
fit_params.add_parameter(Fit_Param<T_real>(STR_COMPTON_F_TAIL, (T_real)0.0, (T_real)1.0, (T_real)0.1, (T_real)0.1, E_Bound_Type::LIMITED_LO));
fit_params.add_parameter(Fit_Param<T_real>(STR_COMPTON_F_TAIL, (T_real)0.0, (T_real)1.0, (T_real)0.1, (T_real)0.1, E_Bound_Type::LIMITED_LO_HI));
fit_params.add_parameter(Fit_Param<T_real>(STR_COMPTON_GAMMA, (T_real)0.1, (T_real)10., (T_real)1.0, (T_real)0.1, E_Bound_Type::FIXED));
fit_params.add_parameter(Fit_Param<T_real>(STR_COMPTON_HI_F_TAIL, (T_real)0.0, (T_real)1.0, (T_real)0.013, (T_real)0.0000001, E_Bound_Type::LIMITED_LO_HI));
fit_params.add_parameter(Fit_Param<T_real>(STR_COMPTON_HI_GAMMA, (T_real)0.1, (T_real)3., (T_real)1.0, (T_real)0.01, E_Bound_Type::FIXED));
Expand Down Expand Up @@ -199,7 +199,7 @@ void Gaussian_Model<T_real>::set_fit_params_preset(Fit_Params_Preset preset)
_fit_parameters[STR_COMPTON_FWHM_CORR].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_AMPLITUDE].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_F_STEP].bound_type = E_Bound_Type::FIXED;
_fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO;
_fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_GAMMA].bound_type = E_Bound_Type::FIXED;
_fit_parameters[STR_COMPTON_HI_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_HI_GAMMA].bound_type = E_Bound_Type::FIXED;
Expand Down Expand Up @@ -237,12 +237,12 @@ void Gaussian_Model<T_real>::set_fit_params_preset(Fit_Params_Preset preset)
_fit_parameters[STR_COMPTON_FWHM_CORR].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_AMPLITUDE].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_F_STEP].bound_type = E_Bound_Type::FIXED;
_fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO;
_fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_GAMMA].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_HI_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_HI_GAMMA].bound_type = E_Bound_Type::LIMITED_LO_HI;

_fit_parameters[STR_SNIP_WIDTH].bound_type = E_Bound_Type::FIT;
_fit_parameters[STR_SNIP_WIDTH].bound_type = E_Bound_Type::LIMITED_LO_HI;

_fit_parameters[STR_F_STEP_OFFSET].bound_type = E_Bound_Type::FIXED;
_fit_parameters[STR_F_STEP_LINEAR].bound_type = E_Bound_Type::FIXED;
Expand Down Expand Up @@ -313,7 +313,7 @@ void Gaussian_Model<T_real>::set_fit_params_preset(Fit_Params_Preset preset)
_fit_parameters[STR_COMPTON_FWHM_CORR].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_AMPLITUDE].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_F_STEP].bound_type = E_Bound_Type::FIXED;
_fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO;
_fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_GAMMA].bound_type = E_Bound_Type::FIXED;
_fit_parameters[STR_COMPTON_HI_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI;
_fit_parameters[STR_COMPTON_HI_GAMMA].bound_type = E_Bound_Type::FIXED;
Expand Down Expand Up @@ -586,7 +586,7 @@ const Spectra<T_real> Gaussian_Model<T_real>::model_spectrum_element(const Fit_P
for (size_t idx = 0; idx < energy_ratios.size(); idx++)
{
const Element_Energy_Ratio<T_real>& er_struct = energy_ratios.at(idx);
T_real sigma = std::sqrt(std::pow((fitp->at(STR_FWHM_OFFSET).value / (T_real)2.3548), (T_real)2.0) + (er_struct.energy) * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value);
T_real sigma = std::sqrt( std::pow((fitp->at(STR_FWHM_OFFSET).value / (T_real)2.3548), (T_real)2.0) + (er_struct.energy * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value ) );
T_real f_step = std::abs<T_real>( er_struct.mu_fraction * ( fitp->at(STR_F_STEP_OFFSET).value + (fitp->at(STR_F_STEP_LINEAR).value * er_struct.energy)));
T_real f_tail = std::abs<T_real>( fitp->at(STR_F_TAIL_OFFSET).value + (fitp->at(STR_F_TAIL_LINEAR).value * er_struct.mu_fraction));
T_real kb_f_tail = std::abs<T_real>( fitp->at(STR_KB_F_TAIL_OFFSET).value + (fitp->at(STR_KB_F_TAIL_LINEAR).value * er_struct.mu_fraction));
Expand Down Expand Up @@ -725,15 +725,26 @@ const ArrayTr<T_real> Gaussian_Model<T_real>::peak(T_real gain, T_real sigma, co
template<typename T_real>
const ArrayTr<T_real> Gaussian_Model<T_real>::step(T_real gain, T_real sigma, const ArrayTr<T_real>& delta_energy, T_real peak_E) const
{
ArrayTr<T_real>counts = delta_energy.unaryExpr([gain, sigma, peak_E](T_real v) { return gain / (T_real)2.0 / peak_E * std::erfc(v / ((T_real)(M_SQRT2)*sigma)); });
/*
ArrayTr<T_real>counts = delta_energy.unaryExpr([gain, sigma, peak_E](T_real v)
{
T_real val = ((T_real)(M_SQRT2) * sigma);
val = std::erfc(v / val);
val = peak_E * val;
val = (T_real)2.0 / val;
return gain / val;
});
return counts;
*/
return delta_energy.unaryExpr([gain, sigma, peak_E](T_real v) { return gain / (T_real)2.0 / peak_E * std::erfc(v / ((T_real)(M_SQRT2)*sigma)); });

/*
ArrayTr<T_real>counts(delta_energy.size());
for (unsigned int i = 0; i < delta_energy.size(); i++)
{
counts[i] = gain / (T_real)2.0 / peak_E * std::erfc((T_real)delta_energy[i] / ((T_real)(M_SQRT2) * sigma));
}
*/
return counts;

}

Expand Down Expand Up @@ -764,7 +775,7 @@ const ArrayTr<T_real> Gaussian_Model<T_real>::elastic_peak(const Fit_Parameters<
{
ArrayTr<T_real> counts(ev.size());
counts.setZero();
T_real sigma = std::sqrt( std::pow( (fitp->at(STR_FWHM_OFFSET).value / (T_real)2.3548), (T_real)2.0 ) + fitp->at(STR_COHERENT_SCT_ENERGY).value * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value );
T_real sigma = std::sqrt( std::pow( (fitp->at(STR_FWHM_OFFSET).value / (T_real)2.3548), (T_real)2.0 ) + (fitp->at(STR_COHERENT_SCT_ENERGY).value * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value ) );
if(false == std::isfinite(sigma))
{
logE << "sigma = " << sigma << "\n";
Expand Down Expand Up @@ -798,7 +809,7 @@ const ArrayTr<T_real> Gaussian_Model<T_real>::compton_peak(const Fit_Parameters<

T_real compton_E = fitp->at(STR_COHERENT_SCT_ENERGY).value/((T_real)1.0 +(fitp->at(STR_COHERENT_SCT_ENERGY).value / (T_real)511.0 ) * ((T_real)1.0 -std::cos( fitp->at(STR_COMPTON_ANGLE).value * (T_real)2.0 * (T_real)(M_PI) / (T_real)360.0 )));

T_real sigma = std::sqrt( std::pow( (fitp->at(STR_FWHM_OFFSET).value/(T_real)2.3548), (T_real)62.0) + compton_E * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value );
T_real sigma = std::sqrt( std::pow( (fitp->at(STR_FWHM_OFFSET).value/(T_real)2.3548), (T_real)2.0) + (compton_E * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value ) );
if(false == std::isfinite(sigma))
{
logE << "sigma = " << sigma << "\n";
Expand Down
32 changes: 20 additions & 12 deletions src/fitting/optimizers/mpfit_optimizer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,7 +321,7 @@ void MPFit_Optimizer<T_real>::set_options(std::unordered_map<std::string, T_real


template<typename T_real>
void MPFit_Optimizer<T_real>::_fill_limits(Fit_Parameters<T_real> *fit_params , std::vector<struct mp_par<T_real> > &par)
bool MPFit_Optimizer<T_real>::_fill_limits(Fit_Parameters<T_real> *fit_params , std::vector<struct mp_par<T_real> > &par)
{
for (auto itr = fit_params->begin(); itr != fit_params->end(); itr++)
{
Expand All @@ -346,13 +346,11 @@ void MPFit_Optimizer<T_real>::_fill_limits(Fit_Parameters<T_real> *fit_params ,
|| fit.bound_type == E_Bound_Type::LIMITED_LO
|| fit.bound_type == E_Bound_Type::LIMITED_LO_HI)
{
if (fit.max_val == fit.min_val)
{
fit.max_val += (T_real)1.0;
fit.min_val -= (T_real)1.0;
(*fit_params)[itr->first].max_val += (T_real)1.0;
(*fit_params)[itr->first].min_val -= (T_real)1.0;
}
if(fit.min_val == fit.max_val)
{
logE<<"Fit parameter "<<fit.name<<" is limited and has min: "<<fit.min_val<<" == to max: "<<fit.max_val<<"\n";
return false;
}
}


Expand Down Expand Up @@ -410,6 +408,7 @@ void MPFit_Optimizer<T_real>::_fill_limits(Fit_Parameters<T_real> *fit_params ,
par[fit.opt_array_index].deriv_abstol = (T_real)0.00001; // Absolute tolerance for derivative debug printout
}
}
return true;
}

//-----------------------------------------------------------------------------
Expand Down Expand Up @@ -483,7 +482,7 @@ OPTIMIZER_OUTCOME MPFit_Optimizer<T_real>::minimize(Fit_Parameters<T_real>*fit_p
std::vector<T_real> resid(energy_range.count());
std::vector<T_real> covar(fitp_arr.size() * fitp_arr.size());

size_t total_itr = num_itr * (fitp_arr.size() + 1);
size_t total_itr = num_itr;
fill_user_data(ud, fit_params, spectra, elements_to_fit, model, energy_range, status_callback, total_itr, use_weights);

/*
Expand Down Expand Up @@ -521,7 +520,10 @@ OPTIMIZER_OUTCOME MPFit_Optimizer<T_real>::minimize(Fit_Parameters<T_real>*fit_p

_options.maxfev = _options.maxiter * ((int)fitp_arr.size() + 1);

_fill_limits(fit_params, par);
if(false == _fill_limits(fit_params, par))
{
return OPTIMIZER_OUTCOME::FAILED;
}

mp_result<T_real> result;
result.xerror = &perror[0];
Expand Down Expand Up @@ -669,7 +671,10 @@ OPTIMIZER_OUTCOME MPFit_Optimizer<T_real>::minimize_func(Fit_Parameters<T_real>

std::vector<struct mp_par<T_real> > par;
par.resize(fitp_arr.size());
_fill_limits(fit_params, par);
if(false == _fill_limits(fit_params, par))
{
return OPTIMIZER_OUTCOME::FAILED;
}

mp_result<T_real> result;
result.xerror = &perror[0];
Expand Down Expand Up @@ -821,7 +826,10 @@ OPTIMIZER_OUTCOME MPFit_Optimizer<T_real>::minimize_quantification(Fit_Parameter

std::vector<struct mp_par<T_real> > par;
par.resize(fitp_arr.size());
_fill_limits(fit_params, par);
if(false == _fill_limits(fit_params, par))
{
return OPTIMIZER_OUTCOME::FAILED;
}

this->_last_outcome = mpfit(quantification_residuals_mpfit<T_real>, (int)ud.quant_map.size(), (int)fitp_arr.size(), &fitp_arr[0], &par[0], &_options, (void *) &ud, &result);
logI << "\nOutcome: " << optimizer_outcome_to_str(this->_outcome_map[this->_last_outcome]) << "\nNum iter: " << result.niter << "\n Norm of the residue vector: " << *result.resid << "\n";
Expand Down
4 changes: 2 additions & 2 deletions src/fitting/optimizers/mpfit_optimizer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,9 +95,9 @@ class DLL_EXPORT MPFit_Optimizer: public Optimizer<T_real>

private:

void _fill_limits(Fit_Parameters<T_real> *fit_params, std::vector<struct mp_par<T_real> > &par);
bool _fill_limits(Fit_Parameters<T_real> *fit_params, std::vector<struct mp_par<T_real> > &par);

struct mp_config<T_real> _options;
struct mp_config<T_real> _options;

};

Expand Down
2 changes: 1 addition & 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 @@ -1076,7 +1076,7 @@ DLL_EXPORT bool save_parameters_override(std::string path, Params_Override<T_rea
//-----------------------------------------------------------------------------

template<typename T_real>
DLL_EXPORT bool save_fit_parameters_override(std::string path, Fit_Parameters<T_real> fit_params, std::string result)
DLL_EXPORT bool save_fit_parameters_override(const std::string& path, const Fit_Parameters<T_real>& fit_params, const std::string& result)
{

std::ofstream out_stream(path);
Expand Down

0 comments on commit dad709b

Please sign in to comment.