From febd7ead9c1989d09438f30262e21b657eae7cc4 Mon Sep 17 00:00:00 2001 From: William Johnson Date: Thu, 19 Dec 2024 14:15:34 -0800 Subject: [PATCH] Fix energy cal not being correctly applied to HTML report of RelCalcAuto. --- src/RelActCalcAuto.cpp | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/src/RelActCalcAuto.cpp b/src/RelActCalcAuto.cpp index 33fb6508..ef2a331a 100644 --- a/src/RelActCalcAuto.cpp +++ b/src/RelActCalcAuto.cpp @@ -2297,18 +2297,20 @@ struct RelActAutoCostFcn /* : ROOT::Minuit2::FCNBase() */ }//if( we failed to get covariance ) / else }//if( solution.m_status == RelActCalcAuto::RelActAutoSolution::Status::Success ) - solution.m_num_function_eval_total = cost_functor->m_ncalls; + solution.m_num_function_eval_total = static_cast( cost_functor->m_ncalls ); solution.m_final_parameters = parameters; solution.m_energy_cal_adjustments[0] = parameters[0]; solution.m_energy_cal_adjustments[1] = parameters[1]; - shared_ptr new_cal = cost_functor->m_energy_cal; + shared_ptr new_cal = solution.m_foreground->energy_calibration(); if( success && options.fit_energy_cal ) { - auto new_cal = make_shared( *cost_functor->m_energy_cal ); + auto modified_new_cal = make_shared( *new_cal ); + const double offset = parameters[0]; + const double gain_mult = parameters[1]; const size_t num_channel = cost_functor->m_energy_cal->num_channels(); const SpecUtils::EnergyCalType energy_cal_type = cost_functor->m_energy_cal->type(); @@ -2320,27 +2322,30 @@ struct RelActAutoCostFcn /* : ROOT::Minuit2::FCNBase() */ { vector coefs = cost_functor->m_energy_cal->coefficients(); assert( coefs.size() >= 2 ); - coefs[0] += parameters[0]; - coefs[1] *= parameters[1]; + coefs[0] -= offset; + coefs[1] /= gain_mult; const auto &dev_pairs = cost_functor->m_energy_cal->deviation_pairs(); if( energy_cal_type == SpecUtils::EnergyCalType::FullRangeFraction ) - new_cal->set_full_range_fraction( num_channel, coefs, dev_pairs ); + modified_new_cal->set_full_range_fraction( num_channel, coefs, dev_pairs ); else - new_cal->set_polynomial( num_channel, coefs, dev_pairs ); + modified_new_cal->set_polynomial( num_channel, coefs, dev_pairs ); break; }//case polynomial or FRF - case SpecUtils::EnergyCalType::LowerChannelEdge: { - vector lower_energies = *cost_functor->m_energy_cal->channel_energies(); - for( float &energy : lower_energies ) - energy = parameters[0] + (parameters[1] * energy); + assert( new_cal->channel_energies() ); + if( !new_cal->channel_energies() || new_cal->channel_energies()->empty() ) + throw runtime_error( "Invalid lower channel energies???" ); + + vector lower_energies = *new_cal->channel_energies(); + for( float &x : lower_energies ) + x = (x - offset) / gain_mult; - new_cal->set_lower_channel_energy( num_channel, std::move(lower_energies) ); + modified_new_cal->set_lower_channel_energy( num_channel, std::move(lower_energies) ); break; }//case LowerChannelEdge @@ -2349,8 +2354,8 @@ struct RelActAutoCostFcn /* : ROOT::Minuit2::FCNBase() */ break; }//switch( m_energy_cal->type() ) - new_cal = new_cal; - solution.m_spectrum->set_energy_calibration( new_cal ); + new_cal = modified_new_cal; + solution.m_spectrum->set_energy_calibration( modified_new_cal ); }//if( options.fit_energy_cal )