Skip to content

Commit

Permalink
Fix energy cal not being correctly applied to HTML report of RelCalcA…
Browse files Browse the repository at this point in the history
…uto.
  • Loading branch information
wcjohns committed Dec 19, 2024
1 parent 3310bb4 commit febd7ea
Showing 1 changed file with 19 additions and 14 deletions.
33 changes: 19 additions & 14 deletions src/RelActCalcAuto.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<int>( 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<const SpecUtils::EnergyCalibration> new_cal = cost_functor->m_energy_cal;
shared_ptr<const SpecUtils::EnergyCalibration> new_cal = solution.m_foreground->energy_calibration();
if( success && options.fit_energy_cal )
{
auto new_cal = make_shared<SpecUtils::EnergyCalibration>( *cost_functor->m_energy_cal );
auto modified_new_cal = make_shared<SpecUtils::EnergyCalibration>( *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();

Expand All @@ -2320,27 +2322,30 @@ struct RelActAutoCostFcn /* : ROOT::Minuit2::FCNBase() */
{
vector<float> 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<float> 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<float> 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
Expand All @@ -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 )


Expand Down

0 comments on commit febd7ea

Please sign in to comment.