diff --git a/src/Lau1DCubicSpline.cc b/src/Lau1DCubicSpline.cc
index dc5008b..5cf3bf8 100644
--- a/src/Lau1DCubicSpline.cc
+++ b/src/Lau1DCubicSpline.cc
@@ -1,687 +1,717 @@
 
 /*
 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 <fstream>
 #include <iostream>
 #include <algorithm>
 
 #include <nlohmann/json.hpp>
 
 #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,
 				   const SplineType type,
 				   const BoundaryType leftBound,
 				   const BoundaryType rightBound,
 				   const Double_t dydx0, const Double_t dydxn) :
 	nKnots_{xs.size()},
 	x_{xs},
 	y_{ys},
 	type_{type},
 	leftBound_{leftBound},
 	rightBound_{rightBound},
 	dydx0_{dydx0},
 	dydxn_{dydxn}
 {
 	init();
 }
 
 std::unique_ptr<Lau1DCubicSpline> Lau1DCubicSpline::readFromJson(const TString& filename, const TString& splineName)
 {
 	using nlohmann::json;
 
 	std::ifstream in(filename, std::ios_base::in);
 	if ( not in ) {
 		std::cerr << "ERROR in Lau1DCubicSpline::readFromJson : couldn't open file \"" << filename << "\"" << std::endl;
 		return std::unique_ptr<Lau1DCubicSpline>{};
 	}
 
 	json j;
 	in >> j;
 
 	if ( splineName != "" ) {
 		j = j[splineName.Data()];
 	}
 
 	// NB can't use std::make_unique here because the default constructor is private
 	std::unique_ptr<Lau1DCubicSpline> spline{ new Lau1DCubicSpline };
 	j.get_to(*spline);
 	spline->init();
 
 	return spline;
 }
 
 void Lau1DCubicSpline::writeToJson(const TString& filename, const TString& splineName, const bool append) const
 {
 	using nlohmann::json;
 
 	// if we're appending, we need to first load in the exising JSON from the file
 	json j;
 	if ( append ) {
 		// if appending we must have a name, so check that first
 		if ( splineName == "" ) {
 			std::cerr << "ERROR in Lau1DCubicSpline::writeToJson : when appending to a file it is mandatory to give a name to the spline" << std::endl;
 			return;
 		}
 
 		std::ifstream in{filename, std::ios_base::in};
 		if ( not in ) {
 			std::cerr << "ERROR in Lau1DCubicSpline::writeToJson : couldn't open file \"" << filename << "\" for initial reading" << std::endl;
 			return;
 		}
 
 		in >> j;
 	}
 
 	if(splineName == "") {
 		j = *this;
 	} else {
 		j[splineName.Data()] = *this;
 	}
 
 	// now we can overwrite the file
 	std::ofstream out{filename, std::ios_base::out};
 	if(not out) {
 		std::cerr << "ERROR in Lau1DCubicSpline::writeToJson : couldn't open file \"" << filename << "\" for writing. No file will be written!" << std::endl;
 		return;
 	}
 
 	out << j.dump(4);
 	out << std::endl;
 	out.close();
 
 	return;
 }
 
 Double_t Lau1DCubicSpline::evaluate(const 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
 	std::size_t cell{0};
 
 	while( x > x_[cell+1] ) {
 		++cell;
 	}
 
 	// obtain x- and y-values of the neighbouring knots
 	const Double_t xLow  { x_[cell] };
 	const Double_t xHigh { x_[cell+1] };
 	const Double_t yLow  { y_[cell] };
 	const Double_t yHigh { y_[cell+1] };
 
 	if(type_ == SplineType::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
 	const Double_t t { (x - xLow) / (xHigh - xLow) };
 	const Double_t a {     dydx_[cell]   * (xHigh - xLow) - (yHigh - yLow) };
 	const Double_t b { -1.*dydx_[cell+1] * (xHigh - xLow) + (yHigh - yLow) };
 
 	const 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(const SplineType type)
 {
 	if(type_ != type) {
 		type_ = type;
 		this->calcDerivatives();
 	}
 }
 
 void Lau1DCubicSpline::updateBoundaryConditions(const BoundaryType leftBound,
 		const BoundaryType rightBound,
 		const Double_t dydx0, const Double_t dydxn)
 {
 	bool updateDerivatives{false};
 
 	if(leftBound_ != leftBound || rightBound_ != rightBound) {
 		leftBound_ = leftBound;
 		rightBound_ = rightBound;
 		updateDerivatives = true;
 	}
 
 	if(dydx0_ != dydx0) {
 		dydx0_ = dydx0;
 		if(leftBound_ == BoundaryType::Clamped) updateDerivatives = true;
 	}
 
 	if(dydxn_ != dydxn) {
 		dydxn_ = dydxn;
 		if(rightBound_ == BoundaryType::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 SplineType::StandardSpline:
 		case SplineType::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 SplineType::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 SplineType::LinearInterpolation:
 			{
 				result[0] = xH*yL-xL*yH;
 				result[1] = yH-yL;
 
 				for(auto& res : result){res /= h;}
 				break;
 			}
 	}
 
 	if (normalise) {
 		const 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
 {
 	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 for normalisation)
 		return this->evaluate( x[0] ) * par[0];
 	};
 
 	//Make the function
 	TF1* func = new TF1("FittedFunc",fitf,x_.front(),x_.back(),1);
 
 	func -> SetParNames("Normalisation");
 	if(normalise)
 	{
 		func->SetParameter(0, 1./func->Integral(x_.front(),x_.back()));
 	}
 	else{func->SetParameter(0,1.);}
 
 	return func;
 }
 
 TF1* Lau1DCubicSpline::normaliseToTH1(TH1* hist, const 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, const LauOutputLevel printLevel, const bool doWL)
 {
 	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
 	// - WL= do a weighted likelihood fit (not chi2)
 	TString fitOptions {"N 0 S"};
 	if(doWL){fitOptions += " 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 SplineType::StandardSpline :
 			this->calcDerivativesStandard();
 			break;
 		case SplineType::AkimaSpline :
 			this->calcDerivativesAkima();
 			break;
 		case SplineType::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
 	switch ( leftBound_ ) {
 		case BoundaryType::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]));
 				break;
 			}
 
 		case BoundaryType::NotAKnot :
 			{
 				// define the width, h, and the 'slope', delta, of the first cell
 				const Double_t h1{x_[1]-x_[0]};
 				const Double_t h2{x_[2]-x_[0]};
 				const Double_t delta1{(y_[1]-y_[0])/h1};
 				const Double_t 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);
 				break;
 			}
 
 		case BoundaryType::Clamped :
 			{
 				b_[0] = 1.;
 				c_[0] = 0.;
 				d_[0] = dydx0_;
 				break;
 			}
 	}
 
 	// set right boundary condition
 	switch ( rightBound_ ) {
 		case BoundaryType::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]));
 				break;
 			}
 		case BoundaryType::NotAKnot :
 			{
 				// define the width, h, and the 'slope', delta, of the last cell
 				const Double_t hnm1{x_[nKnots_-1]-x_[nKnots_-2]};
 				const Double_t hnm2(x_[nKnots_-2]-x_[nKnots_-3]);
 				const Double_t deltanm1{(y_[nKnots_-1]-y_[nKnots_-2])/hnm1};
 				const Double_t 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);
 				break;
 			}
 
 		case BoundaryType::Clamped :
 			{
 				a_[nKnots_-1] = 0.;
 				b_[nKnots_-1] = 1.;
 				d_[nKnots_-1] = dydxn_;
 				break;
 			}
 	}
 
 	// 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.0}, an{0.0}, anp1{0.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] );
 	}
 }
 
 Double_t Lau1DCubicSpline::findMaximum() const
 {
-	//as a starting point, find the highest knot
-	Double_t max = *std::max_element(y_.begin(), y_.end());
-	//if linear interpolation is used, the maximum is the highest knot
-	if(type_ == SplineType::LinearInterpolation){return max;}
+	// as a starting point, find the highest knot
+	Double_t max { *std::max_element(y_.begin(), y_.end()) };
 
-	//Look for any local maxima greater than max
-	for(size_t nKnot = 0; nKnot < nKnots_ - 1; ++nKnot)
+	// if linear interpolation is used, the maximum is the highest knot
+	if ( type_ == SplineType::LinearInterpolation ) {
+		return max;
+	}
+
+	// look for any local maxima greater than max
+	for( std::size_t iSegment {0}; iSegment < (nKnots_ - 1); ++iSegment )
 	{
-		std::array<Double_t,4> coeffs = this->getCoefficients(nKnot);
-		//first derivative is a parabola
-		Double_t a = coeffs[3]*3, b = coeffs[2]*2, c = coeffs[1];
-		Double_t descriminant = b*b - 4*a*c;
-		if(descriminant < 0){continue;}
-		Double_t x1,x2;
-		x1 = (-b + TMath::Sqrt(descriminant))/(2*a);
-		if(x1 > x_[nKnot] and x1 < x_[nKnot+1] and this->evaluate(x1) > max)
-		{max = this->evaluate(x1);}
-		x2 = (-b - TMath::Sqrt(descriminant))/(2*a);
-		if(x2 > x_[nKnot] and x2 < x_[nKnot+1] and this->evaluate(x2) > max)
-		{max = this->evaluate(x2);}
+		const std::array<Double_t,4> coeffs { this->getCoefficients(iSegment, false) };
+
+		// first derivative is a parabola
+		const Double_t a { coeffs[3]*3.0 };
+		const Double_t b { coeffs[2]*2.0 };
+		const Double_t c { coeffs[1] };
+		const Double_t discriminant { b*b - 4.0*a*c };
+
+		// if no real roots, skip this segment
+		if ( discriminant < 0.0) {
+			continue;
+		}
+
+		const Double_t x1 { (-b + TMath::Sqrt(discriminant))/(2.0*a) };
+		if ( x1 > x_[iSegment] and x1 < x_[iSegment+1] ) {
+			// this root is within the segment - check the 2nd derivative
+			const Double_t secondDeriv { 2.0 * a * x1 + b };
+			if ( secondDeriv < 0.0 ) {
+				// evaluate the function and update max if appropriate
+				const Double_t y1 { this->evaluate(x1) };
+				if ( y1 > max ) {
+					max = y1;
+				}
+			}
+		}
+
+		const Double_t x2 { (-b - TMath::Sqrt(discriminant))/(2.0*a) };
+		if ( x2 > x_[iSegment] and x2 < x_[iSegment+1] ) {
+			// this root is within the segment - check the 2nd derivative
+			const Double_t secondDeriv { 2.0 * a * x2 + b };
+			if ( secondDeriv < 0.0 ) {
+				// evaluate the function and update max if appropriate
+				const Double_t y2 { this->evaluate(x2) };
+				if ( y2 > max ) {
+					max = y2;
+				}
+			}
+		}
 	}
+
 	return max;
 }