diff --git a/src/Lau1DCubicSpline.cc b/src/Lau1DCubicSpline.cc
index c9e5cb3..a54ea22 100644
--- a/src/Lau1DCubicSpline.cc
+++ b/src/Lau1DCubicSpline.cc
@@ -1,721 +1,725 @@
 
 /*
 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 "LauJsonTools.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)
 {
 	// NB deliberately not using uniform initialisation here because of this issue:
 	// https://json.nlohmann.me/home/faq/#brace-initialization-yields-arrays
 	const nlohmann::json j = LauJsonTools::readJsonFile( fileName.Data(), splineName.Data(), LauJsonTools::JsonType::Object );
 	if ( j.is_null() ) {
 		if ( splineName != "" ) {
 			std::cerr << "ERROR in Lau1DCubicSpline::readFromJson : unable to retrieve JSON object from element \"" << splineName << "\" in file \"" << fileName << "\"" << std::endl;
 		} else {
 			std::cerr << "ERROR in Lau1DCubicSpline::readFromJson : unable to retrieve JSON object from root element of file \"" << fileName << "\"" << std::endl;
 		}
 		return {};
 	}
 
 	// 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 int indent) const
 {
 	const nlohmann::json j = *this;
 
 	const bool writeOK { LauJsonTools::writeJsonFile( j, fileName.Data(), splineName.Data(), append, indent ) };
 	if ( ! writeOK ) {
 		std::cerr << "ERROR in Lau1DCubicSpline::writeToJson : could not successfully write to file \"" << fileName << std::endl;
 	}
 }
 
 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::ExtraVerbose :
 		case LauOutputLevel::Verbose :
 			fitOptions += " V";
 			break;
 		case LauOutputLevel::Quiet :
+		case LauOutputLevel::None :
 			fitOptions += " Q";
 			break;
 		case LauOutputLevel::Standard :
 			break;
 	}
 
 	hist->Fit( "FittedFunc", fitOptions );
 
 	return FittedFunc;
 }
 
 TFitResultPtr Lau1DCubicSpline::fitToTH1(TH1* hist, const LauOutputLevel printLevel, const bool doWL, std::map<Int_t,Double_t> fixedParams)
 {
 	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]);
 	}
 	for(auto& [knot, val] : fixedParams)
 	{
 		FittedFunc.FixParameter(knot, val);
 	}
 
 	// 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::ExtraVerbose :
 		case LauOutputLevel::Verbose :
 			fitOptions += " V";
 			break;
 		case LauOutputLevel::Quiet :
+		case LauOutputLevel::None :
 			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;
 	}
 
 	// look for any local maxima greater than max
 	for( std::size_t iSegment {0}; iSegment < (nKnots_ - 1); ++iSegment )
 	{
 		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;
 }
 
 //! \cond DOXYGEN_IGNORE
 std::ostream& operator<<(std::ostream& out, const Lau1DCubicSpline::SplineType type)
 {
 	switch (type) {
 		case Lau1DCubicSpline::SplineType::StandardSpline :
 			out << "StandardSpline";
 			break;
 		case Lau1DCubicSpline::SplineType::AkimaSpline :
 			out << "AkimaSpline";
 			break;
 		case Lau1DCubicSpline::SplineType::LinearInterpolation :
 			out << "LinearInterpolation";
 			break;
 	}
 	return out;
 }
 
 std::ostream& operator<<(std::ostream& out, const Lau1DCubicSpline::BoundaryType type)
 {
 	switch (type) {
 		case Lau1DCubicSpline::BoundaryType::Clamped:
 			out << "Clamped";
 			break;
 		case Lau1DCubicSpline::BoundaryType::Natural :
 			out << "Natural";
 			break;
 		case Lau1DCubicSpline::BoundaryType::NotAKnot :
 			out << "NotAKnot";
 			break;
 	}
 	return out;
 }
 //! \endcond
diff --git a/src/LauSPlot.cc b/src/LauSPlot.cc
index 81c61df..8b4e4cb 100644
--- a/src/LauSPlot.cc
+++ b/src/LauSPlot.cc
@@ -1,1295 +1,1288 @@
 
 /*
 Copyright 2006 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 LauSPlot.cc
     \brief File containing implementation of LauSPlot class.
 
     Class for defining the SPlot technique based on TSplot from ROOT by the following authors:
     Muriel Pivk, Anna Kreshuk (10/2005).
     (Original copyright notice below)
 */
 
 /**********************************************************************
  *                                                                    *
  * Copyright (c) 2005 ROOT Foundation,  CERN/PH-SFT                   *
  *                                                                    *
  **********************************************************************/
 
 #include <cfloat>
 #include <iostream>
 #include <vector>
 
 #include "TEventList.h"
 #include "TFile.h"
 #include "TLeaf.h"
 #include "TMath.h"
 #include "TObjArray.h"
 #include "TSystem.h"
 #include "TTree.h"
 #include "TVirtualFitter.h"
 
 #include "LauSPlot.hh"
 
 extern void Yields(Int_t &, Double_t *, Double_t &f, Double_t *x, Int_t iflag);
 
 ClassImp(LauSPlot)
 
 
 LauSPlot::LauSPlot(const TString& fileName, const TString& treeName,
 		Int_t firstExpt,
 		Int_t nExpt,
 		const NameSet& variableNames,
 		const NumbMap& freeSpecies,
 		const NumbMap& fixdSpecies,
 		const TwoDMap& twodimPDFs,
 		Bool_t sigSplit,
 		Bool_t scfDPSmeared) :
 	fileName_(fileName),
 	inputTreeName_(treeName),
 	cnTreeName_(""),
 	sweightTreeName_(""),
 	file_(0),
 	inputTree_(0),
 	cnTree_(0),
 	sweightTree_(0),
 	eventList_(0),
 	variableNames_(variableNames),
 	freeSpecies_(freeSpecies),
 	fixdSpecies_(fixdSpecies),
 	origFreeSpecies_(freeSpecies),
 	origFixdSpecies_(fixdSpecies),
 	twodimPDFs_(twodimPDFs),
 	signalSplit_(sigSplit),
 	scfDPSmear_(scfDPSmeared),
 	readInput_(kFALSE),
 	definedCNBranches_(kFALSE),
 	definedSWeightBranches_(kFALSE),
 	firstExpt_(firstExpt),
 	nExpt_(nExpt),
 	iExpt_(0),
 	nEvents_(0),
 	nDiscVars_(variableNames.size()),
 	nFreeSpecies_(freeSpecies.size()),
 	nFixdSpecies_(fixdSpecies.size()),
 	nSpecies_(freeSpecies.size()+fixdSpecies.size())
 {
 	this->openInputFileAndTree();
 	this->readInputInfo();
 	this->createSWeightTree();
 	if (nFixdSpecies_>0) {
 		this->createCNTree();
 	}
 }
 
 LauSPlot::~LauSPlot()
 {
 	// seems that closing the file deletes the tree
 	// so only delete if the file is still open
 	if (file_ && file_->IsOpen()) {
 		delete inputTree_;  inputTree_ = 0;
 		delete sweightTree_; sweightTree_ = 0;
 		delete cnTree_; cnTree_ = 0;
 	}
 	delete file_; file_ = 0;
 }
 
 void LauSPlot::openInputFileAndTree()
 {
 	// first check whether we've already opened up the file or not
 	if (!file_) {
 		// if not, first check the filename and if all ok create the file
 		if (fileName_ == "") {
 			std::cerr<<"ERROR in LauSPlot::createFileAndTree : Bad filename supplied, not creating file or tree."<<std::endl;
 			return;
 		}
 		file_ = TFile::Open(fileName_, "update");
 		if (!file_ || file_->IsZombie() || !file_->IsWritable()) {
 			std::cerr<<"ERROR in LauSPlot::createFileAndTree : Problem opening file \""<<fileName_<<"\" for updating, can't do anything."<<std::endl;
 			return;
 		}
 	}
 	// next open the input tree for reading
 	if (!inputTree_) {
 		file_->cd();
 		inputTree_ = dynamic_cast<TTree*>(file_->Get(inputTreeName_));
 		inputTree_->SetDirectory(file_);
 	}
 }
 
 void LauSPlot::readInputInfo()
 {
 	// Read the tree and then setup the maps so we know which leaves to
 	// read from the tree to get the various PDF values
 
 	Bool_t inputOK = this->readInputLeaves();
 	if (!inputOK) {
 		this->readInput(inputOK);
 		return;
 	}
 
 	inputOK = this->checkLeaves();
 	this->readInput(inputOK);
 	return;
 }
 
 Bool_t LauSPlot::readInputLeaves()
 {
 	// Read all the leaves in the tree and store them in the leaves map
 
 	if (!inputTree_) {
 		std::cerr<<"ERROR in LauSPlot::readInputInfo : Invalid pointer to data tree."<<std::endl;
 		return kFALSE;
 	}
 
 	Int_t nBranches = inputTree_ ? static_cast<Int_t>(inputTree_->GetNbranches()) : 0;
 	TObjArray* pLeaves = inputTree_->GetListOfLeaves();
 	if (!pLeaves) {
 		std::cerr<<"ERROR in LauSPlot::readInputInfo : Problem retrieving leaves from the tree."<<std::endl;
 		return kFALSE;
 	}
 	TObjArray& leaves = *pLeaves;
 
 	if (nBranches > leaves.GetSize()) {
 		std::cerr<<"ERROR in LauSPlot::readInputInfo : List of leaves is smaller than number of branches - this is strange!"<<std::endl;
 		return kFALSE;
 	}
 
 	for (Int_t iLeaf = 0; iLeaf < nBranches; ++iLeaf) {
 
 		TLeaf * leaf = dynamic_cast<TLeaf*>(leaves[iLeaf]);
 
 		// we can't deal with arrays
 		Int_t size = leaf->GetNdata();
 		if (size != 1) {
 			std::cerr<<"ERROR in LauSPlot::readInputInfo : Tree has array branches, can't deal with those."<<std::endl;
 			return kFALSE;
 		}
 
 		// find the name of the leaf
 		TString name = leaf->GetName();
 
 		// initialise an entry in the maps to hold the value
 		leaves_[name] = leaf;
 	}
 
 	return kTRUE;
 }
 
 Bool_t LauSPlot::checkLeaves() const
 {
 	// Analyse the names of the leaves to check that we have a leaf for
 	// all bits of information we require, i.e. a likelihood value for
 	// each combination of variable and species
 
 	// If we have 2D PDFs then we have to look for some extra leaves
 	for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
 
 		const TString& specName = twodim_iter->first;
 		const TString& firstVarName = twodim_iter->second.first;
 		const TString& secondVarName = twodim_iter->second.second;
 
 		TString expectedName(specName);
 		expectedName += firstVarName;
 		expectedName += secondVarName;
 		expectedName += "Like";
 
 		Bool_t found(kFALSE);
 		for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 			const TString& leafName = leaf_iter->first;
 			const TLeaf* leaf = leaf_iter->second;
 
 			if ((leafName == expectedName) && (leaf != 0)) {
 				found = kTRUE;
 				break;
 			}
 		}
 
 		if (!found) {
 			std::cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<std::endl;
 			return kFALSE;
 		}
 	}
 
 	// Search for all other "standard" leaves, i.e. <species><var>Like
 
 	for (NumbMap::const_iterator fixd_iter = fixdSpecies_.begin(); fixd_iter != fixdSpecies_.end(); ++fixd_iter) {
 
 		const TString& specName = fixd_iter->first;
 
 		// if the signal is split we need to do a dedicated search
 		// for sigTM and sigSCF, so skip the signal here
 		if ( specName == "sig" && this->signalSplit() ) {
 			continue;
 		}
 
 		for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
 			const TString& varName = (*vars_iter);
 			TString expectedName(specName);
 			expectedName += varName;
 			expectedName += "Like";
 			Bool_t found(kFALSE);
 			for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 				const TString& leafName = leaf_iter->first;
 				const TLeaf* leaf = leaf_iter->second;
 
 				if ((leafName == expectedName) && (leaf != 0)) {
 					found = kTRUE;
 					break;
 				}
 			}
 
 			if (!found) {
 				std::cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<std::endl;
 				return kFALSE;
 			}
 		}
 	}
 
 	for (NumbMap::const_iterator free_iter = freeSpecies_.begin(); free_iter != freeSpecies_.end(); ++free_iter) {
 
 		const TString& specName = free_iter->first;
 
 		// if the signal is split we need to do a dedicated search
 		// for sigTM and sigSCF, so skip the signal here
 		if ( specName == "sig" && this->signalSplit() ) {
 			continue;
 		}
 
 		for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
 			const TString& varName = (*vars_iter);
 			TString expectedName(specName);
 			expectedName += varName;
 			expectedName += "Like";
 			Bool_t found(kFALSE);
 			for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 				const TString& leafName = leaf_iter->first;
 				const TLeaf* leaf = leaf_iter->second;
 
 				if ((leafName == expectedName) && (leaf != 0)) {
 					found = kTRUE;
 					break;
 				}
 			}
 
 			if (!found) {
 				std::cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<std::endl;
 				return kFALSE;
 			}
 		}
 	}
 
 	if ( this->signalSplit() ) {
 
 		// now need to search for the sigTM and sigSCF leaves
 
 		for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
 			const TString& varName = (*vars_iter);
 			TString expectedName("sigTM");
 			expectedName += varName;
 			expectedName += "Like";
 			Bool_t found(kFALSE);
 			for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 				const TString& leafName = leaf_iter->first;
 				const TLeaf* leaf = leaf_iter->second;
 
 				if ((leafName == expectedName) && (leaf != 0)) {
 					found = kTRUE;
 					break;
 				}
 			}
 
 			if (!found) {
 				std::cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<std::endl;
 				return kFALSE;
 			}
 
 			expectedName = "sigSCF";
 			expectedName += varName;
 			expectedName += "Like";
 			found = kFALSE;
 			for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 				const TString& leafName = leaf_iter->first;
 				const TLeaf* leaf = leaf_iter->second;
 
 				if ((leafName == expectedName) && (leaf != 0)) {
 					found = kTRUE;
 					break;
 				}
 			}
 
 			if (!found) {
 				std::cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<std::endl;
 				return kFALSE;
 			}
 
 			expectedName = "sigSCFFrac";
 			found = kFALSE;
 			for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 				const TString& leafName = leaf_iter->first;
 				const TLeaf* leaf = leaf_iter->second;
 
 				if ((leafName == expectedName) && (leaf != 0)) {
 					found = kTRUE;
 					break;
 				}
 			}
 
 			if (!found) {
 				std::cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<std::endl;
 				return kFALSE;
 			}
 		}
 	}
 
 	return kTRUE;
 }
 
 void LauSPlot::createCNTree()
 {
 	// check whether we've already created the tree
 	if (!cnTree_) {
 		// if not change to the file's directory and create the tree
 		cnTreeName_ = inputTreeName_;
 		cnTreeName_ += "_cNCoeffs";
 		file_->cd();
 		cnTree_ = new TTree(cnTreeName_, cnTreeName_);
 		cnTree_->SetDirectory(file_);
 		this->definedCNBranches(kFALSE);
 	}
 }
 
 void LauSPlot::createSWeightTree()
 {
 	// check whether we've already created the tree
 	if (!sweightTree_) {
 		// if not change to the file's directory and create the tree
 		sweightTreeName_ = inputTreeName_;
 		sweightTreeName_ += "_sWeights";
 		file_->cd();
 		sweightTree_ = new TTree(sweightTreeName_, sweightTreeName_);
 		sweightTree_->SetDirectory(file_);
 		this->definedSWeightBranches(kFALSE);
 	}
 }
 
 void LauSPlot::defineCNBranches()
 {
 	if (this->definedCNBranches()) {
 		std::cerr<<"ERROR in LauSPlot::defineCNBranches : Already defined branches, not doing it again."<<std::endl;
 		return;
 	}
 	if (cN_.empty()) {
 		std::cerr<<"ERROR in LauSPlot::defineCNBranches : No entries in the cN container, can't define branches."<<std::endl;
 		return;
 	}
 	if (!cnTree_) {
 		std::cerr<<"ERROR in LauSPlot::defineCNBranches : Invalid pointer to the tree, can't define branches."<<std::endl;
 		return;
 	}
 
 	// In the cN tree there is one entry per experiment, so need to know the experiment number.
 	cnTree_->Branch("iExpt", &iExpt_, "iExpt/I");
 
 	for (std::map<TString,NumbMap>::iterator excl_iter = cN_.begin(); excl_iter != cN_.end(); ++excl_iter) {
 		const TString& exclName = excl_iter->first;
 		NumbMap& species = excl_iter->second;
 		for (NumbMap::iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
 			const TString& specName = spec_iter->first;
 			Double_t * pointer = &(spec_iter->second);
 			TString name(specName); name += "_cN";
 			if (exclName == "none") {
 				name += "_all";
 			} else {
 				name += "_no";
 				name += exclName;
 			}
 			TString thirdPart(name);  thirdPart += "/D";
 			cnTree_->Branch(name, pointer, thirdPart);
 		}
 	}
 	this->definedCNBranches(kTRUE);
 }
 
 void LauSPlot::defineSWeightBranches()
 {
 	if (this->definedSWeightBranches()) {
 		std::cerr<<"ERROR in LauSPlot::defineSWeightBranches : Already defined branches, not doing it again."<<std::endl;
 		return;
 	}
 	if (sWeights_.empty()) {
 		std::cerr<<"ERROR in LauSPlot::defineSWeightBranches : No entries in the sWeights container, can't define branches."<<std::endl;
 		return;
 	}
 	if (!sweightTree_) {
 		std::cerr<<"ERROR in LauSPlot::defineSWeightBranches : Invalid pointer to the tree, can't define branches."<<std::endl;
 		return;
 	}
 
 	for (std::map<TString,NumbMap>::iterator excl_iter = sWeightsCurrent_.begin(); excl_iter != sWeightsCurrent_.end(); ++excl_iter) {
 		const TString& exclName = excl_iter->first;
 		NumbMap& species = excl_iter->second;
 		for (NumbMap::iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
 			const TString& specName = spec_iter->first;
 			Double_t * pointer = &(spec_iter->second);
 			TString name(specName); name += "_sWeight";
 			if (exclName == "none") {
 				name += "_all";
 			} else {
 				name += "_no";
 				name += exclName;
 			}
 			TString thirdPart(name);  thirdPart += "/D";
 			sweightTree_->Branch(name, pointer, thirdPart);
 		}
 	}
 	this->definedSWeightBranches(kTRUE);
 }
 
 void LauSPlot::setExperiment(Int_t iExpt)
 {
 	if (!inputTree_) {
 		std::cerr<<"ERROR in LauSPlot::setExperiment : Invalid pointer to data tree."<<std::endl;
 		return;
 	}
 
 	// create the event list if we haven't already done so
 	if (!eventList_) {
 		eventList_ = new TEventList("splotelist","splotelist");
 		eventList_->SetDirectory(file_);
 	}
 
 	// fill the event list with this experiment's events
 	TString listName(eventList_->GetName());
 	listName.Prepend(">>");
 	TString selection("iExpt==");
 	selection += iExpt;
 	std::cout<<"LauSPlot::setExperiment : Setting tree to experiment number "<<iExpt<<"."<<std::endl;
 	inputTree_->Draw(listName,selection);
 
 	// find how many events there are
 	nEvents_ = eventList_->GetN();
 	std::cout<<"                          Found "<<nEvents_<<" events."<<std::endl;
 
 	// make sure we have enough space in the per-event value vectors
 	pdfTot_.clear(); pdfTot_.resize(nEvents_);
 	discPdf_.clear(); discPdf_.resize(nEvents_);
 	scfFrac_.clear(); scfFrac_.resize(nEvents_);
 	sWeights_.clear(); sWeights_.resize(nEvents_);
 
 	// read the info for this experiment
 	this->readExpt();
 }
 
 void LauSPlot::readExpt()
 {
 	for (Int_t iEvt(0); iEvt < nEvents_; ++iEvt) {
 		// Find which entry from the full tree contains the requested event
 		Long64_t iEntry = eventList_ ? eventList_->GetEntry(iEvt) : iEvt;
 		if (iEntry<0) { // this shouldn't happen, but just in case...
 			std::cerr<<"ERROR in LauSPlot::readExpt : Problem retrieving event."<<std::endl;
 			gSystem->Exit(EXIT_FAILURE);
 		}
 
 		// Then retrieve that entry from the tree
 		inputTree_->GetEntry(iEntry);
 
 		// If needed retrieve the SCF fraction values
 		if ( this->signalSplit() ) {
 			for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 				const TString& leafName = leaf_iter->first;
 				const TLeaf* leaf = leaf_iter->second;
 				if ((leafName == "sigSCFFrac") && (leaf != 0)) {
 					scfFrac_[iEvt] = leaf->GetValue();
 					break;
 				}
 			}
 		}
 
 		// Copy the leaf values into discPdf_
 		for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
 			const TString& specName = twodim_iter->first;
 			const TString& firstVarName = twodim_iter->second.first;
 			const TString& secondVarName = twodim_iter->second.second;
 
 			TString varName = firstVarName + secondVarName;
 
 			TString expectedName(specName);
 			expectedName += varName;
 			expectedName += "Like";
 
 			for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 				const TString& leafName = leaf_iter->first;
 				const TLeaf* leaf = leaf_iter->second;
 				if ((leafName == expectedName) && (leaf != 0)) {
 					discPdf_[iEvt][specName][varName] = leaf->GetValue();
 					break;
 				}
 			}
 		}
 
 		Bool_t needSignalSearch(kFALSE);
 		for (NumbMap::const_iterator fixd_iter = fixdSpecies_.begin(); fixd_iter != fixdSpecies_.end(); ++fixd_iter) {
 			const TString& specName = fixd_iter->first;
 
 			// if the signal is split we need to do a dedicated search
 			// for sigTM and sigSCF, so skip the signal here
 			if ( specName == "sig" && this->signalSplit() ) {
 				needSignalSearch = kTRUE;
 				continue;
 			}
 
 			for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
 				const TString& varName = (*vars_iter);
 				TString expectedName(specName);
 				expectedName += varName;
 				expectedName += "Like";
 				for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 					const TString& leafName = leaf_iter->first;
 					const TLeaf* leaf = leaf_iter->second;
 					if ((leafName == expectedName) && (leaf != 0)) {
 						discPdf_[iEvt][specName][varName] = leaf->GetValue();
 						break;
 					}
 				}
 			}
 		}
 
 		for (NumbMap::const_iterator free_iter = freeSpecies_.begin(); free_iter != freeSpecies_.end(); ++free_iter) {
 			const TString& specName = free_iter->first;
 
 			// if the signal is split we need to do a dedicated search
 			// for sigTM and sigSCF, so skip the signal here
 			if ( specName == "sig" && this->signalSplit() ) {
 				needSignalSearch = kTRUE;
 				continue;
 			}
 
 			for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
 				const TString& varName = (*vars_iter);
 				TString expectedName(specName);
 				expectedName += varName;
 				expectedName += "Like";
 				for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 					const TString& leafName = leaf_iter->first;
 					const TLeaf* leaf = leaf_iter->second;
 					if ((leafName == expectedName) && (leaf != 0)) {
 						discPdf_[iEvt][specName][varName] = leaf->GetValue();
 						break;
 					}
 				}
 			}
 		}
 
 		if ( needSignalSearch ) {
 			for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
 				const TString& varName = (*vars_iter);
 				TString specName("sigTM");
 				TString expectedName(specName);
 				expectedName += varName;
 				expectedName += "Like";
 				for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 					const TString& leafName = leaf_iter->first;
 					const TLeaf* leaf = leaf_iter->second;
 					if ((leafName == expectedName) && (leaf != 0)) {
 						discPdf_[iEvt][specName][varName] = leaf->GetValue();
 						break;
 					}
 				}
 				specName = "sigSCF";
 				expectedName = specName;
 				expectedName += varName;
 				expectedName += "Like";
 				for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
 					const TString& leafName = leaf_iter->first;
 					const TLeaf* leaf = leaf_iter->second;
 					if ((leafName == expectedName) && (leaf != 0)) {
 						discPdf_[iEvt][specName][varName] = leaf->GetValue();
 						break;
 					}
 				}
 			}
 		}
 	}
 }
 
 void LauSPlot::runCalculations(const LauOutputLevel printLevel)
 {
 	// Calculates the sWeights and cN coeffs
 
 	if (!this->readInput()) {
 		std::cerr<<"ERROR in LauSPlot::runCalculations : The input ntuple has not been successfully read, can't calculate anything."<<std::endl;
 		return;
 	}
 
 	if (freeSpecies_.empty()) {
 		std::cerr<<"ERROR in LauSPlot::runCalculations : Numbers of events in species have not been set."<<std::endl;
 		return;
 	}
 
 	// Make sure that global fitter is MINUIT
 	this->checkFitter();
 
 	// Loop over experiments
 	for (iExpt_ = firstExpt_; iExpt_ < (firstExpt_+nExpt_); ++iExpt_) {
 
 		this->setExperiment(iExpt_);
 
 		if (nEvents_ < 1) {
 			std::cerr<<"ERROR in LauSPlot::runCalculations : Zero events in experiment "<<iExpt_<<", skipping..."<<std::endl;
 			continue;
 		}
 
 		// Now loop over the PDFs to exclude, including the case where none are to be excluded.
 		NameSet excludePdf;
 		if (variableNames_.size()<2) {
 			excludePdf.insert("none");
 		} else {
 			excludePdf = variableNames_;
 			excludePdf.insert("none");
 		}
 
 		for (NameSet::const_iterator excl_iter = excludePdf.begin(); excl_iter != excludePdf.end(); ++excl_iter) {
 
 			const TString& exclName = (*excl_iter);
 
 			std::cout<<"LauSPlot::runCalculations : Calculating sWeights, excluding PDF: "<<exclName<<"."<<std::endl;
 
 			// Calculate the per-event total PDF values for each species.
 			this->calcTotPDFValues(exclName);
 
 			// Reset the fitter
 			this->initialiseFitter(printLevel);
 
 			// Set the parameters to their initial values
 			this->setFitParameters();
 
 			// Run the fit
 			this->runFit();
 
 			// Get the fitted parameter values back
 			this->retrieveFittedParameters(printLevel);
 
 			// Get the covariance matrix
 			Bool_t covMatOK = this->calcCovMatrix();
 			Double_t * covmat(0);
 			if (!covMatOK) {
 				TVirtualFitter * fitter = TVirtualFitter::GetFitter();
 				covmat = fitter->GetCovarianceMatrix();
 			}
 			if (printLevel == LauOutputLevel::Verbose) {
 				this->printCovMatrixElements(covmat);
 			}
 
 			// calculate the cN and sWeights from the covariance matrix
 			if (nFixdSpecies_ > 0) {
 				this->calcCNCoeffs(exclName, covmat);
 			}
 			this->calcSWeights(exclName, covmat);
 
 			// print verbose info if required
 			if (printLevel == LauOutputLevel::Verbose) {
 				this->printSumOfWeights(exclName);
 			}
 		}
 
 		// Finally fill all the branches for this experiment
 		if (nFixdSpecies_ > 0) {
 			this->fillCNBranches();
 		}
 		this->fillSWeightBranches();
 	}
 }
 
 void LauSPlot::checkFitter() const
 {
 	TString minuitName("TFitter");
 	TVirtualFitter * fitter = TVirtualFitter::GetFitter();
 	if (fitter) {
 		TString fitterName(fitter->IsA()->GetName());
 		if (fitterName != minuitName) {
 			delete fitter; fitter = 0;
 		}
 	}
 	fitter = TVirtualFitter::Fitter(0);
 }
 
 void LauSPlot::initialiseFitter(const LauOutputLevel printLevel)
 {
 	TVirtualFitter * fitter = TVirtualFitter::GetFitter();
 	fitter->Clear();
 	fitter->SetFCN(Yields);
 	fitter->SetObjectFit(this);
 
 	// Set the print level
-	Double_t arglist[10];
-	switch ( printLevel ) {
-		case LauOutputLevel::Quiet :
-			arglist[0]=-1;
-			break;
-		case LauOutputLevel::Standard :
-			arglist[0]=0;
-			break;
-		case LauOutputLevel::Verbose :
-			arglist[0]=1;
-			break;
+	std::array<Double_t,1> arglist;
+	arglist[0] = static_cast<Double_t>(printLevel);
+	fitter->ExecuteCommand("SET PRINT", arglist.data(), 1);
+	if ( printLevel == LauOutputLevel::None ) {
+		fitter->ExecuteCommand("SET NOWARNINGS", arglist.data(), 0);
 	}
-	fitter->ExecuteCommand("SET PRINT", arglist, 1);
 
 	// Need to set the "SET ERR" command to +0.5 for +/-1 sigma errors
 	// for maximum likelihood fit.
 	arglist[0] = 0.5;
-	fitter->ExecuteCommand("SET ERR", arglist, 1);
+	fitter->ExecuteCommand("SET ERR", arglist.data(), 1);
 }
 
 void LauSPlot::setFitParameters() const
 {
 	TVirtualFitter * fitter = TVirtualFitter::GetFitter();
 
 	// must add the parameters in the same order as they are stored in pdfTot_
 	Int_t ispecies(0);
 	const NumbMap& species = pdfTot_.front();
 
 	for (NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
 		const TString& name(spec_iter->first);
 
 		// starting parameters should be the original values,
 		// not those that came out of the last fit
 		NumbMap::const_iterator free_iter = origFreeSpecies_.find(name);
 		NumbMap::const_iterator fixd_iter = origFixdSpecies_.find(name);
 		Bool_t fixed = fixd_iter != origFixdSpecies_.end();
 
 		Double_t value(0.0);
 		if (fixed) {
 			value = fixd_iter->second;
 		} else {
 			value = free_iter->second;
 		}
 		fitter->SetParameter(ispecies, name, value, 1, 0, 0);
 		if (fixed) {
 			fitter->FixParameter(ispecies);
 		}
 		++ispecies;
 	}
 }
 
 void LauSPlot::runFit()
 {
 	TVirtualFitter * fitter = TVirtualFitter::GetFitter();
 
 	Double_t arglist[10];
 	arglist[0] = 1000*nFreeSpecies_; // maximum iterations
 	arglist[1] = 0.05;               // tolerance : min EDM = 0.001*tolerance
 
 	// Execute MIGRAD
 	Int_t fitStatus = fitter->ExecuteCommand("MIGRAD", arglist, 2);
 
 	if (fitStatus != 0) {
 		std::cerr<<"ERROR in LauSPlot::runFit : Error during MIGRAD minimisation."<<std::endl;
 	} else {
 		// Execute HESSE
 		fitStatus = fitter->ExecuteCommand("HESSE", arglist, 1);
 		if (fitStatus != 0) {
 			std::cerr<<"ERROR in LauSPlot::runFit : Error during HESSE error calculation."<<std::endl;
 		}
 	}
 }
 
 void LauSPlot::printCovMatrixElements(const Double_t * covmat) const
 {
 	Int_t a(0);
 	std::cout<<std::endl;
 	for (NumbMap::const_iterator iter_a = freeSpecies_.begin(); iter_a != freeSpecies_.end(); ++iter_a) {
 		Int_t b(0);
 		for (NumbMap::const_iterator iter_b = freeSpecies_.begin(); iter_b != freeSpecies_.end(); ++iter_b) {
 			if (covmat) {
 				std::cout<<"CovMat Elem: ["<<iter_a->first<<","<<iter_b->first<<"] = "<<covmat[a*nSpecies_+b]<<std::endl;
 			} else {
 				std::cout<<"CovMat Elem: ["<<iter_a->first<<","<<iter_b->first<<"] = "<<covMat_(a,b)<<std::endl;
 			}
 			++b;
 		}
 		++a;
 	}
 	std::cout<<std::endl;
 
 	if (!covmat) {
 		covMat_.Print();
 	}
 }
 
 void LauSPlot::retrieveFittedParameters(const LauOutputLevel printLevel)
 {
 	TVirtualFitter * fitter = TVirtualFitter::GetFitter();
 
 	// remember fit parameters are in same order as in pdfTot_
 	Int_t ispecies(0);
 	const NumbMap& species = pdfTot_.front();
 
 	for (NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
 		const TString& name(spec_iter->first);
 		NumbMap::iterator free_iter = freeSpecies_.find(name);
 		if (free_iter != freeSpecies_.end()) {
 			free_iter->second = fitter->GetParameter(ispecies);
 			if ( printLevel != LauOutputLevel::Quiet ) {
 				std::cout<<"Estimated # of events in species "<<name<<" = "<<free_iter->second<<std::endl;
 			}
 		}
 		++ispecies;
 	}
 	if ( printLevel != LauOutputLevel::Quiet ) {
 		std::cout<<std::endl;
 	}
 }
 
 void LauSPlot::printSumOfWeights(const TString& exclName) const
 {
 	for (NumbMap::const_iterator spec_iter = freeSpecies_.begin(); spec_iter != freeSpecies_.end(); ++spec_iter) {
 		const TString& specName = spec_iter->first;
 		Double_t sumweight(0.0);
 		for (Int_t iEvt(0); iEvt < nEvents_; ++iEvt) {
 			const NumbMap& specWeightMap = sWeights_[iEvt].find(exclName)->second;
 			Double_t weight = specWeightMap.find(specName)->second;
 			sumweight += weight;
 		}
 		std::cout<<"Sum of sWeights for species \""<<specName<<"\" = "<<sumweight<<std::endl;
 	}
 	std::cout<<std::endl;
 }
 
 Bool_t LauSPlot::calcCovMatrix()
 {
 	// Calculate our inverse covariance matrix from the various PDFs
 
 	TMatrixD invMatrix(nFreeSpecies_,nFreeSpecies_);
 
 	// First calculate the denominator, which is common to all elements
 	std::vector<Double_t> denominator(nEvents_);
 	for (Int_t iEvt(0); iEvt<nEvents_; ++iEvt) {
 		denominator[iEvt] = 0.0;
 		for (NumbMap::const_iterator spec_iter = freeSpecies_.begin(); spec_iter != freeSpecies_.end(); ++spec_iter) {
 			const TString& specName = spec_iter->first;
 			Double_t specNumEvents = spec_iter->second;
 			denominator[iEvt] += specNumEvents * pdfTot_[iEvt][specName];
 		}
 		for (NumbMap::const_iterator spec_iter = fixdSpecies_.begin(); spec_iter != fixdSpecies_.end(); ++spec_iter) {
 			const TString& specName = spec_iter->first;
 			Double_t specNumEvents = spec_iter->second;
 			denominator[iEvt] += specNumEvents * pdfTot_[iEvt][specName];
 		}
 		// Square to get the final denominator
 		denominator[iEvt] *= denominator[iEvt];
 	}
 
 	// Then calculate each element
 	Int_t i(0);
 	for (NumbMap::const_iterator spec_iter_i = freeSpecies_.begin(); spec_iter_i != freeSpecies_.end(); ++spec_iter_i) {
 		const TString& specName_i = spec_iter_i->first;
 		Int_t j(0);
 		for (NumbMap::const_iterator spec_iter_j = freeSpecies_.begin(); spec_iter_j != freeSpecies_.end(); ++spec_iter_j) {
 			const TString& specName_j = spec_iter_j->first;
 			invMatrix(i,j) = 0.0;
 			for (Int_t iEvt(0); iEvt<nEvents_; ++iEvt) {
 				Double_t numerator = pdfTot_[iEvt][specName_i] * pdfTot_[iEvt][specName_j];
 				invMatrix(i,j) += numerator/denominator[iEvt];
 			}
 			++j;
 		}
 		++i;
 	}
 
 	// Check for a singular matrix
 	if (invMatrix.Determinant() < 1e-15) {
 		std::cerr<<"ERROR in LauSPlot::calcCovMatrix : Calculated inverse covariance matrix is singular, can't invert it."<<std::endl;
 		return kFALSE;
 	}
 
 	// Invert and store in the covariance matrix
 	covMat_.ResizeTo(nFreeSpecies_,nFreeSpecies_);
 	covMat_ = TMatrixD(TMatrixD::kInverted, invMatrix);
 
 	return kTRUE;
 }
 
 void LauSPlot::calcTotPDFValues(const TString& exclName)
 {
 	for (Int_t iEvt(0); iEvt<nEvents_; iEvt++) {
 
 		Bool_t needSignalSearch(kFALSE);
 
 		for (NumbMap::const_iterator spec_iter = fixdSpecies_.begin(); spec_iter != fixdSpecies_.end(); ++spec_iter) {
 			const TString& specName = spec_iter->first;
 
 			// if the signal is split we need to do a dedicated search
 			// for sigTM and sigSCF, so skip the signal here
 			if ( specName == "sig" && this->signalSplit() ) {
 				needSignalSearch = kTRUE;
 				continue;
 			}
 
 			pdfTot_[iEvt][specName] = 1.0;
 
 			// loop through the 2D histo list
 			NameSet skipList;
 			for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
 				// if the entry doesn't refer to this
 				// species then skip on
 				if ( specName != twodim_iter->first ) {
 					continue;
 				}
 
 				// retrieve the two variable names
 				const TString& firstVarName = twodim_iter->second.first;
 				const TString& secondVarName = twodim_iter->second.second;
 				if ( firstVarName != exclName && secondVarName != exclName ) {
 					// if neither name is the one being excluded then...
 					// add them both to the skip list
 					skipList.insert( firstVarName );
 					skipList.insert( secondVarName );
 					// and multiply the total by the combined PDF value
 					TString varName = firstVarName + secondVarName;
 					pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
 				}
 			}
 
 			// loop through all the variables
 			for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
 				const TString& varName = (*var_iter);
 				// if the variable isn't the one being excluded
 				if (exclName != varName) {
 					// and it's not involved in a 2D PDF
 					if ( skipList.find(varName) == skipList.end() ) {
 						// multiply the total by its PDF value
 						pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
 					}
 				}
 			}
 		}
 
 		for (NumbMap::const_iterator spec_iter = freeSpecies_.begin(); spec_iter != freeSpecies_.end(); ++spec_iter) {
 			const TString& specName = spec_iter->first;
 
 			// if the signal is split we need to do a dedicated search
 			// for sigTM and sigSCF, so skip the signal here
 			if ( specName == "sig" && this->signalSplit() ) {
 				needSignalSearch = kTRUE;
 				continue;
 			}
 
 			pdfTot_[iEvt][specName] = 1.0;
 
 			// loop through the 2D histo list
 			NameSet skipList;
 			for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
 				// if the entry doesn't refer to this
 				// species then skip on
 				if ( specName != twodim_iter->first ) {
 					continue;
 				}
 
 				// retrieve the two variable names
 				const TString& firstVarName = twodim_iter->second.first;
 				const TString& secondVarName = twodim_iter->second.second;
 				if ( firstVarName != exclName && secondVarName != exclName ) {
 					// if neither name is the one being excluded then...
 					// add them both to the skip list
 					skipList.insert( firstVarName );
 					skipList.insert( secondVarName );
 					// and multiply the total by the combined PDF value
 					TString varName = firstVarName + secondVarName;
 					pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
 				}
 			}
 
 			// loop through all the variables
 			for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
 				const TString& varName = (*var_iter);
 				// if the variable isn't the one being excluded
 				if (exclName != varName) {
 					// and it's not involved in a 2D PDF
 					if ( skipList.find(varName) == skipList.end() ) {
 						// multiply the total by its PDF value
 						pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
 					}
 				}
 			}
 		}
 
 		if ( needSignalSearch ) {
 
 			// loop through the 2D histo list
 			TString specName("sigTM");
 			Double_t tmPDFVal(1.0);
 			NameSet skipList;
 			for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
 				// if the entry doesn't refer to this
 				// species then skip on
 				if ( specName != twodim_iter->first ) {
 					continue;
 				}
 
 				// retrieve the two variable names
 				const TString& firstVarName = twodim_iter->second.first;
 				const TString& secondVarName = twodim_iter->second.second;
 				if ( firstVarName != exclName && secondVarName != exclName ) {
 					// if neither name is the one being excluded then...
 					// add them both to the skip list
 					skipList.insert( firstVarName );
 					skipList.insert( secondVarName );
 					// and multiply the total by the combined PDF value
 					TString varName = firstVarName + secondVarName;
 					tmPDFVal *= discPdf_[iEvt][specName][varName];
 				}
 			}
 
 			// loop through all the variables
 			for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
 				const TString& varName = (*var_iter);
 				// if the variable isn't the one being excluded
 				if (exclName != varName) {
 					// and it's not involved in a 2D PDF
 					if ( skipList.find(varName) == skipList.end() ) {
 						// multiply the total by its PDF value
 						tmPDFVal *= discPdf_[iEvt][specName][varName];
 					}
 				}
 			}
 
 			tmPDFVal *= (1.0 - scfFrac_[iEvt]);
 
 
 			// loop through the 2D histo list
 			specName = "sigSCF";
 			Double_t scfPDFVal(1.0);
 			skipList.clear();
 			for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
 				// if the entry doesn't refer to this
 				// species then skip on
 				if ( specName != twodim_iter->first ) {
 					continue;
 				}
 
 				// retrieve the two variable names
 				const TString& firstVarName = twodim_iter->second.first;
 				const TString& secondVarName = twodim_iter->second.second;
 				if ( firstVarName != exclName && secondVarName != exclName ) {
 					// if neither name is the one being excluded then...
 					// add them both to the skip list
 					skipList.insert( firstVarName );
 					skipList.insert( secondVarName );
 					// and multiply the total by the combined PDF value
 					TString varName = firstVarName + secondVarName;
 					scfPDFVal *= discPdf_[iEvt][specName][varName];
 				}
 			}
 
 			// loop through all the variables
 			for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
 				const TString& varName = (*var_iter);
 				// if the variable isn't the one being excluded
 				if (exclName != varName) {
 					// and it's not involved in a 2D PDF
 					if ( skipList.find(varName) == skipList.end() ) {
 						// multiply the total by its PDF value
 						scfPDFVal *= discPdf_[iEvt][specName][varName];
 					}
 				}
 			}
 
 			if ( exclName == "DP" || !this->scfDPSmear() ) {
 				scfPDFVal *= scfFrac_[iEvt];
 			}
 
 			pdfTot_[iEvt]["sig"] = tmPDFVal + scfPDFVal;
 		}
 	}
 }
 
 void LauSPlot::calcCNCoeffs(const TString& exclName, const Double_t *covmat)
 {
 	// Computes the cN for the extended sPlots from the covariance matrix
 
 	if (nFixdSpecies_ <= 0) {
 		return;
 	}
 
 	Int_t species_n(0);
 	for (NumbMap::const_iterator iter_n = freeSpecies_.begin(); iter_n != freeSpecies_.end(); ++iter_n) {
 		Int_t species_j(0);
 		const TString& specName = iter_n->first;
 		Double_t value = iter_n->second;
 		cN_[exclName][specName] = value;
 		for (NumbMap::const_iterator iter_j = freeSpecies_.begin(); iter_j != freeSpecies_.end(); ++iter_j) {
 			if (covmat) {
 				cN_[exclName][specName] -= covmat[species_n*nSpecies_+species_j];
 			} else {
 				cN_[exclName][specName] -= covMat_(species_n,species_j);
 			}
 			++species_j;
 		}
 		++species_n;
 	}
 }
 
 void LauSPlot::calcSWeights(const TString& exclName, Double_t *covmat)
 {
 	// Computes the sWeights from the covariance matrix
 	// NB for the extended sPlot the sum in the denominator is still over all species,
 	// while that in the numerator is only over the free species.
 	// Similarly the sWeights can only be calculated for the free species.
 
 	Double_t numerator(0.0), denominator(0.0);
 	for (Int_t iEvent = 0; iEvent < nEvents_; ++iEvent) {
 		denominator = 0.0;
 		for (NumbMap::const_iterator free_iter = freeSpecies_.begin(); free_iter != freeSpecies_.end(); ++free_iter) {
 			denominator += free_iter->second * pdfTot_[iEvent][free_iter->first];
 		}
 		for (NumbMap::const_iterator fixd_iter = fixdSpecies_.begin(); fixd_iter != fixdSpecies_.end(); ++fixd_iter) {
 			denominator += fixd_iter->second * pdfTot_[iEvent][fixd_iter->first];
 		}
 		Int_t species_n(0);
 		for (NumbMap::const_iterator iter_n = freeSpecies_.begin(); iter_n != freeSpecies_.end(); ++iter_n) {
 			numerator = 0.0;
 			Int_t species_j(0);
 			for (NumbMap::const_iterator iter_j = freeSpecies_.begin(); iter_j != freeSpecies_.end(); ++iter_j) {
 				if (covmat) {
 					numerator += covmat[species_n*nSpecies_+species_j] * pdfTot_[iEvent][iter_j->first];
 				} else {
 					numerator += covMat_(species_n,species_j) * pdfTot_[iEvent][iter_j->first];
 				}
 				++species_j;
 			}
 			sWeights_[iEvent][exclName][iter_n->first] = numerator/denominator;
 			++species_n;
 		}
 	}
 }
 
 void LauSPlot::fillCNBranches()
 {
 	if (!cnTree_) {
 		std::cerr<<"ERROR in LauSPlot::fillCNBranches : Tree not created, cannot fill branches."<<std::endl;
 		return;
 	} else if (!this->definedCNBranches()) {
 		this->defineCNBranches();
 	}
 	cnTree_->Fill();
 }
 
 void LauSPlot::fillSWeightBranches()
 {
 	if (!sweightTree_) {
 		std::cerr<<"ERROR in LauSPlot::fillSWeightBranches : Tree not created, cannot fill branches."<<std::endl;
 		return;
 	} else if (sWeights_.empty()) {
 		std::cerr<<"ERROR in LauSPlot::fillSWeightBranches : No sWeights calculated, can't fill branches."<<std::endl;
 		return;
 	} else if (!this->definedSWeightBranches()) {
 		this->copyEventWeights(0);
 		this->defineSWeightBranches();
 	}
 
 	for (Int_t iEvent = 0; iEvent < nEvents_; ++iEvent) {
 		this->copyEventWeights(iEvent);
 		sweightTree_->Fill();
 	}
 }
 
 void LauSPlot::copyEventWeights(Int_t iEvent)
 {
 	for (std::map<TString,NumbMap>::const_iterator excl_iter = sWeights_[iEvent].begin(); excl_iter != sWeights_[iEvent].end(); ++excl_iter) {
 		const TString& exclName = excl_iter->first;
 		const NumbMap& species = excl_iter->second;
 		for (NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
 			const TString& specName = spec_iter->first;
 			sWeightsCurrent_[exclName][specName] = spec_iter->second;
 		}
 	}
 }
 
 void LauSPlot::writeOutResults()
 {
 	// write out the results
 
 	// remove the transient objects that we don't want (re-)written to the file
 	if (eventList_) {
 		delete eventList_; eventList_ = 0;
 	}
 	if (inputTree_) {
 		delete inputTree_; inputTree_ = 0;
 	}
 
 	// first add the input tree as a friend of the output tree
 	this->addFriendTree();
 
 	// then write everything to the file and clean up
 	file_->cd();
 	file_->Write();
 	file_->Close();
 	delete file_; file_ = 0;
 }
 
 void LauSPlot::addFriendTree()
 {
 	if (!sweightTree_) {
 		std::cerr<<"ERROR in LauSPlot::addFriendTree : Tree not created, cannot add friend."<<std::endl;
 		return;
 	}
 	sweightTree_->AddFriend(inputTreeName_,fileName_);
 }
 
 void Yields(Int_t &, Double_t *, Double_t &f, Double_t *x, Int_t /*iflag*/)
 {
 	// FCN-function for Minuit
 
 	f = 0.0;
 
 	TVirtualFitter *fitter = TVirtualFitter::GetFitter();
 	LauSPlot* theModel = dynamic_cast<LauSPlot*>(fitter->GetObjectFit());
 	const std::vector<LauSPlot::NumbMap>& pdfTot = theModel->totalPdf();
 
 	Double_t ntot(0.0);
 	for (std::vector<LauSPlot::NumbMap>::const_iterator evt_iter = pdfTot.begin(); evt_iter != pdfTot.end(); ++evt_iter) {  // loop over events
 		const LauSPlot::NumbMap& species = (*evt_iter);
 		Int_t ispecies(0);
 		Double_t lik(0.0);
 		ntot = 0.0;
 		for (LauSPlot::NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {  // loop over species
 			lik += x[ispecies] * spec_iter->second;
 			ntot += x[ispecies];
 			++ispecies;
 		}
 		if (lik < 0.0) {
 			// make f the worst possible value to force MINUIT
 			// out of this region of parameter space
 			f = -DBL_MAX;
 			break;
 		} else {
 			f += TMath::Log(lik);
 		}
 	}
 
 	// extended likelihood
 	f = (ntot-f);
 }