diff --git a/src/core/defines.h b/src/core/defines.h index 36b9f98f..35b98066 100644 --- a/src/core/defines.h +++ b/src/core/defines.h @@ -111,16 +111,19 @@ T_real parse_input_real(std::string value) } } -#define CALIBRATION_CURVE_SIZE 92 +constexpr auto CALIBRATION_CURVE_SIZE = 92; //namespace keys //{ -#define AVOGADRO 6.02204531e23 -#define HC_ANGSTROMS 12398.52 -#define RE 2.817938070e-13 // in cm -#define ENERGY_RES_OFFSET 150.0 -#define ENERGY_RES_SQRT 12.0 +constexpr auto AVOGADRO = 6.02204531e23; +constexpr auto HC_ANGSTROMS = 12398.52; +constexpr auto RE = 2.817938070e-13; // in cm +constexpr auto ENERGY_RES_OFFSET = 150.0; +constexpr auto ENERGY_RES_SQRT = 12.0; +constexpr auto MIN_COUNTS_LIMIT_LOG = -11.0; +constexpr auto MAX_COUNTS_LIMIT_LOG = 20.0; +constexpr auto STEP_COUNTS_LIMIT_LOG = 0.0005; /** * @brief String defines for fit parameters string value pair. diff --git a/src/core/process_whole.cpp b/src/core/process_whole.cpp index 67aeeab6..0a52a15d 100644 --- a/src/core/process_whole.cpp +++ b/src/core/process_whole.cpp @@ -386,7 +386,7 @@ void load_and_fit_quatification_datasets(data_struct::Analysis_Job* anal { if (false == override_params->fit_params.contains(itr2.first)) { - fit_params.add_parameter(Fit_Param(itr2.first, 1.0e-11, 20.0, quantification_standard->element_counts.at(fit_itr.first).at(itr2.first), 0.0005, E_Bound_Type::LIMITED_LO_HI)); + fit_params.add_parameter(Fit_Param(itr2.first, MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, quantification_standard->element_counts.at(fit_itr.first).at(itr2.first), STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI)); } } fitting::routines::Matrix_Optimized_Fit_Routine* f_routine = (fitting::routines::Matrix_Optimized_Fit_Routine*)fit_routine; diff --git a/src/data_struct/fit_parameters.cpp b/src/data_struct/fit_parameters.cpp index 1848cca1..7e200108 100644 --- a/src/data_struct/fit_parameters.cpp +++ b/src/data_struct/fit_parameters.cpp @@ -260,18 +260,34 @@ void Fit_Parameters::update_values(const Fit_Parameters *overri { itr.second.value = override_fit_params->at(itr.first).value; } + else + { + logW << itr.first << " value = " << itr.second.value << " . Not updating value\n"; + } if( std::isfinite(override_fit_params->at(itr.first).min_val) ) { itr.second.min_val = override_fit_params->at(itr.first).min_val; } + else + { + logW << itr.first << " min_val = " << itr.second.min_val << " . Not updating min_val\n"; + } if( std::isfinite(override_fit_params->at(itr.first).max_val) ) { itr.second.max_val = override_fit_params->at(itr.first).max_val; } + else + { + logW << itr.first << " max_val = " << itr.second.max_val << " . Not updating max_val\n"; + } if( override_fit_params->at(itr.first).bound_type != E_Bound_Type::NOT_INIT) { itr.second.bound_type = override_fit_params->at(itr.first).bound_type; } + else + { + logW << itr.first << " bound_type == E_Bound_Type::NOT_INIT. Not updating bound_type\n"; + } } } } diff --git a/src/data_struct/spectra.h b/src/data_struct/spectra.h index 1ed7e7e8..f8c030ac 100644 --- a/src/data_struct/spectra.h +++ b/src/data_struct/spectra.h @@ -172,21 +172,37 @@ class Spectra : public ArrayTr<_T> { _elapsed_livetime += val; } + else + { + logW << " Spectra _elapsed_livetime = " << val << " . Not updating!\n"; + } val = spectra.elapsed_realtime(); if(std::isfinite(val)) { _elapsed_realtime += val; } + else + { + logW << " Spectra _elapsed_realtime = " << val << " . Not updating!\n"; + } val = spectra.input_counts(); if(std::isfinite(val)) { _input_counts += val; } + else + { + logW << " Spectra _input_counts = " << val << " . Not updating!\n"; + } val = spectra.output_counts(); if(std::isfinite(val)) { _output_counts += val; } + else + { + logW << " Spectra _output_counts = " << val << " . Not updating!\n"; + } } // Note: T_real is different template typename, this is to convert double to float or the other way around. @@ -204,21 +220,37 @@ class Spectra : public ArrayTr<_T> { _elapsed_livetime += val; } + else + { + logW << " Spectra _elapsed_livetime = " << val << " . Not updating!\n"; + } val = spectra->elapsed_realtime(); if (std::isfinite(val)) { _elapsed_realtime += val; } + else + { + logW << " Spectra _elapsed_realtime = " << val << " . Not updating!\n"; + } val = spectra->input_counts(); if (std::isfinite(val)) { _input_counts += val; } + else + { + logW << " Spectra _input_counts = " << val << " . Not updating!\n"; + } val = spectra->output_counts(); if (std::isfinite(val)) { _output_counts += val; } + else + { + logW << " Spectra _output_counts = " << val << " . Not updating!\n"; + } } } diff --git a/src/fitting/models/gaussian_model.cpp b/src/fitting/models/gaussian_model.cpp index 9a39e69e..9600c582 100644 --- a/src/fitting/models/gaussian_model.cpp +++ b/src/fitting/models/gaussian_model.cpp @@ -104,11 +104,11 @@ Fit_Parameters Gaussian_Model::_generate_default_fit_parameters( fit_params.add_parameter(Fit_Param(STR_FWHM_FANOPRIME, (T_real)0.000001, (T_real)0.05, (T_real)0.00012, (T_real)0.000001, E_Bound_Type::LIMITED_LO_HI)); fit_params.add_parameter(Fit_Param(STR_COHERENT_SCT_ENERGY, (T_real)9.4, (T_real)10.4, (T_real)9.99, (T_real)0.001, E_Bound_Type::LIMITED_LO_HI)); - fit_params.add_parameter(Fit_Param(STR_COHERENT_SCT_AMPLITUDE, (T_real)0.000001, (T_real)10.0, (T_real)5.0, (T_real)0.0005, E_Bound_Type::LIMITED_LO_HI)); + fit_params.add_parameter(Fit_Param(STR_COHERENT_SCT_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(STR_COMPTON_ANGLE, (T_real)0.0, (T_real)359.0, (T_real)90.0, (T_real)0.1, E_Bound_Type::LIMITED_LO_HI)); fit_params.add_parameter(Fit_Param(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(STR_COMPTON_AMPLITUDE, (T_real)0.00001, (T_real)10.00, (T_real)5.0, (T_real)0.0005, E_Bound_Type::LIMITED_LO_HI)); + fit_params.add_parameter(Fit_Param(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(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(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(STR_COMPTON_GAMMA, (T_real)0.1, (T_real)10., (T_real)1.0, (T_real)0.1, E_Bound_Type::FIXED)); @@ -575,6 +575,7 @@ const Spectra Gaussian_Model::model_spectrum_element(const Fit_P if (false == std::isfinite(pre_faktor)) { + logE << "Prefactor = " << pre_faktor << "\n"; spectra_model = (ArrayTr)(spectra_model).unaryExpr([](T_real v) { return std::numeric_limits::quiet_NaN(); }); return spectra_model; } @@ -725,13 +726,14 @@ const ArrayTr Gaussian_Model::peak(T_real gain, T_real sigma, co template const ArrayTr Gaussian_Model::step(T_real gain, T_real sigma, const ArrayTr& delta_energy, T_real peak_E) const { - //delta_energy = delta_energy.unaryExpr([sigma](T_real v) { return (T_real)std::erfc(v /( M_SQRT2*sigma)); } ); - //return (gain / (T_real)2.0 / peak_E * delta_energy ); + ArrayTrcounts = 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)); }); + /* ArrayTrcounts(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; } @@ -739,10 +741,20 @@ const ArrayTr Gaussian_Model::step(T_real gain, T_real sigma, co // ---------------------------------------------------------------------------- template -const ArrayTr Gaussian_Model::tail(T_real gain, T_real sigma, ArrayTr delta_energy, T_real gamma) const +const ArrayTr Gaussian_Model::tail(T_real gain, T_real sigma, const ArrayTr &delta_energy, T_real gamma) const { - delta_energy = delta_energy.unaryExpr([sigma,gamma](T_real v) { return (v < (T_real)0.0) ? std::exp(v/ (gamma * sigma)) * std::erfc(v / ((T_real)(M_SQRT2)*sigma) + ((T_real)1.0/(gamma*(T_real)(M_SQRT2)))) : std::erfc(v / ((T_real)(M_SQRT2)*sigma) + ((T_real)1.0/(gamma*(T_real)(M_SQRT2)))); } ); - return( gain / (T_real)2.0 / gamma / sigma / exp((T_real)-0.5/pow(gamma, (T_real)2.0)) * delta_energy); + T_real val = pow(gamma, (T_real)2.0); + val = exp((T_real)-0.5 / val); + val = sigma / val; + val = gamma / val; + val = (T_real)2.0 / val; + val = gain / val; + + return delta_energy.unaryExpr([val,sigma,gamma](T_real v) { return (v < (T_real)0.0) ? + std::exp(v/ (gamma * sigma)) * val * std::erfc(v / ((T_real)(M_SQRT2)*sigma) + ((T_real)1.0/(gamma*(T_real)(M_SQRT2)))) + : + val * std::erfc(v / ((T_real)(M_SQRT2)*sigma) + ((T_real)1.0/(gamma*(T_real)(M_SQRT2)))); } ); + } // ---------------------------------------------------------------------------- @@ -755,6 +767,7 @@ const ArrayTr Gaussian_Model::elastic_peak(const Fit_Parameters< 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"; counts = (ArrayTr)(counts).unaryExpr([](T_real v) { return std::numeric_limits::quiet_NaN(); }); return counts; } @@ -787,6 +800,7 @@ const ArrayTr Gaussian_Model::compton_peak(const Fit_Parameters< 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 ); if(false == std::isfinite(sigma)) { + logE << "sigma = " << sigma << "\n"; counts = (ArrayTr)(counts).unaryExpr([](T_real v) { return std::numeric_limits::quiet_NaN(); }); return counts; } diff --git a/src/fitting/models/gaussian_model.h b/src/fitting/models/gaussian_model.h index 7760452d..b7e1de19 100644 --- a/src/fitting/models/gaussian_model.h +++ b/src/fitting/models/gaussian_model.h @@ -114,7 +114,7 @@ class DLL_EXPORT Gaussian_Model: public Base_Model */ const ArrayTr step(T_real gain, T_real sigma, const ArrayTr& delta_energy, T_real peak_E) const; - const ArrayTr tail(T_real gain, T_real sigma, ArrayTr delta_energy, T_real gamma) const; + const ArrayTr tail(T_real gain, T_real sigma, const ArrayTr& delta_energy, T_real gamma) const; virtual const ArrayTr elastic_peak(const Fit_Parameters* const fitp, const ArrayTr& ev, T_real gain) const; diff --git a/src/fitting/optimizers/lmfit_optimizer.cpp b/src/fitting/optimizers/lmfit_optimizer.cpp index ea08568e..9d5f1292 100644 --- a/src/fitting/optimizers/lmfit_optimizer.cpp +++ b/src/fitting/optimizers/lmfit_optimizer.cpp @@ -185,7 +185,7 @@ void quantification_residuals_lmfit( const T_real *par, int m_dat, const void *d { if (std::isfinite(result_map[itr.first]) == false) { - fvec[idx] = itr.second.e_cal_ratio; + fvec[idx] = itr.second.e_cal_ratio * 10.0; } else { diff --git a/src/fitting/optimizers/mpfit_optimizer.cpp b/src/fitting/optimizers/mpfit_optimizer.cpp index 863d17db..0dcb9425 100644 --- a/src/fitting/optimizers/mpfit_optimizer.cpp +++ b/src/fitting/optimizers/mpfit_optimizer.cpp @@ -84,22 +84,16 @@ int residuals_mpfit(int m, int params_size, T_real *params, T_real *dy, T_real * // Add background ud->spectra_model += ud->spectra_background; // Remove nan's and inf's - ud->spectra_model = (ArrayTr)ud->spectra_model.unaryExpr([ud](T_real v) { return std::isfinite(v) ? v : ud->normalizer; }); + //ud->spectra_model = (ArrayTr)ud->spectra_model.unaryExpr([ud](T_real v) { return std::isfinite(v) ? v : ud->normalizer; }); + T_real sum = 0.0; //Calculate residuals for (int i=0; ispectra[i] / ud->norm_arr[i]; - //T_real n_model = ud->spectra_model[i] / ud->norm_arr[i]; - //dy[i] = pow((n_raw - n_model), (T_real)2.0) * ud->weights[i]; - - T_real n_raw = ud->spectra[i] / ud->normalizer; - T_real n_model = ud->spectra_model[i] / ud->normalizer; - dy[i] = pow((n_raw - n_model), (T_real)2.0) * ud->weights[i]; - + dy[i] = pow((ud->spectra[i] - ud->spectra_model[i]), (T_real)2.0) * ud->weights[i]; + sum += dy[i]; - //dy[i] = pow((ud->spectra[i] - ud->spectra_model[i]), (T_real)2.0) * ud->weights[i]; - if (std::isfinite(dy[i]) == false) + if (std::isfinite(dy[i]) == false) { if(first) { @@ -119,7 +113,7 @@ int residuals_mpfit(int m, int params_size, T_real *params, T_real *dy, T_real * dy[i] = ud->normalizer; } } - + logI << "f = " << sum << "\n"; ud->cur_itr++; if (ud->status_callback != nullptr) { @@ -202,7 +196,7 @@ int quantification_residuals_mpfit(int m, int params_size, T_real *params, T_rea { if (std::isfinite(result_map[itr.first]) == false) { - dy[idx] = itr.second.e_cal_ratio; + dy[idx] = itr.second.e_cal_ratio * 10.0; } else { @@ -391,7 +385,7 @@ void MPFit_Optimizer::_fill_limits(Fit_Parameters *fit_params , break; } - par[fit.opt_array_index].step = fit.step_size; // 0 = auto ,Step size for finite difference + par[fit.opt_array_index].step = 0; // 0 = auto ,Step size for finite difference par[fit.opt_array_index].parname = 0; par[fit.opt_array_index].relstep = 0; // Relative step size for finite difference par[fit.opt_array_index].side = 0; // Sidedness of finite difference derivative diff --git a/src/fitting/routines/hybrid_param_nnls_fit_routine.cpp b/src/fitting/routines/hybrid_param_nnls_fit_routine.cpp index c8b31507..60210647 100644 --- a/src/fitting/routines/hybrid_param_nnls_fit_routine.cpp +++ b/src/fitting/routines/hybrid_param_nnls_fit_routine.cpp @@ -179,7 +179,7 @@ OPTIMIZER_OUTCOME Hybrid_Param_NNLS_Fit_Routine::fit_spectra_parameters( out_fit_params.append_and_update(fit_params); for (const auto& itr : *elements_to_fit) { - out_fit_params.add_parameter(Fit_Param(itr.first, 1.0e-11, 20.0, log10(out_counts[itr.first]), 0.0005, E_Bound_Type::LIMITED_LO_HI)); + out_fit_params.add_parameter(Fit_Param(itr.first, MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, log10(out_counts[itr.first]), STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI)); } } } diff --git a/src/fitting/routines/matrix_optimized_fit_routine.cpp b/src/fitting/routines/matrix_optimized_fit_routine.cpp index 16621645..6d0daed5 100644 --- a/src/fitting/routines/matrix_optimized_fit_routine.cpp +++ b/src/fitting/routines/matrix_optimized_fit_routine.cpp @@ -125,7 +125,7 @@ std::unordered_map> Matrix_Optimized_Fit_Routine fp(itr.first, (T_real)-10.0, (T_real)20.0, (T_real)0.0, (T_real)0.0005, E_Bound_Type::LIMITED_LO_HI); + Fit_Param fp(itr.first, MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, (T_real)0.0, STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI); fit_parameters[itr.first] = fp; } else @@ -189,7 +189,7 @@ std::unordered_map> Matrix_Optimized_Fit_Routine fp(itr.first, (T_real)-10.0, (T_real)20.0, (T_real)0.0, (T_real)0.0005, E_Bound_Type::LIMITED_LO_HI); + Fit_Param fp(itr.first, MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, (T_real)0.0, STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI); fit_parameters[itr.first] = fp; } else @@ -222,7 +222,7 @@ std::unordered_map> Matrix_Optimized_Fit_Routine fp(itr.first, (T_real)-10.0, (T_real)20.0, (T_real)0.0, (T_real)0.0005, E_Bound_Type::LIMITED_LO_HI); + Fit_Param fp(itr.first, MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, (T_real)0.0, STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI); fit_parameters[itr.first] = fp; } else diff --git a/src/fitting/routines/param_optimized_fit_routine.cpp b/src/fitting/routines/param_optimized_fit_routine.cpp index c75ecc08..5b7099cf 100644 --- a/src/fitting/routines/param_optimized_fit_routine.cpp +++ b/src/fitting/routines/param_optimized_fit_routine.cpp @@ -106,7 +106,7 @@ void Param_Optimized_Fit_Routine::_add_elements_to_fit_parameters(Fit_Pa //if element counts is not in fit params structure, add it if( false == fit_params->contains(el_itr.first) ) { - Fit_Param fp(element->full_name(), (T_real)-11.0, 20, e_guess, (T_real)0.0005, E_Bound_Type::LIMITED_LO_HI); + Fit_Param fp(element->full_name(), MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, e_guess, STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI); (*fit_params)[el_itr.first] = fp; } if(spectra != nullptr && energies.size() > 0) diff --git a/src/io/file/mca_io.h b/src/io/file/mca_io.h index eb8c6e16..8d94e7f6 100644 --- a/src/io/file/mca_io.h +++ b/src/io/file/mca_io.h @@ -104,6 +104,7 @@ DLL_EXPORT bool load_integrated_spectra(std::string path, data_struct::Spectra* fit_route, float val = out_counts.at(itr.first); if (val != 0.0) { - fit_params[itr.first] = data_struct::Fit_Param(itr.first, 1.0e-11, 20.0 log10(val), 0.0005, E_Bound_Type::LIMITED_LO_HI); + fit_params[itr.first] = data_struct::Fit_Param(itr.first, MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, log10(val), STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI); } else { - fit_params[itr.first] = data_struct::Fit_Param(itr.first, 1.0e-11, 20.0, val, 0.0005, E_Bound_Type::LIMITED_LO_HI); + fit_params[itr.first] = data_struct::Fit_Param(itr.first, MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, val, STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI); } } diff --git a/src/quantification/models/quantification_model.cpp b/src/quantification/models/quantification_model.cpp index fd29259a..fa2c5992 100644 --- a/src/quantification/models/quantification_model.cpp +++ b/src/quantification/models/quantification_model.cpp @@ -259,6 +259,7 @@ std::unordered_map Quantification_Model::model_cali T_real val = p * itr.second.absorption * itr.second.transmission_Be * itr.second.transmission_Ge * itr.second.yield * ((T_real)1. - itr.second.transmission_through_Si_detector) * itr.second.transmission_through_air; if(false == std::isfinite(val)) { + logW << "val = " << val << ". Defaulting to 0.0\n"; val = 0; } result_map.emplace(std::pair(itr.first, val)); @@ -278,6 +279,7 @@ void Quantification_Model::model_calibrationcurve(std::vector