diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.cc b/examples/ProgOpts/Test_Dpipi_ProgOpts.cc
index fe9739c..ab84e7b 100644
--- a/examples/ProgOpts/Test_Dpipi_ProgOpts.cc
+++ b/examples/ProgOpts/Test_Dpipi_ProgOpts.cc
@@ -1,148 +1,149 @@
 
 #include <iostream>
 
 #include "boost/program_options.hpp"
 
 #include "Test_Dpipi_ProgOpts.hh"
 
 namespace po = boost::program_options;
 
 void validate( boost::any& v, const std::vector<std::string>& values, Command* /*target_type*/, int)
 {
 	// Make sure no previous assignment to 'v' has been made
 	po::validators::check_first_occurrence(v);
 
 	// Extract the first string from 'values'. If there is more than
 	// one string, it's an error, and exception will be thrown.
 	const std::string& s { po::validators::get_single_string(values) };
 
 	if ( s == "gen" ) {
 		v = boost::any( Command::Generate );
 	} else if ( s == "fit" ) {
 		v = boost::any( Command::Fit );
 	} else if ( s == "simfit" ) {
 		v = boost::any( Command::SimFit );
 	} else {
 		throw po::validation_error(po::validation_error::invalid_option_value, "command", s, 3);
 	}
 }
 
 void validate( boost::any& v, const std::vector<std::string>& values, LauTimeDepFitModel::CPEigenvalue* /*target_type*/, int)
 {
 	// Make sure no previous assignment to 'v' has been made
 	po::validators::check_first_occurrence(v);
 
 	// Extract the first string from 'values'. If there is more than
 	// one string, it's an error, and exception will be thrown.
 	const std::string& s { po::validators::get_single_string(values) };
 
 	if ( s == "QFS" ) {
 		v = boost::any( LauTimeDepFitModel::CPEigenvalue::QFS );
 	} else if ( s == "CPEven" ) {
 		v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPEven );
 	} else if ( s == "CPOdd" ) {
 		v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPOdd );
 	} else {
 		throw po::validation_error(po::validation_error::invalid_option_value);
 	}
 }
 
 namespace LauDecayTime {
 void validate( boost::any& v, const std::vector<std::string>& values, EfficiencyMethod* /*target_type*/, int)
 {
 	// Make sure no previous assignment to 'v' has been made
 	po::validators::check_first_occurrence(v);
 
 	// Extract the first string from 'values'. If there is more than
 	// one string, it's an error, and exception will be thrown.
 	const std::string& s { po::validators::get_single_string(values) };
 
 	if ( s == "uniform" || s == "flat" ) {
 		v = boost::any( EfficiencyMethod::Uniform );
 	} else if ( s == "binned" || s == "hist" ) {
 		v = boost::any( EfficiencyMethod::Binned );
 	} else if ( s == "spline" ) {
 		v = boost::any( EfficiencyMethod::Spline );
 	} else {
 		throw po::validation_error(po::validation_error::invalid_option_value);
 	}
 }
 };
 
 TestDpipi_ProgramSettings::TestDpipi_ProgramSettings(const int argc, const  char ** argv)
 {
 	po::options_description main_desc{"Main options"};
 	main_desc.add_options()
 		("command", po::value<Command>(&command)->required(), "main command: gen, fit, or simfit")
 		;
 	po::positional_options_description p;
 	p.add("command", 1);
 
 	po::options_description common_desc{"Common options"};
 	common_desc.add_options()
 		("help", "produce help message")
 		("dtype", po::value<LauTimeDepFitModel::CPEigenvalue>(&dType)->default_value(LauTimeDepFitModel::CPEigenvalue::QFS,"QFS"), "type of D decay: QFS, CPOdd, or CPEven")
 		("dta-model", po::value<LauDecayTime::EfficiencyMethod>(&timeEffModel)->default_value(LauDecayTime::EfficiencyMethod::Uniform,"uniform"), "decay-time acceptance model: uniform/flat, binned/hist, spline")
 		("dtr", po::value<Bool_t>(&timeResolution)->default_value(kTRUE), "enable/disable decay-time resolution")
 		("dtr-perevent", po::value<Bool_t>(&perEventTimeErr)->default_value(kFALSE), "enable/disable use of per-event decay-time error (requires decay-time resolution to be enabled to take effect)")
 		("seed", po::value<UInt_t>(&RNGseed)->default_value(0), "set the seed for the RNG; if not set, the time is used to generate a seed")
 		("dir", po::value<std::string>(&directory)->default_value("","<none>"), "set the directory used to find the nTuples within the root file. Defaults to no directory")
 		;
 
 	po::options_description gen_desc{"Generation options"};
 	gen_desc.add_options()
 		("firstExptGen", po::value<UInt_t>(&firstExptGen)->default_value(0), "first experiment to generate")
 		("nExptGen", po::value<UInt_t>(&nExptGen)->default_value(1), "number of experiments to generate")
 		;
 
 	po::options_description fit_desc{"Fitting options"};
 	fit_desc.add_options()
 		("firstExpt", po::value<UInt_t>(&firstExptFit)->default_value(0), "first experiment to fit")
 		("nExpt", po::value<UInt_t>(&nExptFit)->default_value(1), "number of experiments to fit")
 		("iFit", po::value<UInt_t>(&iFit)->default_value(0), "fit ID - used to distinguish multiple fits to a given dataset (used to disentangle multiple minima)")
 		("fixTau", po::value<Bool_t>(&fixLifetime)->default_value(kTRUE), "enable/disable floating of B lifetime parameter")
 		("fixDm", po::value<Bool_t>(&fixDeltaM)->default_value(kTRUE), "enable/disable floating of B mixing frequency parameter")
 		("fixPhiMix", po::value<Bool_t>(&fixPhiMix)->default_value(kFALSE), "enable/disable floating of B mixing phase parameter(s)")
 		("fixSplineKnots", po::value<Bool_t>(&fixSplineKnots)->default_value(kFALSE), "enable/disable floating of the decay-time-acceptance spline")
 		("useSinCos", po::value<Bool_t>(&useSinCos)->default_value(kTRUE), "enable/disable using sine and cosine of B mixing phase as the floating parameters rather than the phase itself")
 		("useAveDeltaCalibVals", po::value<Bool_t>(&useAveDeltaCalibVals)->default_value(kTRUE), "enable/disable fitting calib values as their average +/- delta values rather than having separate values for B and Bbar")
+		("fitBack", po::value<Bool_t>(&fitBack)->default_value(kFALSE), "enable/disable refitting some generated data from a model")
 		;
 
 	po::options_description simfit_desc{"Simultaneous fitting options"};
 	simfit_desc.add_options()
 		("port", po::value<UInt_t>(&port)->default_value(0), "port number on which sim fit coordinator is listening")
 		;
 
 	po::options_description all_desc;
 	all_desc.add(main_desc).add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc);
 
 	po::variables_map vm;
 	try {
 		po::store(po::command_line_parser(argc, argv).
 				options(all_desc).
 				positional(p).
 				run(),
 				vm);
 		po::notify(vm);
 	}
 	catch ( boost::wrapexcept<po::required_option>& e ) {
 		std::cout << argv[0] << " requires a command, one of 'gen', 'fit', or 'simfit'\n"
 		          << "Append --help to those commands to see help on the related options" << std::endl;
 		parsedOK = kFALSE;
 		return;
 	}
 	catch ( po::validation_error& e ) {
 		std::cerr << e.what() << std::endl;
 		parsedOK = kFALSE;
 		return;
 	}
 
 	if ( vm.count("help") ) {
 		po::options_description help_desc;
 		help_desc.add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc);
 
 		std::cout << help_desc << std::endl;
 		helpRequested = kTRUE;
 		return;
 	}
 }
diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.hh b/examples/ProgOpts/Test_Dpipi_ProgOpts.hh
index 807d115..1462789 100644
--- a/examples/ProgOpts/Test_Dpipi_ProgOpts.hh
+++ b/examples/ProgOpts/Test_Dpipi_ProgOpts.hh
@@ -1,42 +1,43 @@
 #ifndef TEST_DPIPI_PROGOPTS
 #define TEST_DPIPI_PROGOPTS
 
 #include "Rtypes.h"
 
 #include "LauDecayTime.hh"
 #include "LauTimeDepFitModel.hh"
 
 #include "Command.hh"
 #include <string>
 
 struct TestDpipi_ProgramSettings {
 	TestDpipi_ProgramSettings(const int argc, const  char ** argv);
 
 	Bool_t parsedOK{kTRUE};
 	Bool_t helpRequested{kFALSE};
 
 	Command command{Command::Generate};
 	LauTimeDepFitModel::CPEigenvalue dType{LauTimeDepFitModel::CPEigenvalue::QFS};
 	LauDecayTime::EfficiencyMethod timeEffModel{LauDecayTime::EfficiencyMethod::Uniform};
 
 	UInt_t firstExptGen{0};
 	UInt_t nExptGen{1};
 	UInt_t firstExptFit{0};
 	UInt_t nExptFit{1};
 	UInt_t iFit{0};
 	UInt_t port{0};
 	UInt_t RNGseed{0};
 
 	Bool_t timeResolution{kTRUE};
 	Bool_t perEventTimeErr{kFALSE};
 	Bool_t fixLifetime{kTRUE};
 	Bool_t fixDeltaM{kTRUE};
 	Bool_t fixPhiMix{kFALSE};
 	Bool_t fixSplineKnots{kFALSE};
 	Bool_t useSinCos{kTRUE};
 	Bool_t useAveDeltaCalibVals{kTRUE};
 
 	std::string directory{""};
+	Bool_t fitBack{kFALSE};
 };
 
 #endif
