test_energy_calibration.cpp

/* SpecUtils: a library to parse, save, and manipulate gamma spectrum data files.
 
 Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC
 (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
 Government retains certain rights in this software.
 For questions contact William Johnson via email at wcjohns@sandia.gov, or
 alternative emails of interspec@sandia.gov.
 
 This library is free software; you can redistribute it and/or
 modify it under the terms of the GNU Lesser General Public
 License as published by the Free Software Foundation; either
 version 2.1 of the License, or (at your option) any later version.
 
 This library is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public
 License along with this library; if not, write to the Free Software
 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include <random>
#include <string>
#include <vector>
#include <limits>
#include <sstream>
#include <ostream>

#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
#include "doctest.h"

#include "SpecUtils/SpecFile.h"
#include "SpecUtils/EnergyCalibration.h"

#ifdef _WIN32
#undef min
#undef max
#endif

using namespace std;
using namespace SpecUtils;

string print_vec( const vector<float> &info )
{
  stringstream out;
  out << "{";
  for( size_t i = 0; i < info.size(); ++i )
    out << (i?",":"") << info[i];
  out << "}";
  return out.str();
}


bool is_similar( const vector<float> &lhs, const vector<float> &rhs )
{
  const size_t minsize = std::min( lhs.size(), rhs.size() );
  for( size_t i = 0; i < minsize; ++i )
  {
    const float larger = std::max(fabs(lhs[i]),fabs(rhs[i]));
    if( fabs(lhs[i]-rhs[i]) > 1.0E-5 * larger )
      return false;
  }
  
  for( size_t i = minsize; i < lhs.size(); ++i )
    if( fabs(lhs[i]) > numeric_limits<float>::epsilon() )
      return false;
  for( size_t i = minsize; i < rhs.size(); ++i )
    if( fabs(rhs[i]) > numeric_limits<float>::epsilon() )
      return false;
  
  return true;
}//bool is_similar( const vector<float> &lhs, const vector<float> &rhs )


TEST_CASE( "testCalibration" )
{
  using namespace SpecUtils;

  const size_t nbin = 1024;
  vector<float> frf_coefs;
  frf_coefs.push_back( 0.0f );
  frf_coefs.push_back( 3072.0f );
  frf_coefs.push_back( 0.0f );
  
  vector<float> poly_coefs = fullrangefraction_coef_to_polynomial( frf_coefs, nbin );
  vector<float> new_frf_coefs = polynomial_coef_to_fullrangefraction( poly_coefs, nbin );
  
  CHECK_MESSAGE( is_similar( frf_coefs, new_frf_coefs ), \
                       "Full Width Fraction coefficnets didnt make round trip: " \
                       << print_vec(frf_coefs) << "--->" << print_vec(new_frf_coefs) );

  vector<pair<float,float>> dev_pairs;
  auto frf_binning = fullrangefraction_binning( frf_coefs, nbin, dev_pairs );
  
  REQUIRE_MESSAGE( !!frf_binning, \
                         "Failed to make Full Width Fraction Binning for " \
                         << print_vec(frf_coefs) );
  REQUIRE_MESSAGE( frf_binning->size() == nbin, \
                         "Full range fraction returned " << frf_binning->size() \
                         << " intead of the expected " << nbin );
  
  for( size_t i = 0; i < nbin; ++i )
  {
    const float bin = static_cast<float>(i);
    const float lowerbinenergy = fullrangefraction_energy( bin, frf_coefs, nbin, dev_pairs );
    const float expected = frf_binning->at(i);
    const float larger = std::max(fabs(lowerbinenergy),fabs(expected));
    REQUIRE_MESSAGE( fabs(lowerbinenergy-expected) <= 1.0E-5*larger, \
                           "fullrangefraction_energy disagreed with fullrangefraction_binning" \
                           " (starting) at bin " << i \
                           << " got " << lowerbinenergy << " and " << expected \
                           << "respectively  for coefs=" << print_vec(frf_coefs) );
  }
  
  
  auto poly_binning = polynomial_binning( poly_coefs, nbin, dev_pairs );
  REQUIRE_MESSAGE( !!poly_binning, \
                        "Failed to make Polynomial Binning for " \
                        << print_vec(poly_coefs) );
  
  REQUIRE_MESSAGE( poly_binning->size() == nbin, \
                        "Polynomial binning returned " << poly_binning->size() \
                        << " intead of the expected " << nbin );
  
  for( size_t i = 0; i < nbin; ++i )
  {
    //const float bin = static_cast<float>(i);
    const float fwf = frf_binning->at(i);
    const float poly_eqn_energy = polynomial_energy( i, poly_coefs, dev_pairs );
    const float poly = poly_binning->at(i);
    const float larger = std::max(fabs(fwf),fabs(poly));
    
    REQUIRE_MESSAGE( fabs(fwf-poly) <= 1.0E-5*larger, "" \
                           "Lower channel energies for FWF and Polynomial coefficnets arent equal" \
                           " (starting) at bin " << i \
                           << " got " << fwf << " and " << poly \
                           << "respectively  for coefs=" << print_vec(frf_coefs) \
                           << " giving values: " << fwf << " and " << poly << " respectively" );
    
    REQUIRE_MESSAGE( fabs(poly_eqn_energy-poly) <= 1.0E-5*larger, "" \
                          "Lower channel energy for polynomial_energy and Poly binning arent equal" \
                          " (starting) at bin " << i \
                          << " got " << poly_eqn_energy << " and " << poly \
                          << "respectively  for coefs=" << print_vec(frf_coefs) \
                          << " giving values: " << poly_eqn_energy << " and " << poly
                          << " respectively" );
  }
  //Need to do a lot more tests here...

//  Need to test SpecFile::calibrationIsValid(...)
  
//  TEST_MESSAGE( "Input Directory:" << indir );
//  CHECK( false );
}//TEST_CASE( testCalibration )

TEST_CASE( "testFullRangeFractionFindEnergy" )
{
  // \TODO: Further tests to make:
  //        - Add deviation pairs
  //        - Test channels less than zero or greater than nbin
  //        - Test with 64k channels
  //        - Test with invalid calibration
  //        - Test with 2, 4 and 5 coefficents
  using namespace SpecUtils;

  const size_t nbin = 1024;
  vector<float> fwf_coefs;
  fwf_coefs.push_back( -1.926107f );
  fwf_coefs.push_back( 3020.178f );
  fwf_coefs.push_back( -8.720629f );
  
  vector<pair<float,float>> dev_pairs;
  const float accuracy = 0.001f;
  
  float binnum;
  const float energies[] = { 1121.68f, 1450.87f, 1480.65f };
  const size_t nenergies = sizeof(energies)/sizeof(energies[0]);
  
  for( size_t i = 0; i < nenergies; ++i )
  {
    const float energy = energies[i];
    CHECK_NOTHROW( binnum = find_fullrangefraction_channel( energy, fwf_coefs, nbin, dev_pairs, accuracy ) );
    const float binenergy = SpecUtils::fullrangefraction_energy( binnum, fwf_coefs, nbin, dev_pairs );



    REQUIRE_MESSAGE( fabs(binenergy-energy) < 0.1, "Found bin " << binnum \
                           << " for energy " << energy \
                           << " but found bin actually cooresponds to " \
                           << binenergy << " keV" );
  }//for( size_t i = 0; i < nenergies; ++i )
  


}//TEST_CASE( testFullRangeFractionFindEnergy )


TEST_CASE( "testFullRangeFractionFindEnergyFiveCoeffs" )
{
  // Test case from a real-world issue: coefficients with 5 terms including the low-energy term.
  // With 5 FRF coefficients, the 5th term is C_4/(1+60*x) where x = channel/nbin.
  // This term adds energy at low channels (31.1 keV when x=0), making the minimum energy at
  // channel 0 equal to about 8.3 keV for this calibration.
  //
  // When searching for an energy below the valid range, the function now extrapolates
  // using just the gain term, returning a negative channel number.
  using namespace SpecUtils;

  const size_t nbin = 1024;
  const vector<float> fwf_coefs = { -22.7999992f, 3013.77002f, 828.0f, -173.0f, 31.1f };
  const vector<pair<float,float>> dev_pairs;
  const float accuracy = 0.001f;

  // Verify the energy range of this calibration
  const double energy_at_channel_0 = fullrangefraction_energy( 0.0, fwf_coefs, nbin, dev_pairs );
  const double energy_at_channel_nbin = fullrangefraction_energy( static_cast<double>(nbin), fwf_coefs, nbin, dev_pairs );

  // Energy at channel 0 should be about 8.3 keV (= -22.8 + 31.1 from the 5th coef term)
  CHECK( energy_at_channel_0 > 8.0 );
  CHECK( energy_at_channel_0 < 9.0 );
  
  // The search energy is below the minimum valid energy
  const double search_energy = 5.7976100417681184;
  CHECK( search_energy < energy_at_channel_0 );
  
  // The function should now return a negative channel (extrapolated using gain)
  // instead of throwing an exception
  double channel = 0;
  CHECK_NOTHROW( channel = find_fullrangefraction_channel( search_energy, fwf_coefs, nbin, dev_pairs, accuracy ) );
  
  // The returned channel should be negative (below channel 0)
  CHECK( channel < 0.0 );
  
  // The extrapolation uses: channel = (energy - energy_at_0) * nbin / C_1
  const double expected_channel = (search_energy - energy_at_channel_0) * nbin / fwf_coefs[1];
  CHECK( fabs(channel - expected_channel) < 0.001 );
  
  // Test searching for a valid energy (within the calibration range)
  const double valid_energy = 100.0;
  CHECK_NOTHROW( channel = find_fullrangefraction_channel( valid_energy, fwf_coefs, nbin, dev_pairs, accuracy ) );
  const double found_energy = fullrangefraction_energy( channel, fwf_coefs, nbin, dev_pairs );
  CHECK( fabs(found_energy - valid_energy) < 0.1 );
  
  // Test searching for an energy above the valid range
  const double high_energy = energy_at_channel_nbin + 100.0;
  CHECK_NOTHROW( channel = find_fullrangefraction_channel( high_energy, fwf_coefs, nbin, dev_pairs, accuracy ) );
  CHECK( channel > nbin ); // Should be above the last channel
  
  // The extrapolation uses: channel = nbin + (energy - energy_at_nbin) * nbin / C_1
  const double expected_high_channel = nbin + (high_energy - energy_at_channel_nbin) * nbin / fwf_coefs[1];
  CHECK( fabs(channel - expected_high_channel) < 0.001 );
}//TEST_CASE( testFullRangeFractionFindEnergyFiveCoeffs )


TEST_CASE( "testPolynomialFindEnergy" )
{
  // \TODO: Further tests to make:
  //        - Add more deviation pairs
  //        - Test channels less than zero or greater than nbin
  //        - Test with 64k channels
  //        - Test with invalid calibration
  //        - Test with 2, 4 and 5 coefficents
  const size_t nbin = 1024;
  
  vector<float> poly_coefs = { -1.926107f, 2.9493925f, -0.00000831663990020752f };
  vector<pair<float,float>> dev_pairs = { {0.0f,0.0f}, {1460.0f,-10.0f}, {2614.0f,0.0f} };
  const float accuracy = 0.001f;
  
  float binnum;
  const float energies[] = { -100.0f, -10.0f, 511.0f, 1121.68f, 1450.87f, 1480.65f, 60000.0f };
  
  for( const float energy : energies )
  {
    CHECK_NOTHROW( binnum = find_polynomial_channel( energy, poly_coefs, nbin, dev_pairs, accuracy ) );
    const float binenergy = polynomial_energy( binnum, poly_coefs, dev_pairs );
    
    //Note: this doesnt test the case for multiple solution that the wanted solution is returned,
    //      it just checks the solution is correct.
    
    REQUIRE_MESSAGE( fabs(binenergy-energy) < 0.1, "Found bin " << binnum \
                           << " for energy " << energy \
                           << " but found bin actually cooresponds to " \
                           << binenergy << " keV" );
  }//for( loop over energies )

}//TEST_CASE( testPolynomialFindEnergy )


TEST_CASE( "testPolynomialFindEnergyLinearSimple" )
{
  const float energies[] = { -100.1f, -10.0f, 511.005f, 1121.68f, 1450.87f, 1480.65f, 60000.0f };
  
  for( const float energy : energies )
  {
    float binnum;
    CHECK_NOTHROW( binnum = find_polynomial_channel( energy, {0.0f, 1.0f}, 1024, {}, 0.001f ) );
    REQUIRE_MESSAGE( fabs(binnum - energy) < 0.1, "Found bin " << binnum \
                           << " for energy " << energy \
                           << " but found bin actually cooresponds to " \
                           << binnum << " keV" );
  }//for( loop over energies )
}//TEST_CASE( testPolynomialFindEnergyLinearSimple )


TEST_CASE( "testPolynomialFindEnergyRand" )
{
  const size_t nbin = 1024;
  const vector<float> poly_coefs = { -10.0f, 3.0f, -1.0f/(4.0f*nbin) };
  const vector<pair<float,float>> dev_pairs = { {0.0f,0.0f}, {661.0f,-19.0f}, {1460.0f,-10.0f}, {2614.0f,0.0f} };
  const float accuracy = 0.001f;
  
  std::random_device rd;
  std::mt19937 mt(rd());
  std::uniform_real_distribution<float> dist(-4.0*nbin, 4.0*nbin);
  
  //We will loop over bin-nummbers and make sure find channels returns the same channel-number as
  //  wanted
  for( size_t i = 0; i < 10000; ++i )
  {
    const float channel = dist(mt);
    const float channel_energy = polynomial_energy( channel, poly_coefs, dev_pairs );
    float found_channel;
    CHECK_NOTHROW( found_channel = find_polynomial_channel( channel_energy, poly_coefs, nbin, dev_pairs, accuracy ) );
    REQUIRE_MESSAGE( fabs(channel - found_channel) < 0.01, "Found channel " << found_channel \
                           << " for channel_energy " << channel_energy \
                           << " but actually wanted channel " << channel );
  }//for( loop over energies )
}//TEST_CASE( testPolynomialFindEnergyRand )


TEST_CASE( "testEnergyCalibrationLowerChannel" )
{
  const size_t nbin = 1024;
  vector<float> lower_channel( nbin + 1, 0.0f );
  for( size_t i = 0; i <= nbin; ++i )
    lower_channel[i] = i;
  
  EnergyCalibration cal;
  CHECK_NOTHROW( cal.set_lower_channel_energy( nbin, lower_channel ) );
  
  CHECK_THROWS_AS( cal.channel_for_energy(nbin+2), std::exception );
  CHECK_THROWS_AS( cal.channel_for_energy( -1 ), std::exception );
  CHECK_THROWS_AS( cal.energy_for_channel( nbin+2 ), std::exception );
  CHECK_THROWS_AS( cal.energy_for_channel( -1 ), std::exception );
  
  std::random_device rd;
  std::mt19937 mt(rd());
  std::uniform_real_distribution<float> dist( 0.0, nbin );
  
  for( size_t i = 0; i < 2000; ++i )
  {
    const float energy = dist(mt);
    
    float found_channel;
    CHECK_NOTHROW( found_channel = cal.channel_for_energy( energy ) );
    REQUIRE_MESSAGE( fabs(found_channel - energy) < 0.001, "Found channel " << found_channel \
                           << " for energy " << energy );
    
    const float channel = dist(mt);
    const float found_energy = cal.energy_for_channel( channel );
    REQUIRE_MESSAGE( fabs(found_energy - channel) < 0.001, "Found energy " << found_energy \
                           << " for channel " << channel );
  }//for( loop over energies )
  
}//TEST_CASE( testEnergyCalibrationLowerChannel )


TEST_CASE( "testCALpFile" )
{
  string calp_contents = R"(#PeakEasy CALp File Ver:  4.00
Offset (keV)           :  1.50000e+00
Gain (keV / Chan)      :  3.00000e+00
2nd Order Coef         :  0.00000e+00
3rd Order Coef         :  0.00000e+00
4th Order Coef         :  0.00000e+00
Deviation Pairs        :  5
7.70000e+01 -1.00000e+00
1.22000e+02 -5.00000e+00
2.39000e+02 -5.00000e+00
6.61000e+02 -2.90000e+01
2.61400e+03  0.00000e+00
#END)";
  
  stringstream input( calp_contents );


  string det_name;
  size_t num_channels = 1024;
  shared_ptr<EnergyCalibration> cal;
  
  CHECK_NOTHROW( cal = SpecUtils::energy_cal_from_CALp_file( input, num_channels, det_name ) );
  REQUIRE_MESSAGE( !!cal, "Failed to read basic CALp file" );
  
  CHECK( cal->valid() );
  CHECK_EQ( static_cast<int>(cal->type()), static_cast<int>(SpecUtils::EnergyCalType::Polynomial) );
  CHECK_EQ( cal->num_channels(), num_channels );
  CHECK_EQ( cal->deviation_pairs().size(), 5 );
  CHECK_EQ( cal->coefficients().size(), 2 );

  REQUIRE( cal->coefficients().size() >= 2 );
  CHECK_EQ( cal->coefficients()[0], 1.5 );
  CHECK_EQ( cal->coefficients()[1], 3.0 );

  REQUIRE( cal->deviation_pairs().size() == 5 );
  CHECK_EQ( cal->deviation_pairs()[0].first, 77.0f );
  CHECK_EQ( cal->deviation_pairs()[0].second, -1.0f );
  CHECK_EQ( cal->deviation_pairs()[1].first, 122.0f );
  CHECK_EQ( cal->deviation_pairs()[1].second, -5.0f );
  //...
  CHECK_EQ( cal->deviation_pairs()[4].first, 2614.0f );
  CHECK_EQ( cal->deviation_pairs()[4].second, 0.0f );

  // Test an invalid calibration specified in the CALp file
  calp_contents = R"(#PeakEasy CALp File Ver:  4.00
Offset (keV)           :  1.50000e+00
Gain (keV / Chan)      :  -3.00000e+00
2nd Order Coef         :  0.00000e+00
3rd Order Coef         :  0.00000e+00
4th Order Coef         :  0.00000e+00
#END)";
  input.str( calp_contents );
  cal = nullptr;
  CHECK_THROWS_AS( cal = SpecUtils::energy_cal_from_CALp_file( input, num_channels, det_name ), std::exception );
  CHECK( !cal );

  // Test an empty CALp file
  calp_contents = R"()";
  input.str( calp_contents );
  input.clear();
  CHECK_THROWS_AS( SpecUtils::energy_cal_from_CALp_file( input, num_channels, det_name ), std::exception );

  // TODO: add tests a multiple named detector demo, and then tests for `SpecFile::set_energy_calibration_from_CALp_file(...)`
  //       could/should also add tests for other calibration types

}//TEST_CASE( testCALpFile )


TEST_CASE( "fit_poly_energy_cal_from_points" )
{
  using namespace std;
  using namespace SpecUtils;

  SUBCASE( "Empty input throws" )
  {
    vector<pair<float,float>> empty_pairs;
    CHECK_THROWS_AS( fit_poly_energy_cal_from_points( empty_pairs, 2 ), std::exception );
  }

  SUBCASE( "Zero max_orders throws" )
  {
    vector<pair<float,float>> pairs = {{0.0f, 0.0f}, {100.0f, 300.0f}};
    CHECK_THROWS_AS( fit_poly_energy_cal_from_points( pairs, 0 ), std::exception );
  }

  SUBCASE( "max_orders exceeds number of points throws" )
  {
    vector<pair<float,float>> pairs = {{0.0f, 0.0f}, {100.0f, 300.0f}};
    CHECK_THROWS_AS( fit_poly_energy_cal_from_points( pairs, 3 ), std::exception );
  }

  SUBCASE( "Single point fit (linear through origin)" )
  {
    // Single point: channel=100, energy=300 keV
    // Should return [offset=0, gain=3.0]
    vector<pair<float,float>> pairs = {{100.0f, 300.0f}};

    vector<float> coeffs;
    CHECK_NOTHROW( coeffs = fit_poly_energy_cal_from_points( pairs, 1 ) );

    REQUIRE_EQ( coeffs.size(), 2 );
    CHECK_EQ( coeffs[0], doctest::Approx(0.0f).epsilon(0.0001) );  // offset
    CHECK_EQ( coeffs[1], doctest::Approx(3.0f).epsilon(0.0001) );  // gain
  }

  SUBCASE( "Two point linear fit" )
  {
    // Two points: (0, 10) and (100, 310)
    // Should fit: E = 10 + 3*ch
    vector<pair<float,float>> pairs = {{0.0f, 10.0f}, {100.0f, 310.0f}};

    vector<float> coeffs;
    CHECK_NOTHROW( coeffs = fit_poly_energy_cal_from_points( pairs, 2 ) );

    REQUIRE_EQ( coeffs.size(), 2 );
    CHECK_EQ( coeffs[0], doctest::Approx(10.0f).epsilon(0.001) );  // offset
    CHECK_EQ( coeffs[1], doctest::Approx(3.0f).epsilon(0.001) );   // gain
  }

  SUBCASE( "Three point linear fit (overdetermined)" )
  {
    // Three points on a line: E = 5 + 2*ch
    vector<pair<float,float>> pairs = {
      {0.0f, 5.0f},
      {50.0f, 105.0f},
      {100.0f, 205.0f}
    };

    vector<float> coeffs;
    CHECK_NOTHROW( coeffs = fit_poly_energy_cal_from_points( pairs, 2 ) );

    REQUIRE_EQ( coeffs.size(), 2 );
    CHECK_EQ( coeffs[0], doctest::Approx(5.0f).epsilon(0.001) );   // offset
    CHECK_EQ( coeffs[1], doctest::Approx(2.0f).epsilon(0.001) );   // gain
  }

  SUBCASE( "Three point quadratic fit (exact)" )
  {
    // Three points on parabola: E = 1 + 2*ch + 0.01*ch^2
    // At ch=0: E=1, ch=50: E=126, ch=100: E=301
    vector<pair<float,float>> pairs = {
      {0.0f, 1.0f},
      {50.0f, 126.0f},
      {100.0f, 301.0f}
    };

    vector<float> coeffs;
    CHECK_NOTHROW( coeffs = fit_poly_energy_cal_from_points( pairs, 3 ) );

    REQUIRE_EQ( coeffs.size(), 3 );
    CHECK_EQ( coeffs[0], doctest::Approx(1.0f).epsilon(0.001) );    // offset
    CHECK_EQ( coeffs[1], doctest::Approx(2.0f).epsilon(0.001) );    // linear
    CHECK_EQ( coeffs[2], doctest::Approx(0.01f).epsilon(0.001) );   // quadratic
  }

  SUBCASE( "Five point quadratic fit (overdetermined)" )
  {
    // Five points approximately on: E = 10 + 3*ch + 0.001*ch^2
    vector<pair<float,float>> pairs = {
      {0.0f, 10.0f},
      {100.0f, 320.0f},
      {200.0f, 650.0f},
      {300.0f, 1000.0f},
      {400.0f, 1370.0f}
    };

    vector<float> coeffs;
    CHECK_NOTHROW( coeffs = fit_poly_energy_cal_from_points( pairs, 3 ) );

    REQUIRE_EQ( coeffs.size(), 3 );
    CHECK_EQ( coeffs[0], doctest::Approx(10.0f).epsilon(0.01) );    // offset
    CHECK_EQ( coeffs[1], doctest::Approx(3.0f).epsilon(0.01) );     // linear
    CHECK_EQ( coeffs[2], doctest::Approx(0.001f).epsilon(0.0001) ); // quadratic
  }

  SUBCASE( "Typical detector calibration" )
  {
    // Realistic gamma calibration points from common peaks
    // Ba-133 @ 81 keV, Co-60 @ 1173 and 1333 keV
    vector<pair<float,float>> pairs = {
      {26.8f, 81.0f},      // Ba-133 81 keV
      {388.5f, 1173.5f},   // Co-60 1173 keV
      {441.5f, 1332.5f}    // Co-60 1333 keV
    };

    vector<float> coeffs;
    CHECK_NOTHROW( coeffs = fit_poly_energy_cal_from_points( pairs, 2 ) );

    REQUIRE_EQ( coeffs.size(), 2 );

    // Verify calibration is reasonable for a typical NaI detector
    CHECK( coeffs[0] > -10.0f );   // offset around 0 keV
    CHECK( coeffs[0] < 10.0f );
    CHECK( coeffs[1] > 2.5f );     // gain around 3 keV/channel
    CHECK( coeffs[1] < 3.5f );

    // Verify the fit reproduces the input points reasonably
    for( const auto &pair : pairs )
    {
      const float ch = pair.first;
      const float energy_expected = pair.second;
      const float energy_fitted = coeffs[0] + coeffs[1] * ch;
      CHECK_EQ( energy_fitted, doctest::Approx(energy_expected).epsilon(0.01) );
    }
  }

  SUBCASE( "Fractional channel numbers" )
  {
    // Test with non-integer channel numbers (as might come from peak fitting)
    vector<pair<float,float>> pairs = {
      {10.5f, 31.5f},
      {50.25f, 150.75f},
      {100.75f, 302.25f}
    };

    vector<float> coeffs;
    CHECK_NOTHROW( coeffs = fit_poly_energy_cal_from_points( pairs, 2 ) );

    REQUIRE_EQ( coeffs.size(), 2 );
    CHECK_EQ( coeffs[1], doctest::Approx(3.0f).epsilon(0.01) );  // gain = 3
  }

  SUBCASE( "All zeros throws (singular matrix)" )
  {
    vector<pair<float,float>> pairs = {
      {0.0f, 0.0f},
      {0.0f, 0.0f},
      {0.0f, 0.0f}
    };

    CHECK_THROWS_AS( fit_poly_energy_cal_from_points( pairs, 2 ), std::exception );
  }

  SUBCASE( "Negative gain detection" )
  {
    // Points that would result in negative gain
    vector<pair<float,float>> pairs = {
      {0.0f, 300.0f},
      {100.0f, 0.0f}
    };

    vector<float> coeffs;
    // Function should not throw, but PHD parser validation will catch negative gain
    CHECK_NOTHROW( coeffs = fit_poly_energy_cal_from_points( pairs, 2 ) );
    REQUIRE_EQ( coeffs.size(), 2 );
    CHECK( coeffs[1] < 0.0f );  // Verify we got negative gain
  }
}//TEST_CASE( "fit_poly_energy_cal_from_points" )


TEST_CASE( "testCachedDeviationPairSplines" )
{
  // Verify that channel_for_energy and energy_for_channel produce correct roundtrip results
  // when deviation pairs are present, and that the cached-spline member functions give the
  // same results as the free functions that build splines on the fly.
  using namespace SpecUtils;

  const size_t nbin = 1024;
  const vector<float> poly_coefs = { -1.926107f, 2.9493925f, -0.00000831663990020752f };
  const vector<pair<float,float>> dev_pairs = { {0.0f, 0.0f}, {1460.0f, -10.0f}, {2614.0f, 0.0f} };
  const float accuracy = 0.001f;

  EnergyCalibration cal;
  cal.set_polynomial( nbin, poly_coefs, dev_pairs );

  // Test roundtrip: energy -> channel -> energy should be close to original
  const vector<double> test_energies = { 10.0, 100.0, 500.0, 1000.0, 1460.0, 2000.0, 2600.0 };

  for( const double energy : test_energies )
  {
    const double channel = cal.channel_for_energy( energy );
    const double roundtrip_energy = cal.energy_for_channel( channel );
    CHECK_MESSAGE( fabs( roundtrip_energy - energy ) < 0.01,
      "Roundtrip failed for energy " << energy << " keV: got " << roundtrip_energy );
  }

  // Compare cached member function results against free function results
  for( const double energy : test_energies )
  {
    const double cached_channel = cal.channel_for_energy( energy );
    const double free_channel = find_polynomial_channel( energy, poly_coefs, nbin, dev_pairs, accuracy );
    CHECK_MESSAGE( fabs( cached_channel - free_channel ) < 1.0e-6,
      "channel_for_energy mismatch at " << energy << " keV: cached=" << cached_channel
      << " free=" << free_channel );
  }

  // Compare energy_for_channel: cached vs free function
  const vector<double> test_channels = { 0.0, 100.0, 500.0, 800.0, 1023.0 };
  for( const double channel : test_channels )
  {
    const double cached_energy = cal.energy_for_channel( channel );
    const double free_energy = polynomial_energy( channel, poly_coefs, dev_pairs );
    CHECK_MESSAGE( fabs( cached_energy - free_energy ) < 1.0e-6,
      "energy_for_channel mismatch at channel " << channel << ": cached=" << cached_energy
      << " free=" << free_energy );
  }
}//TEST_CASE( "testCachedDeviationPairSplines" )


TEST_CASE( "testCachedSplineUpdateOnSetPolynomial" )
{
  // Verify that cached splines are properly updated when setting new calibrations.
  using namespace SpecUtils;

  const size_t nbin = 1024;
  EnergyCalibration cal;

  // Set polynomial with deviation pairs
  const vector<float> coefs1 = { 0.0f, 3.0f };
  const vector<pair<float,float>> dev1 = { {0.0f, 0.0f}, {1500.0f, 5.0f}, {3000.0f, 0.0f} };
  cal.set_polynomial( nbin, coefs1, dev1 );

  const double ch1 = cal.channel_for_energy( 1500.0 );
  const double e1 = cal.energy_for_channel( 500.0 );

  // Set new polynomial with different deviation pairs
  const vector<float> coefs2 = { 0.0f, 3.0f };
  const vector<pair<float,float>> dev2 = { {0.0f, 0.0f}, {1500.0f, -5.0f}, {3000.0f, 0.0f} };
  cal.set_polynomial( nbin, coefs2, dev2 );

  const double ch2 = cal.channel_for_energy( 1500.0 );
  const double e2 = cal.energy_for_channel( 500.0 );

  // Results should differ since deviation pairs changed
  CHECK( fabs( ch1 - ch2 ) > 0.1 );
  CHECK( fabs( e1 - e2 ) > 0.1 );

  // Set polynomial with empty deviation pairs
  const vector<pair<float,float>> no_dev;
  cal.set_polynomial( nbin, coefs2, no_dev );

  const double ch3 = cal.channel_for_energy( 1500.0 );
  const double e3 = cal.energy_for_channel( 500.0 );

  // Without deviation pairs, energy_for_channel should be a simple polynomial
  CHECK( fabs( e3 - 1500.0 ) < 0.001 ); // 0.0 + 3.0*500 = 1500
  CHECK( fabs( ch3 - 500.0 ) < 0.001 ); // inverse: (1500 - 0) / 3.0 = 500
}//TEST_CASE( "testCachedSplineUpdateOnSetPolynomial" )


TEST_CASE( "testCachedSplinesFRF" )
{
  // Test the same caching behavior with Full Range Fraction calibration.
  using namespace SpecUtils;

  const size_t nbin = 1024;
  const vector<float> frf_coefs = { 0.0f, 3072.0f }; // gives 0 to 3072 keV
  const vector<pair<float,float>> dev_pairs = { {0.0f, 0.0f}, {1500.0f, 8.0f}, {3000.0f, 0.0f} };

  EnergyCalibration cal;
  cal.set_full_range_fraction( nbin, frf_coefs, dev_pairs );

  // Roundtrip test
  const vector<double> test_energies = { 50.0, 500.0, 1000.0, 1500.0, 2500.0, 3000.0 };
  for( const double energy : test_energies )
  {
    const double channel = cal.channel_for_energy( energy );
    const double roundtrip = cal.energy_for_channel( channel );
    CHECK_MESSAGE( fabs( roundtrip - energy ) < 0.01,
      "FRF roundtrip failed at " << energy << " keV: got " << roundtrip );
  }

  // Compare cached vs free function
  for( const double energy : test_energies )
  {
    const double cached = cal.channel_for_energy( energy );
    const double free_fn = find_fullrangefraction_channel( energy, frf_coefs, nbin, dev_pairs, 0.001 );
    CHECK_MESSAGE( fabs( cached - free_fn ) < 1.0e-6,
      "FRF channel_for_energy mismatch at " << energy << " keV" );
  }
}//TEST_CASE( "testCachedSplinesFRF" )


TEST_CASE( "testCachedSplinesLowerChannelEdge" )
{
  // Verify that set_lower_channel_energy clears splines and works correctly.
  using namespace SpecUtils;

  EnergyCalibration cal;

  // First set polynomial with deviation pairs
  const size_t nbin = 128;
  const vector<float> coefs = { 0.0f, 3.0f };
  const vector<pair<float,float>> dev_pairs = { {0.0f, 0.0f}, {200.0f, 5.0f}, {400.0f, 0.0f} };
  cal.set_polynomial( nbin, coefs, dev_pairs );

  // Now set lower channel edge (should clear deviation pair splines)
  vector<float> energies( nbin + 1 );
  for( size_t i = 0; i <= nbin; ++i )
    energies[i] = static_cast<float>( i * 3.0 );
  cal.set_lower_channel_energy( nbin, energies );

  CHECK( cal.type() == EnergyCalType::LowerChannelEdge );

  // channel_for_energy and energy_for_channel should work without issues
  const double ch = cal.channel_for_energy( 150.0 );
  CHECK( fabs( ch - 50.0 ) < 0.001 ); // 150 / 3.0 = 50

  const double e = cal.energy_for_channel( 50.0 );
  CHECK( fabs( e - 150.0 ) < 0.001 );
}//TEST_CASE( "testCachedSplinesLowerChannelEdge" )