@@ -146,6 +146,7 @@ void Gaussian_Model<T_real>::set_fit_params_preset(Fit_Params_Preset preset)
146146
147147 switch (preset)
148148 {
149+ case Fit_Params_Preset::NOT_SET:
149150 case Fit_Params_Preset::MATRIX_BATCH_FIT: // matrix batch fit
150151 _fit_parameters[STR_ENERGY_OFFSET].bound_type = E_Bound_Type::FIXED;
151152 _fit_parameters[STR_ENERGY_SLOPE].bound_type = E_Bound_Type::FIXED;
@@ -580,7 +581,7 @@ const Spectra<T_real> Gaussian_Model<T_real>::model_spectrum_element(const Fit_P
580581 for (size_t idx = 0 ; idx < energy_ratios.size (); idx++)
581582 {
582583 const Element_Energy_Ratio<T_real>& er_struct = energy_ratios.at (idx);
583- 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 ) );
584+ 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)3.58 * fitp->at (STR_FWHM_FANOPRIME).value ) );
584585 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 )));
585586 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 ));
586587 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 ));
@@ -719,44 +720,16 @@ const ArrayTr<T_real> Gaussian_Model<T_real>::peak(T_real gain, T_real sigma, co
719720template <typename T_real>
720721const 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
721722{
722- /*
723- ArrayTr<T_real>counts = delta_energy.unaryExpr([gain, sigma, peak_E](T_real v)
724- {
725- T_real val = ((T_real)(M_SQRT2) * sigma);
726- val = std::erfc(v / val);
727- val = peak_E * val;
728- val = (T_real)2.0 / val;
729- return gain / val;
730- });
731- return counts;
732- */
733723 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)); });
734-
735- /*
736- ArrayTr<T_real>counts(delta_energy.size());
737- for (unsigned int i = 0; i < delta_energy.size(); i++)
738- {
739- counts[i] = gain / (T_real)2.0 / peak_E * std::erfc((T_real)delta_energy[i] / ((T_real)(M_SQRT2) * sigma));
740- }
741- */
742-
743724}
744725
745726// ----------------------------------------------------------------------------
746727
747728template <typename T_real>
748729const ArrayTr<T_real> Gaussian_Model<T_real>::tail(T_real gain, T_real sigma, const ArrayTr<T_real> &delta_energy, T_real gamma) const
749730{
750-
751- T_real val = pow (gamma, (T_real)2.0 );
752- val = exp ((T_real)-0.5 / val);
753- val = sigma / val;
754- val = gamma / val;
755- val = (T_real)2.0 / val;
756- val = gain / val;
757-
758731
759- // T_real val = gain / ( (T_real)2.0 * gamma * sigma * exp( (T_real)-0.5 / pow(gamma, (T_real)2.0) ) );
732+ T_real val = gain / ( (T_real)2.0 * gamma * sigma * exp ( (T_real)-0.5 / pow (gamma, (T_real)2.0 ) ) );
760733 return delta_energy.unaryExpr ([val,sigma,gamma](T_real v) { return (v < (T_real)0.0 ) ?
761734 std::exp (v/ (gamma * sigma)) * val * std::erfc (v / ( ((T_real)(M_SQRT2)*sigma) + ((T_real)1.0 /(gamma*(T_real)(M_SQRT2))) ) )
762735 :
@@ -771,7 +744,7 @@ const ArrayTr<T_real> Gaussian_Model<T_real>::elastic_peak(const Fit_Parameters<
771744{
772745 ArrayTr<T_real> counts (ev.size ());
773746 counts.setZero ();
774- 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 ) );
747+ 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)3.58 * fitp->at (STR_FWHM_FANOPRIME).value ) );
775748 if (false == std::isfinite (sigma))
776749 {
777750 logE << " sigma = " << sigma << " \n " ;
@@ -786,10 +759,7 @@ const ArrayTr<T_real> Gaussian_Model<T_real>::elastic_peak(const Fit_Parameters<
786759
787760 fvalue = fvalue * std::pow ((T_real)10.0 , fitp->at (STR_COHERENT_SCT_AMPLITUDE).value );
788761
789- // Spectra value = fvalue * this->peak(gain, *sigma, delta_energy);
790- // counts = counts + value;
791762 counts += ( fvalue * this ->peak (gain, sigma, delta_energy) );
792- // //counts += fvalue * (gain / ( sigma * (T_real)(SQRT_2xPI) ) * Eigen::exp((T_real)-0.5 * Eigen::pow((delta_energy / sigma), (T_real)2.0) ) );
793763
794764 return counts;
795765}
@@ -805,7 +775,7 @@ const ArrayTr<T_real> Gaussian_Model<T_real>::compton_peak(const Fit_Parameters<
805775
806776 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 )));
807777
808- 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 ) );
778+ 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)3.58 * fitp->at (STR_FWHM_FANOPRIME).value ) );
809779 if (false == std::isfinite (sigma))
810780 {
811781 logE << " sigma = " << sigma << " \n " ;
0 commit comments