diff --git a/inc/LauSplineDecayTimeEfficiency.hh b/inc/LauSplineDecayTimeEfficiency.hh
index 47be645..627991a 100644
--- a/inc/LauSplineDecayTimeEfficiency.hh
+++ b/inc/LauSplineDecayTimeEfficiency.hh
@@ -1,418 +1,423 @@
 
 /*
 Copyright 2021 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
 /*! \file LauSplineDecayTimeEfficiency.hh
     \brief File containing declaration of LauSplineDecayTimeEfficiency class.
 */
 
 /*! \class LauSplineDecayTimeEfficiency
     \brief Class for defining a spline-interpolated model for decay time efficiency
 */
 
 #ifndef LAU_SPLINE_DECAYTIME_EFFICIENCY
 #define LAU_SPLINE_DECAYTIME_EFFICIENCY
 
 #include <array>
 #include <memory>
 #include <vector>
 
 #include "Rtypes.h"
 #include "TFitResult.h"
 #include "TMath.h"
 #include "TSystem.h"
 
 #include "Lau1DCubicSpline.hh"
 #include "LauAbsDecayTimeEfficiency.hh"
 #include "LauParameter.hh"
 
 class TH1;
 
 //! Defines the allowed orders of spline polynomials
 enum class LauSplineOrder : std::size_t {
 	Cubic = 3, //!< 3rd order
 	Sixth = 6  //!< 6th order
 };
 
 template <LauSplineOrder Order>
 class LauSplineDecayTimeEfficiency : public LauAbsDecayTimeEfficiency {
 	public:
 		//! Constructor for a Cubic spline
 		/*!
 		    \param [in] prefix the prefix to use when forming the names of the knot parameters
 		    \param [in] effSpline the cubic spline to be used to represent the efficiency variation
 		*/
 		template <LauSplineOrder SpOrder = Order, typename std::enable_if_t<SpOrder == LauSplineOrder::Cubic,bool> = true>
 		LauSplineDecayTimeEfficiency( const TString& prefix, std::unique_ptr<Lau1DCubicSpline> effSpline ) :
 			prefix_{prefix},
 			effSpline_{std::move(effSpline)}
 		{
 			this->initialSetup();
 		}
 
 		//! Constructor for a Sixth order spline
 		/*!
 		    \param [in] prefix the prefix to use when forming the names of the knot parameters
 		    \param [in] effSpline the cubic spline being used to represent the efficiency variation in another channel
 		    \param [in] correctionSpline the cubic spline to be used to represent the correction to effSpline to obtain the efficiency variation for this channel
 		*/
 		template <LauSplineOrder SpOrder = Order, typename std::enable_if_t<SpOrder == LauSplineOrder::Sixth,bool> = true>
 		LauSplineDecayTimeEfficiency( const TString& prefix, const LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>& effSpline, std::unique_ptr<Lau1DCubicSpline> correctionSpline ) :
 			prefix_{prefix},
 			effSpline_{std::move(correctionSpline)},
 			otherSpline_{std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>>(effSpline)}
 		{
 			this->initialSetup();
 		}
 
 		//! Destructor
 		~LauSplineDecayTimeEfficiency() = default;
 		//! Move constructor (default)
 		LauSplineDecayTimeEfficiency(LauSplineDecayTimeEfficiency<Order>&& other) = default;
 		//! Move assignment operator (default)
 		LauSplineDecayTimeEfficiency& operator=(LauSplineDecayTimeEfficiency<Order>&& other) = default;
 		//! Copy constructor (custom)
 		LauSplineDecayTimeEfficiency(const LauSplineDecayTimeEfficiency<Order>& other);
 		//! Copy assignment operator (deleted)
 		LauSplineDecayTimeEfficiency& operator=(const LauSplineDecayTimeEfficiency<Order>& other) = delete;
 
 
 		//! The number of coefficients of each spline segment
 		static constexpr std::size_t nCoeffs { static_cast<std::underlying_type_t<LauSplineOrder>>(Order) + 1 };
 
 		//! Retrieve the efficiency for a given value of the decay time
 		/*!
 		    \param abscissa the value of the decay time
 		    \return the efficiency
 		*/
 		Double_t getEfficiency( const Double_t abscissa ) const override;
 
 		//! Retrieve the parameters of the efficiency model so that they can be loaded into a fit
 		/*!
 		    \return the parameters of the efficiency model
 		*/
 		std::vector<LauAbsRValue*> getParameters() override { return params_; }
 
 		//! Initialise the parameter cache
 		/*!
 			Must be called prior to starting fitting or generation
 		*/
 		void initialise() override;
 
 		//! Propagate any updates to parameters and recalculate information as neeeded
 		/*!
 			Should be called at each fit iteration
 		*/
 		void propagateParUpdates() override;
 
 		//! Retrieve the number of segments in the spline
 		std::size_t nSegments() const { return effSpline_->getnKnots() - 1; }
 
 		//! Retrieve the positions of the spline knots
 		const std::vector<Double_t>& knotPositions() const { return effSpline_->getXValues(); }
 
 		//! Fix all knots
 		void fixKnots();
 
 		//! Float all knots
 		void floatKnots();
 
 		//! Fix or float a specific knot
 		/*!
 		    \param knotIndex the index of the knot to fix/float
 		    \param fixed true if knot to be fixed, false if knot to be floated
 		*/
 		void fixKnot( const std::size_t knotIndex, const bool fixed );
 
 		//! Retrieve the polynomial coefficients for the given spline segment
 		/*!
 		    \param segmentIndex the index of the spline segment
 		    \return the polynomial coefficients
 		*/
 		std::array<Double_t,nCoeffs> getCoefficients( const std::size_t segmentIndex ) const;
 
 		//! Retrieve whether any of the parameters are floating in the fit
 		bool anythingFloating() const { return anyKnotFloating_; }
 
 		//! Retrieve whether any of the parameters have changed in the latest fit iteration
 		const bool& anythingChanged() const override { return anyKnotChanged_; }
 
 		//! Fit our spline to a TH1
 		/*!
 			Useful to obtain reasonable starting values and uncertainties for the knot parameters
 
 			\param [in] hist the histogram to be fit to
 			\param [in] printLevel the level of printout desired from fit
 		*/
 		void fitToTH1( TH1* hist, const LauOutputLevel printLevel = LauOutputLevel::Quiet );
 
 	private:
 		//! Perform the initial setup - shared by the various constructors
 		void initialSetup( const LauSplineDecayTimeEfficiency<Order>* other = nullptr );
 
 		//! Update the cached parameter values
 		void updateParameterCache();
 
 		//! The prefix for the parameter names
 		const TString prefix_;
 
 		//! The spline
 		std::unique_ptr<Lau1DCubicSpline> effSpline_;
 
 		//! The other spline
 		std::unique_ptr<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>> otherSpline_;
 
 		//! The spline parameters
 		std::vector<std::unique_ptr<LauParameter>> ownedParams_;
 
 		//! The spline parameters (for giving access)
 		std::vector<LauAbsRValue*> params_;
 
 		//! The spline values
 		std::vector<Double_t> values_;
 
 		//! Are any of our knot parameters floating in the fit?
 		bool ourKnotsFloating_{false};
 
 		//! Are any of the other spline's knot parameters floating in the fit?
 		bool otherKnotsFloating_{false};
 
 		//! Are any of the knot parameters floating in the fit?
 		bool anyKnotFloating_{false};
 
 		//! Have any of our knot parameters changed in the latest fit iteration?
 		bool ourKnotsChanged_{false};
 
 		//! Have any of the other spline's knot parameters changed in the latest fit iteration?
 		bool otherKnotsChanged_{false};
 
 		//! Have any of the knot parameters changed in the lastest fit iteration?
 		bool anyKnotChanged_{false};
 
 		ClassDefOverride(LauSplineDecayTimeEfficiency, 0)
 };
 
 templateClassImp(LauSplineDecayTimeEfficiency);
 
 template <LauSplineOrder Order>
 void LauSplineDecayTimeEfficiency<Order>::initialSetup( const LauSplineDecayTimeEfficiency<Order>* other )
 {
 	if ( not effSpline_ ) {
 		std::cerr<<"ERROR in LauSplineDecayTimeEfficiency::initialSetup : supplied efficiency spline pointer is null"<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	if constexpr ( Order == LauSplineOrder::Sixth ) {
 		std::vector<Double_t> ourKnots { effSpline_->getXValues() };
 		std::vector<Double_t> otherKnots { otherSpline_->knotPositions() };
 
 		if ( ourKnots != otherKnots ) {
 			std::cerr<<"ERROR in LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>::initialSetup : mismatch in knot positions"<<std::endl;
 			gSystem->Exit(EXIT_FAILURE);
 		}
 
 		// TODO - any other checks we need to do?
 	}
 
 	values_ = effSpline_->getYValues();
 
 	const std::size_t nPars { values_.size() };
 	ownedParams_.reserve( nPars );
 
 	if ( other ) {
 		for ( auto& ptr : other->ownedParams_ ) {
 			// TODO - maybe need a way to give these a unique prefix?
 			ownedParams_.emplace_back( ptr->createClone() );
 		}
 	} else {
 		for( std::size_t index{0}; index < nPars; ++index ) {
 			constexpr Double_t maxVal { (Order == LauSplineOrder::Cubic) ? 1.0 : 2.0 };
 			ownedParams_.emplace_back( std::make_unique<LauParameter>( Form( "%s_knot_%lu", prefix_.Data(), index ), values_[index], 0.0, maxVal, kTRUE ) );
 		}
 	}
 
 	if constexpr ( Order == LauSplineOrder::Sixth ) {
 		params_ = otherSpline_->getParameters();
 		params_.reserve( 2*nPars );
 	} else {
 		params_.reserve( nPars );
 	}
 
 	for ( auto& ptr : ownedParams_ ) {
 		params_.push_back( ptr.get() );
 	}
 }
 
 template <LauSplineOrder Order>
 LauSplineDecayTimeEfficiency<Order>::LauSplineDecayTimeEfficiency(const LauSplineDecayTimeEfficiency<Order>& other) :
 	prefix_{ other.prefix_ },
 	effSpline_{ std::make_unique<Lau1DCubicSpline>( * other.effSpline_ ) },
 	otherSpline_{ other.otherSpline_ ? std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>>( *other.otherSpline_ ) : nullptr }
 {
 	this->initialSetup( &other );
 }
 
 template <LauSplineOrder Order>
 void LauSplineDecayTimeEfficiency<Order>::fixKnots()
 {
 	for ( auto& ptr : ownedParams_ ) {
 		ptr->fixed(true);
 	}
 }
 
 template <LauSplineOrder Order>
 void LauSplineDecayTimeEfficiency<Order>::floatKnots()
 {
 	for ( auto& ptr : ownedParams_ ) {
 		ptr->fixed(false);
 	}
 }
 
 template <LauSplineOrder Order>
 void LauSplineDecayTimeEfficiency<Order>::fixKnot( const std::size_t knotIndex, const bool fixed )
 {
 	if ( knotIndex >= ownedParams_.size() ) {
 		std::cerr<<"ERROR in LauSplineDecayTimeEfficiency::fixKnot : supplied knot index is out of range"<<std::endl;
 		return;
 	}
 	ownedParams_[knotIndex]->fixed(fixed);
 }
 
 template <LauSplineOrder Order>
 Double_t LauSplineDecayTimeEfficiency<Order>::getEfficiency( const Double_t abscissa ) const
 {
 	Double_t eff { effSpline_->evaluate( abscissa ) };
 	if constexpr ( Order == LauSplineOrder::Sixth ) {
 		eff *= otherSpline_->getEfficiency( abscissa );
 	}
 	return eff;
 }
 
 template <LauSplineOrder Order>
 std::array<Double_t,LauSplineDecayTimeEfficiency<Order>::nCoeffs> LauSplineDecayTimeEfficiency<Order>::getCoefficients( const std::size_t segmentIndex ) const
 {
 	if constexpr ( Order == LauSplineOrder::Cubic ) {
 		return effSpline_->getCoefficients(segmentIndex);
 	} else {
 		std::array<Double_t,4> ourCoeffs { effSpline_->getCoefficients(segmentIndex) };
 		std::array<Double_t,4> otherCoeffs { otherSpline_->getCoefficients(segmentIndex) };
 		std::array<Double_t,nCoeffs> coeffs;
 		coeffs[0] = ourCoeffs[0] * otherCoeffs[0];
 		coeffs[1] = ourCoeffs[0] * otherCoeffs[1] + ourCoeffs[1] * otherCoeffs[0];
 		coeffs[2] = ourCoeffs[0] * otherCoeffs[2] + ourCoeffs[1] * otherCoeffs[1] + ourCoeffs[2] * otherCoeffs[0];
 		coeffs[3] = ourCoeffs[0] * otherCoeffs[3] + ourCoeffs[1] * otherCoeffs[2] + ourCoeffs[2] * otherCoeffs[1] + ourCoeffs[3] * otherCoeffs[0];
 		coeffs[4] = ourCoeffs[1] * otherCoeffs[3] + ourCoeffs[2] * otherCoeffs[2] + ourCoeffs[3] * otherCoeffs[1];
 		coeffs[5] = ourCoeffs[2] * otherCoeffs[3] + ourCoeffs[3] * otherCoeffs[2];
 		coeffs[6] = ourCoeffs[3] * otherCoeffs[3];
 		return coeffs;
 	}
 }
 
 template <LauSplineOrder Order>
 void LauSplineDecayTimeEfficiency<Order>::initialise()
 {
 	static auto isFixed = []( const std::unique_ptr<LauParameter>& par ){ return par->fixed(); };
 
 	anyKnotFloating_ = ourKnotsFloating_ = not std::all_of( ownedParams_.begin(), ownedParams_.end(), isFixed );
 
 	if constexpr ( Order == LauSplineOrder::Sixth ) {
 		otherSpline_->initialise();
 		otherKnotsFloating_ = otherSpline_->anythingFloating();
 		anyKnotFloating_ |= otherKnotsFloating_;
 	}
 
 	this->updateParameterCache();
 }
 
 template <LauSplineOrder Order>
 void LauSplineDecayTimeEfficiency<Order>::updateParameterCache()
 {
 	static auto assignValue = []( const std::unique_ptr<LauParameter>&  par ){ return par->unblindValue(); };
 
 	// Update our local cache
 	std::transform( ownedParams_.begin(), ownedParams_.end(), values_.begin(), assignValue );
 
 	// Update the spline itself
 	effSpline_->updateYValues( values_ );
 }
 
 template <LauSplineOrder Order>
 void LauSplineDecayTimeEfficiency<Order>::propagateParUpdates()
 {
 	// If none of the parameters are floating there's nothing to do
 	if ( not anyKnotFloating_ ) {
 		return;
 	}
 
 	// Otherwise, determine which of the floating parameters have changed (if any) and act accordingly
 
 	if ( ourKnotsFloating_ ) {
 		static auto checkEquality = []( const std::unique_ptr<LauParameter>& par, const Double_t cacheVal ){ return par->unblindValue() == cacheVal; };
 
 		anyKnotChanged_ = ourKnotsChanged_ = not std::equal( ownedParams_.begin(), ownedParams_.end(), values_.begin(), checkEquality);
 
 		// Update the acceptance spline if any of the knot values have changed
 		if ( ourKnotsChanged_ ) {
 			this->updateParameterCache();
 		}
 	}
 
 	if constexpr ( Order == LauSplineOrder::Sixth ) {
 		if ( otherKnotsFloating_ ) {
 			otherSpline_->propagateParUpdates();
 			otherKnotsChanged_ = otherSpline_->anythingChanged();
 			anyKnotChanged_ |= otherKnotsChanged_;
 		}
 	}
 }
 
 template <LauSplineOrder Order>
 void LauSplineDecayTimeEfficiency<Order>::fitToTH1( TH1* hist, const LauOutputLevel printLevel )
 {
 	std::cout << "INFO in LauSplineDecayTimeEfficiency::fitToTH1 : fitting " << prefix_ << " spline to supplied histogram \"" << hist->GetName() << "\"\n";
 	std::cout << "                                               : to obtain initial estimates of parameter values and errors" << std::endl;
 
 	TFitResultPtr fitResult { effSpline_->fitToTH1( hist, printLevel ) };
-	if ( fitResult->Status() != 0 || fitResult->CovMatrixStatus() != 3 ) {
-		std::cout << "WARNING in LauSplineDecayTimeEfficiency::fitToTH1 : fit failed with status " << fitResult->Status() << " and covariance matrix status " << fitResult->CovMatrixStatus() << "\n";
-		std::cout << "                                                  : will not update spline parameters based on this fit result" << std::endl;
+	if ( fitResult->Status() != 0 || fitResult->CovMatrixStatus() < 2 ) {
+		std::cerr << "WARNING in LauSplineDecayTimeEfficiency::fitToTH1 : fit failed with status " << fitResult->Status() << " and covariance matrix status " << fitResult->CovMatrixStatus() << "\n";
+		std::cerr << "                                                  : will not update spline parameters based on this fit result" << std::endl;
 		return;
 	}
+	if ( fitResult->CovMatrixStatus() == 2 )
+	{
+		std::cerr << "WARNING in LauSplineDecayTimeEfficiency::fitToTH1 : final covariance matrix not PosDef!" << std::endl;
+		std::cerr << "                                                  : continuing, but double check your fit!" << std::endl;
+	}
 
 	const std::size_t nKnots { ownedParams_.size() };
 	for ( std::size_t knot{0}; knot < nKnots; ++knot ) {
 		const Double_t knotVal { fitResult->Value(knot) };
 		const Double_t knotErr { fitResult->Error(knot) };
 
 		std::cout << "                                               : Setting parameter " << ownedParams_[knot]->name() << " to " << knotVal << " +/- " << knotErr << std::endl;
 
 		ownedParams_[knot]->value( knotVal );
 		ownedParams_[knot]->genValue( knotVal );
 		ownedParams_[knot]->initValue( knotVal );
 		ownedParams_[knot]->error( knotErr );
 	}
 
 	this->updateParameterCache();
 }
 
 #endif
diff --git a/src/Lau1DCubicSpline.cc b/src/Lau1DCubicSpline.cc
index 0cf449c..b3e62cb 100644
--- a/src/Lau1DCubicSpline.cc
+++ b/src/Lau1DCubicSpline.cc
@@ -1,565 +1,566 @@
 
 /*
 Copyright 2015 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
 /*! \file Lau1DCubicSpline.cc
   \brief File containing implementation of Lau1DCubicSpline class.
  */
 
 #include <cmath>
 #include <cstdlib>
 #include <iostream>
 
 #include <TH2.h>
 #include <TSystem.h>
 
 #include "Lau1DCubicSpline.hh"
 #include "LauAbsRValue.hh"
 #include "LauParameter.hh"
 
 ClassImp(Lau1DCubicSpline)
 
 
 Lau1DCubicSpline::Lau1DCubicSpline(const std::vector<Double_t>& xs, const std::vector<Double_t>& ys,
 				   LauSplineType type, 
 				   LauSplineBoundaryType leftBound, 
 				   LauSplineBoundaryType rightBound,
 				   Double_t dydx0, Double_t dydxn) :
 	nKnots_(xs.size()),
 	x_(xs),
 	y_(ys),
 	type_(type),
 	leftBound_(leftBound),
 	rightBound_(rightBound),
 	dydx0_(dydx0),
 	dydxn_(dydxn)
 {
 	init();
 }
 
 Double_t Lau1DCubicSpline::evaluate(Double_t x) const
 {
 	// do not attempt to extrapolate the spline
 	if( x<x_[0] || x>x_[nKnots_-1] ) {
 		std::cerr << "WARNING in Lau1DCubicSpline::evaluate : function is only defined between " << x_[0] << " and " << x_[nKnots_-1] << std::endl;
 		std::cerr << "                                        value at " << x << " returned as 0" << std::endl;
 		return 0.;
 	}
 
 	// first determine which 'cell' of the spline x is in
 	// cell i runs from knot i to knot i+1
 	Int_t cell(0);
 
 	while( x > x_[cell+1] ) {
 		++cell;
 	}
 
 	// obtain x- and y-values of the neighbouring knots
 	Double_t xLow  = x_[cell];
 	Double_t xHigh = x_[cell+1];
 	Double_t yLow  = y_[cell];
 	Double_t yHigh = y_[cell+1];
 
 	if(type_ == Lau1DCubicSpline::LinearInterpolation) {
 		return yHigh*(x-xLow)/(xHigh-xLow) + yLow*(xHigh-x)/(xHigh-xLow);
 	}
 
 	// obtain t, the normalised x-coordinate within the cell, 
 	// and the coefficients a and b, which are defined in cell i as:
 	//
 	//      a_i =  k_i  *(x_i+1 - x_i) - (y_i+1 - y_i),
 	//      b_i = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i)
 	//
 	// where k_i is (by construction) the first derivative at knot i
 	Double_t t = (x - xLow) / (xHigh - xLow);
 	Double_t a =     dydx_[cell]   * (xHigh - xLow) - (yHigh - yLow);
 	Double_t b = -1.*dydx_[cell+1] * (xHigh - xLow) + (yHigh - yLow);
 
 	Double_t retVal = (1 - t) * yLow + t * yHigh + t * (1 - t) * ( a * (1 - t) + b * t );
 
 	return retVal;
 }
 
 void Lau1DCubicSpline::updateYValues(const std::vector<Double_t>& ys)
 {
 	y_ = ys;
 	this->calcDerivatives();
 }
 
 void Lau1DCubicSpline::updateYValues(const std::vector<LauParameter*>& ys)
 {
 	for (std::size_t i=0; i<nKnots_; ++i){
 		y_[i] = ys[i]->unblindValue();
 	}
 	this->calcDerivatives();
 }
 
 void Lau1DCubicSpline::updateYValues(const std::vector<LauAbsRValue*>& ys)
 {
 	for (std::size_t i=0; i<nKnots_; ++i){
 		y_[i] = ys[i]->unblindValue();
 	}
 	this->calcDerivatives();
 }
 
 void Lau1DCubicSpline::updateType(LauSplineType type)
 {
 	if(type_ != type) {
 		type_ = type;
 		this->calcDerivatives();
 	}
 }
 
 void Lau1DCubicSpline::updateBoundaryConditions(LauSplineBoundaryType leftBound,
 		LauSplineBoundaryType rightBound,
 		Double_t dydx0, Double_t dydxn)
 {
 	bool updateDerivatives{false};
 
 	if(leftBound_ != leftBound || rightBound_ != rightBound) {
 		leftBound_ = leftBound;
 		rightBound_ = rightBound;
 		updateDerivatives = true;
 	}
 
 	if(dydx0_ != dydx0) {
 		dydx0_ = dydx0;
 		if(leftBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = true;
 	}
 
 	if(dydxn_ != dydxn) {
 		dydxn_ = dydxn;
 		if(rightBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = true;
 	}
 
 	if(updateDerivatives) {
 		this->calcDerivatives();
 	}
 }
 
 std::array<Double_t,4> Lau1DCubicSpline::getCoefficients(const std::size_t segIndex, const bool normalise) const
 {
 	std::array<Double_t,4> result = {0.,0.,0.,0.};
 	
 	if(segIndex >= nKnots_-1)
 	{
 		std::cerr << "ERROR in Lau1DCubicSpline::getCoefficients requested for too high a knot value" << std::endl;
 		return result;
 	}
 	
 	Double_t xL = x_[segIndex]   , xH = x_[segIndex+1];
 	Double_t yL = y_[segIndex]   , yH = y_[segIndex+1];
 	Double_t h = xH-xL; //This number comes up a lot
 	
 	switch(type_)
 	{
 		case Lau1DCubicSpline::StandardSpline:
 		case Lau1DCubicSpline::AkimaSpline:
 			{
 				Double_t kL = dydx_[segIndex], kH = dydx_[segIndex+1];
 
 				//a and b based on definitions from https://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline
 				Double_t a = kL*h-(yH-yL);
 				Double_t b =-kH*h+(yH-yL);
 
 				Double_t denom = -h*h*h;//The terms have a common demoninator
 
 				result[0] = -b*xL*xL*xH + a*xL*xH*xH + h*h*(xL*yH - xH*yL);
 				result[1] = -a*xH*(2*xL+xH) + b*xL*(xL+2*xH) + h*h*(yL-yH);
 				result[2] = -b*(2*xL+xH) + a*(xL+2*xH);
 				result[3] = -a+b;
 
 				for(auto& res : result){res /= denom;}
 				break;
 			}
 /*		case Lau1DCubicSpline::AkimaSpline: // Double check the Akima description of splines (in evaluate) right now they're the same except for the first derivatives
 			{
 				//using fomulae from https://asmquantmacro.com/2015/09/01/akima-spline-interpolation-in-excel/
 				std::function<Double_t(Int_t)> m = [&](Int_t j) //formula to get the straight line gradient
 				{
 					if(j < 0){return 2*m(j+1)-m(j+2);}
 					if(j >= nKnots_){return 2*m(j-1)-m(j-2);}
 					return (y_[j+1]-y_[j]) / (x_[j+1]-x_[j]);
 				};
 				
 				auto t = [&](Int_t j)
 				{
 					Double_t res = 0.; //originally res was called 't' but that produced a shadow warning
 					Double_t denom = TMath::Abs(m(j+1)-m(j)) + TMath::Abs(m(j-1)-m(j-2));
 					if(denom == 0){res = (m(j)-m(j-1))/2;} //Special case for when denom = 0
 					else
 					{
 						res = TMath::Abs(m(j+1)-m(j))*m(j-1) + TMath::Abs(m(j-1)-m(j-2))*m(j);
 						res /= denom;
 					}
 					return res;
 				};
 
 				//These are the p's to get the spline in the form p_k(x-xL)^k
 				Double_t pDenom = x_[segIndex+1]/x_[segIndex]; //a denominator used for p[2] and p[3]
 				std::array<Double_t,4> p = {y_[segIndex],t(segIndex),0.,0.}; //we'll do the last 2 below
 				
 				p[2] = 3*m(segIndex)-2*t(segIndex)-t(segIndex+1);
 				p[2]/= pDenom;
 
 				p[3] = t(segIndex)+t(segIndex+1)-2*m(segIndex);
 				p[3]/= pDenom*pDenom;
 
 				//Now finally rearranging the p's into the desired results
 				result[0] = p[0]-p[1]*xL+p[2]*xL*xL-p[3]*xL*xL*xL;
 				result[1] = p[1]-2*p[2]*xL+3*p[3]*xL*xL;
 				result[2] = p[2]-3*p[3]*xL;
 				result[3] = p[3];
 
 				break;
 			}*/
 		case Lau1DCubicSpline::LinearInterpolation:
 			{
 				result[0] = xH*yL-xL*yH;
 				result[1] = yH-yL;
 
 				for(auto& res : result){res /= h;}
 				break;
 			}
 	}
 
 	if(normalise)
 	{
 		Double_t integral = this->integral();
 		for(auto& res : result){res /= integral;}
 	}
 
 	return result;
 }
 
 Double_t Lau1DCubicSpline::integral() const
 {
 	Double_t integral{0.0};
 	for(std::size_t iKnot{0}; iKnot < nKnots_ -1; ++iKnot)
 	{
 		const Double_t minAbs { x_[iKnot] };
 		const Double_t maxAbs { x_[iKnot+1] };
 
 		const std::array<Double_t,4> coeffs { this->getCoefficients(iKnot, false) };
 
 		auto integralFunc = [&coeffs](Double_t x){return coeffs[0]*x + coeffs[1]*x*x/2.0 + coeffs[2]*x*x*x/3.0 + coeffs[3]*x*x*x*x/4.0;};
 
 		integral += integralFunc(maxAbs);
 		integral -= integralFunc(minAbs);
 	}
 	return integral;
 }
 
 TF1* Lau1DCubicSpline::makeTF1(const bool normalise) const
 {
 	TString functionString{""}; //make a long piecewise construction of all the spline pieces
 	for(std::size_t segIndex{0}; segIndex < nKnots_-1; ++segIndex)
 	{
 		const std::array<Double_t,4> coeffs { this->getCoefficients(segIndex,normalise) };
 		functionString += Form("(x>%f && x<= %f)*",x_[segIndex],x_[segIndex+1]);//get the bounds of this piece
 		functionString += Form("(%f + %f*x + %f*x^2 + %f*x^3)",coeffs[0],coeffs[1],coeffs[2],coeffs[3]);
 		if(segIndex < nKnots_ -2){functionString += " + \n";}//add to all lines except the last
 	}
 
 	TF1* func = new TF1("SplineFunction", functionString, x_.front(), x_.back());
 
 	return func;
 }
 
 TF1* Lau1DCubicSpline::normaliseToTH1(TH1* hist, LauOutputLevel printLevel) const
 {
 	//first define the fit function
 	auto fitf = [this](Double_t* x, Double_t* par) {//there is only 1 x (the abscissa) and 1 par (a scaling of the entire thing)
 		return this->evaluate( x[0] ) * par[0];
 	};
 	//Make the function
 	TF1* FittedFunc = new TF1("FittedFunc",fitf,x_.front(),x_.back(),1);
 	//Set the parameter name
 	FittedFunc -> SetParNames("Constant");
 	//Set the parameter limits and default value
 	FittedFunc -> SetParLimits(0,0.,10.);
 	FittedFunc -> SetParameter(0,1.);
 
 	// Options to:
 	// - N = do not store the graphics function
 	// - 0 = do not plot the fit result
 	TString fitOptions {"N 0"};
 
 	switch ( printLevel ) {
 		case LauOutputLevel::Verbose :
 			fitOptions += " V";
 			break;
 		case LauOutputLevel::Quiet :
 			fitOptions += " Q";
 			break;
 		case LauOutputLevel::Standard :
 			break;
 	}
 
 	hist->Fit( "FittedFunc", fitOptions );
 
 	return FittedFunc;
 }
 
 TFitResultPtr Lau1DCubicSpline::fitToTH1(TH1* hist, LauOutputLevel printLevel)
 {
 	auto fitf = [this](Double_t* x, Double_t* par) {
 		this -> updateYValues( std::vector<Double_t>(par, par + nKnots_) );
 		return this -> evaluate( x[0] );
 	};
 
 	TF1 FittedFunc("FittedFunc",fitf,x_.front(),x_.back(),nKnots_);
 
 	const Double_t knotMax { hist->GetMaximum() * 1.5 };
 	
 	for(std::size_t knot{0}; knot <= nKnots_ ; ++knot)
 	{
 		FittedFunc.SetParName(knot, Form("Knot%lu",knot));
 		FittedFunc.SetParLimits(knot, 0., knotMax);
 		FittedFunc.SetParameter(knot, y_[knot]);
 	}
 
 	// Options to:
 	// - N = do not store the graphics function
 	// - 0 = do not plot the fit result
 	// - S = save the result of the fit
-	TString fitOptions {"N 0 S"};
+	// - WL= do a weighted likelihood fit (not chi2)
+	TString fitOptions {"N 0 S WL"};
 
 	switch ( printLevel ) {
 		case LauOutputLevel::Verbose :
 			fitOptions += " V";
 			break;
 		case LauOutputLevel::Quiet :
 			fitOptions += " Q";
 			break;
 		case LauOutputLevel::Standard :
 			break;
 	}
 
 	return hist->Fit( "FittedFunc", fitOptions );
 }
 
 void Lau1DCubicSpline::init()
 {
 	if( y_.size() != x_.size()) {
 		std::cerr << "ERROR in Lau1DCubicSpline::init : The number of y-values given does not match the number of x-values" << std::endl;
 		std::cerr << "                                  Found " <<  y_.size() << ", expected " << x_.size() << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	if( nKnots_ < 3 ) {
 		std::cerr << "ERROR in Lau1DCubicSpline::init : The number of knots is too small" << std::endl;
 		std::cerr << "                                  Found " <<  nKnots_ << ", expected at least 3 (to have at least 1 internal knot)" << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	dydx_.assign(nKnots_,0.0);
 
 	a_.assign(nKnots_,0.0);
 	b_.assign(nKnots_,0.0);
 	c_.assign(nKnots_,0.0);
 	d_.assign(nKnots_,0.0);
 
 	this->calcDerivatives();
 }
 
 void Lau1DCubicSpline::calcDerivatives()
 {
 	switch ( type_ ) {
 		case Lau1DCubicSpline::StandardSpline :
 			this->calcDerivativesStandard();
 			break;
 		case Lau1DCubicSpline::AkimaSpline :
 			this->calcDerivativesAkima();
 			break;
 		case Lau1DCubicSpline::LinearInterpolation :
 			//derivatives not needed for linear interpolation
 			break;
 	}
 }
 
 void Lau1DCubicSpline::calcDerivativesStandard()
 {
 	// derivatives are determined such that the second derivative is continuous at internal knots
 
 	// derivatives, k_i, are the solutions to a set of linear equations of the form:
 	// a_i * k_i-1  +  b_i * k+i  +  c_i * k_i+1  =  d_i with a_0 = 0, c_n-1 = 0
 
 	// this is solved using the tridiagonal matrix algorithm as on en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
 
 	// first and last equations give boundary conditions
 	// - for natural boundary, require f''(x) = 0 at end knot
 	// - for 'not a knot' boundary, require f'''(x) continuous at second knot
 	// - for clamped boundary, require predefined value of f'(x) at end knot
 
 	// non-zero values of a_0 and c_n-1 would give cyclic boundary conditions
 	a_[0] = 0.;
 	c_[nKnots_-1] = 0.;
 
 	// set left boundary condition
 	if(leftBound_ == Lau1DCubicSpline::Natural) {
 		b_[0] = 2./(x_[1]-x_[0]);
 		c_[0] = 1./(x_[1]-x_[0]);
 		d_[0] = 3.*(y_[1]-y_[0])/((x_[1]-x_[0])*(x_[1]-x_[0]));
 
 	} else if(leftBound_ == Lau1DCubicSpline::NotAKnot) {
 		// define the width, h, and the 'slope', delta, of the first cell
 		Double_t h1(x_[1]-x_[0]), h2(x_[2]-x_[0]);
 		Double_t delta1((y_[1]-y_[0])/h1), delta2((y_[2]-y_[1])/h2);
 
 		// these coefficients can be determined by requiring f'''_0(x_1) = f'''_1(x_1)
 		// the requirement f''_0(x_1) = f''_1(x_1) has been used to remove the dependence on k_2
 		b_[0] = h2;
 		c_[0] = h1+h2;
 		d_[0] = delta1*(2.*h2*h2 + 3.*h1*h2)/(h1+h2) + delta2*5.*h1*h1/(h1+h2);
 
 	} else { //Clamped
 		b_[0] = 1.;
 		c_[0] = 0.;
 		d_[0] = dydx0_;
 	}
 
 	// set right boundary condition
 	if(rightBound_ == Lau1DCubicSpline::Natural) {
 		a_[nKnots_-1] = 1./(x_[nKnots_-1]-x_[nKnots_-2]);
 		b_[nKnots_-1] = 2./(x_[nKnots_-1]-x_[nKnots_-2]);
 		d_[nKnots_-1] = 3.*(y_[nKnots_-1]-y_[nKnots_-2])/((x_[nKnots_-1]-x_[nKnots_-2])*(x_[nKnots_-1]-x_[nKnots_-2]));
 
 	} else if(rightBound_ == Lau1DCubicSpline::NotAKnot) {
 		// define the width, h, and the 'slope', delta, of the last cell
 		Double_t hnm1(x_[nKnots_-1]-x_[nKnots_-2]), hnm2(x_[nKnots_-2]-x_[nKnots_-3]);
 		Double_t deltanm1((y_[nKnots_-1]-y_[nKnots_-2])/hnm1), deltanm2((y_[nKnots_-2]-y_[nKnots_-3])/hnm2);
 
 		// these coefficients can be determined by requiring f'''_n-3(x_n-2) = f'''_n-2(x_n-2)
 		// the requirement f''_n-3(x_n-2) = f''_n-2(x_n-2) has been used to remove 
 		// the dependence on k_n-3
 		a_[nKnots_-1] = hnm2 + hnm1;
 		b_[nKnots_-1] = hnm1;
 		d_[nKnots_-1] = deltanm2*hnm1*hnm1/(hnm2+hnm1) + deltanm1*(2.*hnm2*hnm2 + 3.*hnm2*hnm1)/(hnm2+hnm1);
 
 	} else { //Clamped
 		a_[nKnots_-1] = 0.;
 		b_[nKnots_-1] = 1.;
 		d_[nKnots_-1] = dydxn_;
 	}
 
 	// the remaining equations ensure that f_i-1''(x_i) = f''_i(x_i) for all internal knots
 	for(std::size_t i{1}; i<nKnots_-1; ++i) {
 		a_[i] = 1./(x_[i]-x_[i-1]);
 		b_[i] = 2./(x_[i]-x_[i-1]) + 2./(x_[i+1]-x_[i]);
 		c_[i] =                      1./(x_[i+1]-x_[i]);
 		d_[i] = 3.*(y_[i]-y_[i-1])/((x_[i]-x_[i-1])*(x_[i]-x_[i-1])) + 3.*(y_[i+1]-y_[i])/((x_[i+1]-x_[i])*(x_[i+1]-x_[i]));
 	}
 
 	// forward sweep - replace c_i and d_i with
 	//
 	// c'_i =  c_i /  b_i                 for i = 0
 	//         c_i / (b_i - a_i * c_i-1)  otherwise
 	//
 	// and 
 	//
 	// d'_i =   d_i                /  b_i                 for i = 0
 	//         (d_i - a_i * d_i-1) / (b_i - a_i * c_i-1)  otherwise
 
 	c_[0] /= b_[0];
 	d_[0] /= b_[0];
 
 	for(std::size_t i{1}; i<nKnots_-1; ++i) {
 		c_[i] =  c_[i]                  / (b_[i] - a_[i]*c_[i-1]);
 		d_[i] = (d_[i] - a_[i]*d_[i-1]) / (b_[i] - a_[i]*c_[i-1]);
 	}
 
 	d_[nKnots_-1] = (d_[nKnots_-1] - a_[nKnots_-1]*d_[nKnots_-2]) / (b_[nKnots_-1] - a_[nKnots_-1]*c_[nKnots_-2]);
 
 	// back substitution - calculate k_i = dy/dx at each knot
 	//
 	// k_i = d'_i                   for i = n-1
 	//       d'_i - c'_i * k_i+1    otherwise
 
 	dydx_[nKnots_-1] = d_[nKnots_-1];
 
 	for(int i {static_cast<int>(nKnots_-2)}; i>=0; --i) {
 		dydx_[i] = d_[i] - c_[i]*dydx_[i+1];
 	}
 }
 
 void Lau1DCubicSpline::calcDerivativesAkima() {
 	//derivatives are calculated according to the Akima method
 	// J.ACM vol. 17 no. 4 pp 589-602
 
 	Double_t am1(0.), an(0.), anp1(0.);
 
 	// a[i] is the slope of the segment from point i-1 to point i
 	//
 	// n.b. segment 0 is before point 0 and segment n is after point n-1
 	//      internal segments are numbered 1 - n-1
 	for(std::size_t i{1}; i<nKnots_; ++i) {
 		a_[i] = (y_[i]-y_[i-1])/(x_[i]-x_[i-1]);
 	}
 
 	// calculate slopes between additional points on each end according to the Akima method
 	// method assumes that the additional points follow a quadratic defined by the last three points
 	// this leads to the relations a[2] - a[1] = a[1] - a[0] = a[0] - a[-1] and a[n-1] - a[n-2] = a[n] - a[n-1] = a[n+1] - a[n]
 	a_[0]         = 2*a_[1]         - a_[2];
 	am1           = 2*a_[0]         - a_[1];
 
 	an            = 2*a_[nKnots_-1] - a_[nKnots_-2];
 	anp1          = 2*an            - a_[nKnots_-1];
 
 	// b[i] is the weight of a[i]   towards dydx[i]
 	// c[i] is the weight of a[i+1] towards dydx[i]
 	// See Appendix A of J.ACM vol. 17 no. 4 pp 589-602 for a justification of these weights
 	b_[0] = TMath::Abs(a_[1] - a_[0]);
 	b_[1] = TMath::Abs(a_[2] - a_[1]);
 
 	c_[0] = TMath::Abs(a_[0] - am1  );
 	c_[1] = TMath::Abs(a_[1] - a_[0]);
 
 	for(std::size_t i{2}; i<nKnots_-2; ++i) {
 		b_[i] = TMath::Abs(a_[i+2] - a_[i+1]);
 		c_[i] = TMath::Abs(a_[i]   - a_[i-1]);
 	}
 
 	b_[nKnots_-2] = TMath::Abs(an            - a_[nKnots_-1]);
 	b_[nKnots_-1] = TMath::Abs(anp1          - an           );
 
 	c_[nKnots_-2] = TMath::Abs(a_[nKnots_-2] - a_[nKnots_-3]);
 	c_[nKnots_-1] = TMath::Abs(a_[nKnots_-1] - a_[nKnots_-2]);
 
 	// dy/dx calculated as the weighted average of a[i] and a[i+1]:
 	// dy/dx_i = ( | a_i+2 - a_i+1 | a_i  +  | a_i - a_i-1 | a_i+1 ) / ( | a_i+2 - a_i+1 |  +  | a_i - a_i-1 | )
 	// in the special case a_i-1 == a_i != a_i+1 == a_i+2 this function is undefined so dy/dx is then defined as (a_i + a_i+1) / 2
 	for(std::size_t i{0}; i<nKnots_-2; ++i) {
 		if(b_[i]==0 && c_[i]==0) {
 			dydx_[i] = ( a_[i] + a_[i+1] ) / 2.;
 		} else {
 			dydx_[i] = ( b_[i] * a_[i] + c_[i] * a_[i+1] ) / ( b_[i] + c_[i] );
 		}
 	}
 
 	if(b_[nKnots_-1]==0 && c_[nKnots_-1]==0) {
 		dydx_[nKnots_-1] = ( a_[nKnots_-1] + an ) / 2.;
 	} else {
 		dydx_[nKnots_-1] = ( b_[nKnots_-1] * a_[nKnots_-1] + c_[nKnots_-1] * an ) / ( b_[nKnots_-1] + c_[nKnots_-1] );
 	}
 }