diff --git a/examples/GenFitTimeDep_Bd2Dpipi.cc b/examples/GenFitTimeDep_Bd2Dpipi.cc
index 2528429..2f84c37 100644
--- a/examples/GenFitTimeDep_Bd2Dpipi.cc
+++ b/examples/GenFitTimeDep_Bd2Dpipi.cc
@@ -1,1054 +1,1021 @@
 
 #include <algorithm>
 #include <filesystem>
 #include <iostream>
 #include <map>
 #include <tuple>
 #include <vector>
 
 #include "TFile.h"
 #include "TH2.h"
 #include "TRandom.h"
 #include "TString.h"
 #include "TSystem.h"
 #include "TF1.h"
 #include "TCanvas.h"
 
 #include "Lau1DHistPdf.hh"
 #include "Lau1DCubicSpline.hh"
 #include "LauBinnedDecayTimeEfficiency.hh"
 #include "LauBkgndDPModel.hh"
 #include "LauCrystalBallPdf.hh"
 #include "LauDaughters.hh"
 #include "LauDecayTime.hh"
 #include "LauDecayTimePdf.hh"
 #include "LauDecayTimePhysicsModel.hh"
 #include "LauDecayTimeResolution.hh"
 #include "LauEffModel.hh"
 #include "LauExponentialPdf.hh"
 #include "LauFlavTag.hh"
 #include "LauGounarisSakuraiRes.hh"
 #include "LauIsobarDynamics.hh"
 #include "LauLHCbNR.hh"
 #include "LauMagPhaseCoeffSet.hh"
 #include "LauRandom.hh"
 #include "LauRealImagCoeffSet.hh"
 #include "LauResonanceMaker.hh"
 #include "LauSplineDecayTimeEfficiency.hh"
 #include "LauSumPdf.hh"
 #include "LauTimeDepFitModel.hh"
 #include "LauVetoes.hh"
 
 #include "TimeDep_Dpipi_ProgOpts.hh"
 
 // Helper function for creating Double-Crystal-Ball PDFs
 LauAbsPdf* makeDCBPDF( const TString& varName, const Double_t varMin, const Double_t varMax, const TString& componentName, const std::map<TString,Double_t> paramVals );
 
 // Helper function for creating sets of DTA splines
 // TODO - not yet used
 using LauCubicSplineDecayTimeEfficiencyPtr = std::unique_ptr<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>>;
 using LauSixthSplineDecayTimeEfficiencyPtr = std::unique_ptr<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>;
 auto createSplineDTASet( const TString& jsonFileName, const TString& controlSplineName, const TString& signalSplineName, const std::vector<TString>& bkgndSplineNames )
 {
 	auto controlSpline = Lau1DCubicSpline::readFromJson( jsonFileName, controlSplineName );
 	LauCubicSplineDecayTimeEfficiencyPtr controlEff { std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>>( controlSplineName, std::move( controlSpline ) ) };
 
 	LauSixthSplineDecayTimeEfficiencyPtr signalEff;
 	if ( signalSplineName != "" ) {
 		auto signalSpline = Lau1DCubicSpline::readFromJson( jsonFileName, signalSplineName );
 		signalEff = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( signalSplineName, *controlEff, std::move( signalSpline ) );
 	}
 
 	std::vector<LauSixthSplineDecayTimeEfficiencyPtr> bkgndEffs;
 	bkgndEffs.reserve( bkgndSplineNames.size() );
 	for ( auto& bkgndSplineName : bkgndSplineNames ) {
 		auto bkgndSpline = Lau1DCubicSpline::readFromJson( jsonFileName, bkgndSplineName );
 		bkgndEffs.push_back( std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( bkgndSplineName, *controlEff, std::move( bkgndSpline ) ) );
 	}
 
 	return std::make_tuple( std::move(controlEff), std::move(signalEff), std::move(bkgndEffs) );
 }
 
 class BgTagDecayFlvAsym {
 	private:
 		const Double_t eff_decayB0_taggedB0       {0.};
 		const Double_t eff_decayB0_taggedB0bar    {0.};
 		const Double_t eff_decayB0bar_taggedB0    {0.};
 		const Double_t eff_decayB0bar_taggedB0bar {0.};
 		const Double_t eff_decayB0_ave            {0.};
 		const Double_t eff_decayB0bar_ave         {0.};
 		const Double_t eff_decayB0_delta          {0.};
 		const Double_t eff_decayB0bar_delta       {0.};
 
 	public:
 		BgTagDecayFlvAsym() = default;
 
 		BgTagDecayFlvAsym( //This constructor assumes a tagging efficincy of 100%
 			const Double_t n_decayB0_taggedB0,
 			const Double_t n_decayB0_taggedB0bar,
 			const Double_t n_decayB0bar_taggedB0,
 			const Double_t n_decayB0bar_taggedB0bar
 		):
 			eff_decayB0_taggedB0      {n_decayB0_taggedB0   /(n_decayB0_taggedB0+n_decayB0_taggedB0bar)},
 			eff_decayB0_taggedB0bar   {n_decayB0_taggedB0bar/(n_decayB0_taggedB0+n_decayB0_taggedB0bar)},
 			eff_decayB0bar_taggedB0   {n_decayB0bar_taggedB0   /(n_decayB0bar_taggedB0+n_decayB0bar_taggedB0bar)},
 			eff_decayB0bar_taggedB0bar{n_decayB0bar_taggedB0bar/(n_decayB0bar_taggedB0+n_decayB0bar_taggedB0bar)},
 
 			eff_decayB0_ave           {0.5*(eff_decayB0_taggedB0    + eff_decayB0_taggedB0bar)},
 			eff_decayB0bar_ave        {0.5*(eff_decayB0bar_taggedB0 + eff_decayB0bar_taggedB0bar)},
 			eff_decayB0_delta         {eff_decayB0_taggedB0    - eff_decayB0_taggedB0bar},
 			eff_decayB0bar_delta      {eff_decayB0bar_taggedB0 - eff_decayB0bar_taggedB0bar}
 		{}
 
 		std::pair<Double_t,Double_t> decayB0()       const {return std::make_pair(eff_decayB0_taggedB0   ,eff_decayB0_taggedB0bar   );}
 		std::pair<Double_t,Double_t> decayB0bar()    const {return std::make_pair(eff_decayB0bar_taggedB0,eff_decayB0bar_taggedB0bar);}
 		std::pair<Double_t,Double_t> aveDeltaB0()    const {return std::make_pair(eff_decayB0_ave   , eff_decayB0_delta)   ;}
 		std::pair<Double_t,Double_t> aveDeltaB0bar() const {return std::make_pair(eff_decayB0bar_ave, eff_decayB0bar_delta);}
 
 	//TODO add constructor for if tagging eff isn't 100%?
 };
 
 
 int main(const int argc, const  char ** argv)
 {
 	const TestDpipi_ProgramSettings settings{argc,argv};
 	if ( ! settings.parsedOK ) {
 		return EXIT_FAILURE;
 	}
 
 	if ( settings.helpRequested ) {
 		return EXIT_SUCCESS;
 	}
 
 	if(settings.dataFit and settings.command == Command::Generate)
 	{
 		std::cerr << "I can't generate data! Wait for Run 3" << std::endl;
 		return EXIT_FAILURE;
 	}
 	if(settings.dataFit and not settings.blindFit and not settings.fixPhiMix)
 	{
 		std::cerr << "We don't have permission to do unblinded fits yet!" << std::endl;
 		return EXIT_FAILURE;
 	}
 	if(settings.run == 0 or settings.run > 2)
 	{
 		std::cerr << "The Run number must be either 1 or 2" << std::endl;
 		return EXIT_FAILURE;
 	}
 
 	const UInt_t runNo = settings.run;
 	LauRandom::setSeed(settings.RNGseed);
 
 	LauDaughters* daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0",    kTRUE);
 	LauDaughters* daughtersB0    = new LauDaughters("B0",     "pi+", "pi-", "D0_bar",kTRUE);
 
 	// efficiency
 	auto vetoes { LauVetoes::readFromJson( "Bd2D0pipi_DP_Vetoes.json", "Bd2D0pipi_Vetoes" ) };
 	LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes.get());
 	LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes.get());
 
 	TFile* effFile(nullptr);
 	TH2* effHist(nullptr);
 
 	if (settings.directory == "B2Dhh_Kpi"){
 		effFile = TFile::Open(Form("%s/%s/eff_map_mprime_thprime_model_d2kpi_run%u.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str(),runNo));
 		effHist = dynamic_cast<TH2*>(effFile->Get(Form("eff_map_mprime_thprime_model_d2kpi_run%u",runNo)));
 	} else if (settings.directory == "B2Dhh_KK"){
 		effFile = TFile::Open(Form("%s/%s/eff_map_mprime_thprime_model_d2kk_run%u.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str(),runNo));
 		effHist = dynamic_cast<TH2*>(effFile->Get(Form("eff_map_mprime_thprime_model_d2kk_run%u",runNo)));
 	} else if (settings.directory == "B2Dhh_pipi"){
 		effFile = TFile::Open(Form("%s/%s/eff_map_mprime_thprime_model_d2pipi_run%u.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str(),runNo));
 		effHist = dynamic_cast<TH2*>(effFile->Get(Form("eff_map_mprime_thprime_model_d2pipi_run%u",runNo)));
 	} else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;}
 
 	effModelB0->setEffHisto(effHist, kTRUE, kFALSE, 0.0, 0.0, kFALSE, kTRUE);
 	effModelB0bar->setEffHisto(effHist, kTRUE, kFALSE, 0.0, 0.0, kFALSE, kTRUE);
 
 	// background types
 	std::map<TString,LauFlavTag::BkgndType> BkgndTypes {
 		{"comb",    LauFlavTag::BkgndType::Combinatorial},
 		{"Bs2DKpi", LauFlavTag::BkgndType::FlavourSpecific},
 		{"Bd2DKpi", LauFlavTag::BkgndType::FlavourSpecific},
 		{"Lb2Dppi", LauFlavTag::BkgndType::FlavourSpecific}
 	};
 	auto notInBkgndList = [&settings](const std::map<TString,LauFlavTag::BkgndType>::value_type& e) {
 		return ( settings.bkgndList.find( e.first.Data() ) == settings.bkgndList.end() );
 	};
 	std::vector<TString> toRemove;
 	for ( const auto& elem : BkgndTypes ) {
 		if ( notInBkgndList(elem) ) {
 			toRemove.push_back( elem.first );
 		}
 	}
 	for ( const auto& name : toRemove ) {
 		BkgndTypes.erase( name );
 	}
 
 	// setup flavour tagging
 	const Bool_t useAveDeltaCalibVals { settings.useAveDeltaCalibVals };
 	const Bool_t useEtaPrime { kFALSE };
 	//LauFlavTag* flavTag = new LauFlavTag(useAveDeltaCalibVals,useEtaPrime);
 	LauFlavTag* flavTag = new LauFlavTag(useAveDeltaCalibVals,useEtaPrime,BkgndTypes);
 	if (settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS) {
 		flavTag->setDecayFlvVarName("decayFlv");
 	}
 	if(not settings.dataFit)
 	{
 		flavTag->setTrueTagVarName("trueTag");
 	}
 
 	TFile* etaFile = TFile::Open(Form("%s/%s/ETAhists.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str()));
 	TH1* etaHist{nullptr};
 	TH1* etaHistB0{nullptr};
 	TH1* etaHistB0bar{nullptr};
 
 	if (settings.directory == "B2Dhh_Kpi"){
 		etaHistB0 = dynamic_cast<TH1*>(etaFile->Get(Form("signalMC_Kpi_Run%u_etaHist_dcyB",runNo)));
 		etaHistB0bar = dynamic_cast<TH1*>(etaFile->Get(Form("signalMC_Kpi_Run%u_etaHist_dcyBbar",runNo)));
 		etaHist = dynamic_cast<TH1*>(etaHistB0->Clone(Form("signalMC_Kpi_Run%u_etaHist",runNo)));
 		etaHist->Add(etaHistB0bar);
 	} else if (settings.directory == "B2Dhh_KK"){
 		etaHist = dynamic_cast<TH1*>(etaFile->Get(Form("signalMC_KK_Run%u_etaHist",runNo)));
 	} else if (settings.directory == "B2Dhh_pipi"){
 		etaHist = dynamic_cast<TH1*>(etaFile->Get(Form("signalMC_pipi_Run%u_etaHist",runNo)));
 	} else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;}
 
 	// Crude check as to whether we're doing perfect vs realistic mis-tag
 	// - in the former case all entries should be in the first bin
 	// If the tagging is perfect then don't interpolate the eta histogram
 	// and also make it perfectly efficient, otherwise do interpolate and
 	// make it 50% efficient
 	Lau1DHistPdf* etaHistPdf = new Lau1DHistPdf( "eta", etaHist, 0.0, 0.5, kTRUE, kFALSE );
 	const Double_t meanEta { etaHistPdf->getMean() };
 
 	const Double_t tagEffVal{1.0};
 	std::pair<Double_t,Double_t> tagEff {tagEffVal, useAveDeltaCalibVals ? 0.0 : tagEffVal};
 
 	// use a null calibration for the time being, so p0 = <eta> and p1 = 1
 	std::pair<Double_t,Double_t> calib0 {meanEta, useAveDeltaCalibVals ? 0.0 : meanEta};
 	std::pair<Double_t,Double_t> calib1 {1.0, useAveDeltaCalibVals ? 0.0 : 1.0};
 
 	if(settings.addTaggers) {
 		flavTag->addTagger(Form("IFT_Run%u",runNo), "IFT_TAGDEC", "IFT_TAGETA", etaHistPdf, tagEff, calib0, calib1);
 		if (settings.floatCalibPars) {
 			flavTag->floatAllCalibPars();
 		}
 	}
 
 	// flavour tagging setup for backgrounds
 	std::map< TString, std::pair<Lau1DHistPdf*,Lau1DHistPdf*> > BkgndEtas;
 	TH1* bkgndEtaHist1{nullptr};
 	Lau1DHistPdf* bkgndEtaHistPDF1{nullptr};
 	TH1* bkgndEtaHistMinus1{nullptr};
 	Lau1DHistPdf* bkgndEtaHistPDFMinus1{nullptr};
 
 	for(auto& [name, _] : BkgndTypes)
 	{
 		const char* decayName = name == "comb" ? "sideband" : name.Data();
 		if (settings.directory == "B2Dhh_Kpi"){
 			bkgndEtaHist1      = dynamic_cast<TH1*>(etaFile->Get( Form("%s_Kpi_Run%u_etaHist_dcyB" ,decayName,runNo) ));
 			bkgndEtaHistPDF1      = new Lau1DHistPdf("eta",bkgndEtaHist1     ,0.0,0.5,kTRUE,kFALSE);
 			bkgndEtaHistMinus1 = dynamic_cast<TH1*>(etaFile->Get( Form("%s_Kpi_Run%u_etaHist_dcyBbar",decayName,runNo) ));
 			bkgndEtaHistPDFMinus1 = new Lau1DHistPdf("eta",bkgndEtaHistMinus1,0.0,0.5,kTRUE,kFALSE);
 		} else if (settings.directory == "B2Dhh_KK"){
 			bkgndEtaHist1      = dynamic_cast<TH1*>(etaFile->Get( Form("%s_KK_Run%u_etaHist" ,decayName,runNo) ));
 			bkgndEtaHistPDF1      = new Lau1DHistPdf("eta",bkgndEtaHist1     ,0.0,0.5,kTRUE,kFALSE);
 		} else if (settings.directory == "B2Dhh_pipi"){
 			bkgndEtaHist1      = dynamic_cast<TH1*>(etaFile->Get( Form("%s_pipi_Run%u_etaHist" ,decayName,runNo) ));
 			bkgndEtaHistPDF1      = new Lau1DHistPdf("eta",bkgndEtaHist1     ,0.0,0.5,kTRUE,kFALSE);
 		} else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;}
 
 		BkgndEtas.emplace( std::make_pair(name, std::make_pair(bkgndEtaHistPDF1,bkgndEtaHistPDFMinus1) ) );
 	}
 
 	for(auto& [name,hists] : BkgndEtas)
 	{
 		if (settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS){
 			if(name == "comb")
 			{
 				const BgTagDecayFlvAsym bgt =
 					runNo == 1 ? BgTagDecayFlvAsym(535.,576.,456.,588.) : BgTagDecayFlvAsym(1514.,1497.,1220.,1575.);
 				if(useAveDeltaCalibVals)
 				{
 					flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo),
 							hists.first  , bgt.aveDeltaB0(),
 							hists.second , bgt.aveDeltaB0bar() );
 				}
 				else
 				{
 					flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo),
 							hists.first  , bgt.decayB0(),
 							hists.second , bgt.decayB0bar());
 				}
 			}
 			else{flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hists.first, tagEff, hists.second, tagEff );}
 			//flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hist, name == "comb" ? std::make_pair(0.5,0.038) : tagEff );
 		} else if (settings.dType == LauTimeDepFitModel::CPEigenvalue::CPEven){
 			if(name == "comb")
 			{
 				BgTagDecayFlvAsym* bgt = nullptr;
 				if(settings.directory == "B2Dhh_KK")
 				{
 					bgt = runNo == 1 ? new BgTagDecayFlvAsym(344.,442.,344.,442.) : new BgTagDecayFlvAsym(710.,814.,710.,814.);
 				}
 				else
 				{
 					bgt = runNo == 1 ? new BgTagDecayFlvAsym(55.,68.,55.,68.)     : new BgTagDecayFlvAsym(210.,241.,210.,241.);
 				}
 
 				std::pair<Double_t,Double_t> effs{bgt->decayB0()};
 				if(useAveDeltaCalibVals){effs = bgt->aveDeltaB0();}
 				flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hists.first, effs);
 				delete bgt;
 			}
 			else{flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hists.first, tagEff);}
 //			else{flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hists.first, tagEff, hists.second, tagEff );} FIXME once we have hists split by decay flv, uncomment this line and remove the above... actually don't FIXME, this is fine as above since the data isn't flavour specific in the CP modes. Keep it combined
 		}
 	}
 
 	// signal dynamics
 	if ( not std::filesystem::exists("Bd2D0pipi_DP_Model_Coeffs.json") ) {
 		std::cerr << "Warning : couldn't find the json file of DP model coeffs; attempting to fetch from eos ..." << std::endl;
 		int r = system( Form("xrdcp %s/%s/Bd2D0pipi_DP_Model_Coeffs.json .", settings.eosRoot.c_str(), settings.eosConfigDir.c_str()) );
 		if (r != 0) { std::cerr << "FATAL : couldn't get the json file from remote!" << std::endl; return 1;}
 	}
 
 	LauIsobarDynamics* sigModelB0bar = new LauIsobarDynamics(daughtersB0bar, effModelB0bar);
 	if (settings.directory == "B2Dhh_Kpi"){
 		sigModelB0bar->setIntFileName(Form("integ_Kpi_Run%u_B0bar.dat",runNo));
 	} else if (settings.directory == "B2Dhh_KK"){
 		sigModelB0bar->setIntFileName(Form("integ_KK_Run%u_B0bar.dat",runNo));
 	} else if (settings.directory == "B2Dhh_pipi"){
 		sigModelB0bar->setIntFileName(Form("integ_pipi_Run%u_B0bar.dat",runNo));
 	}
 	sigModelB0bar->constructModelFromJson( "Bd2D0pipi_DP_Model_Coeffs.json", "Bd2D0pipi_B0bar_Model" );
 
 	LauIsobarDynamics* sigModelB0 = new LauIsobarDynamics(daughtersB0, effModelB0);
 	if (settings.directory == "B2Dhh_Kpi"){
 		sigModelB0->setIntFileName(Form("integ_Kpi_Run%u_B0.dat",runNo));
 	} else if (settings.directory == "B2Dhh_KK"){
 		sigModelB0->setIntFileName(Form("integ_KK_Run%u_B0.dat",runNo));
 	} else if (settings.directory == "B2Dhh_pipi"){
 		sigModelB0->setIntFileName(Form("integ_pipi_Run%u_B0.dat",runNo));
 	}
 	sigModelB0->constructModelFromJson( "Bd2D0pipi_DP_Model_Coeffs.json", "Bd2D0pipi_B0_Model" );
 
 	// fit model
 	LauTimeDepFitModel* fitModel = new LauTimeDepFitModel(sigModelB0bar,sigModelB0,flavTag);
 	fitModel->breakOnBadLikelihood(kFALSE);
 	fitModel->sendFixedYieldsToFitter(kTRUE);
 
 	std::vector<std::unique_ptr<LauAbsCoeffSet>> coeffset { LauAbsCoeffSet::readFromJson( "Bd2D0pipi_DP_Model_Coeffs.json", "Bd2D0pipi_Coeffs" ) };
 
 	for ( auto& coeff : coeffset ) {
 		fitModel->setAmpCoeffSet( std::move(coeff) );
 	}
 
 	// background DP models
 	LauBkgndDPModel* BkgndModel{nullptr};
 	LauBkgndDPModel* BkgndModelbar{nullptr};
 	TFile* BkgndHistFile = TFile::Open(Form("%s/%s/BG_SDPs.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str()));
 	TFile* BkgndHistFile2 = TFile::Open(Form("%s/%s/sidebandSDP.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str()));
 	TH2* BkgndHistB0{nullptr};
 	TH2* BkgndHistB0bar{nullptr};
 	for(auto& [name, _] : BkgndTypes)
 	{
 		//Only have Kpi histos for now for misid
 		if( name=="comb"){
 			if( settings.directory=="B2Dhh_Kpi"){
 				BkgndHistB0 = dynamic_cast<TH2*>(BkgndHistFile2->   Get(Form("sideband_Run%u_Kpi_SDP_decayFlv1",runNo)));
 				BkgndHistB0bar = dynamic_cast<TH2*>(BkgndHistFile2->Get(Form("sideband_Run%u_Kpi_SDP_decayFlvMinus1",runNo)));
 			} else if( settings.directory=="B2Dhh_KK"){
 				BkgndHistB0 = dynamic_cast<TH2*>(BkgndHistFile2->Get(Form("sideband_Run%u_KK_SDP",runNo)));
 			} else if( settings.directory=="B2Dhh_pipi"){
 				BkgndHistB0 = dynamic_cast<TH2*>(BkgndHistFile2->Get(Form("sideband_Run%u_pipi_SDP",runNo)));
 			}
 		} else {
 			TString mode = settings.directory;
 			mode.ReplaceAll("B2Dhh_","");
 			TString histName = Form("%s_%s_Run%u_SDP_decFlv",name.Data(),mode.Data(),runNo);
 			BkgndHistB0 = dynamic_cast<TH2*>(BkgndHistFile->Get(histName+"1"));
 			BkgndHistB0bar = dynamic_cast<TH2*>(BkgndHistFile->Get(histName+"Minus1"));
 		}
 
 		BkgndModel = new LauBkgndDPModel(daughtersB0, vetoes.get());
 		BkgndHistB0->Smooth();
 		BkgndModel->setBkgndSpline(BkgndHistB0, kFALSE, kFALSE, kTRUE );
 
 		//For CP comb we don't split based on decay flavour since it's impossible
 		if(name == "comb" and settings.directory!="B2Dhh_Kpi")
 		{BkgndModelbar=nullptr;}
 		else
 		{
 			BkgndModelbar = new LauBkgndDPModel(daughtersB0bar, vetoes.get());
 			BkgndHistB0bar->Smooth();
 			BkgndModelbar->setBkgndSpline(BkgndHistB0bar, kFALSE, kFALSE, kTRUE );
 		}
 		fitModel->setBkgndDPModels( name, BkgndModel, BkgndModelbar );
 	}
 
 	// decay type and mixing parameter
 	const Bool_t fixPhiMix{ settings.fixPhiMix || settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS };
 	const Bool_t useSinCos{ settings.useSinCos };
 
 	fitModel->setCPEigenvalue( settings.dType );
 	fitModel->setPhiMix( 2.0*LauConstants::beta, fixPhiMix, useSinCos );
 
 	if ( useSinCos ) {
 		if ( settings.blindFit ) {
 			//Blinding strings generated from: https://www.thewordfinder.com/random-sentence-generator/
 			const TString sinBlindString{settings.dataFit ? "Tom realized he could be making a big mistake." : "Don't Look a Gift Horse In The Mouth"};
 			const TString cosBlindString{settings.dataFit ? "He admitted defeat." : "Long In The Tooth"};
 			fitModel->blindPhiMix( sinBlindString, cosBlindString );
 		}
 		std::cout << "Is sin(phiMix) blinded? " << std::boolalpha << fitModel->sinPhiMixBlinded() << std::endl;
 		std::cout << "Is cos(phiMix) blinded? " << std::boolalpha << fitModel->cosPhiMixBlinded() << std::endl;
 	} else {
 		if ( settings.blindFit ) {
 			const TString blindString{"A Cut Above"};
 			fitModel->blindPhiMix( blindString );
 		}
 		std::cout << "Is phiMix blinded? " << std::boolalpha << fitModel->phiMixBlinded() << std::endl;
 	}
 
 	// production asymmetries
 	fitModel->setAsymmetries( 0.0, kTRUE );
 	for(auto& [name, _] : BkgndTypes) {
 		fitModel->setBkgndAsymmetries( name, 0.0, kTRUE );
 	}
 
 	// decay time PDFs
 	const Double_t minDt{0.2};
 	const Double_t maxDt{15.0};
 	const Double_t minDtErr{0.0};
 	const Double_t maxDtErr{0.15};
 
 
 	// physics models
 	auto dtPhysicsModels { LauDecayTimePhysicsModel::readFromJson( "Bd2D0pipi_DecayTime.json", "PhysicsModels" ) };
 
 	auto& dtPhysModel { dtPhysicsModels.at("signal") };
-	auto lifetimePar = dynamic_cast<LauParameter*>( dtPhysModel->findParameter("dt_tau_Bd") );
+	auto lifetimePar = dynamic_cast<LauParameter*>( dtPhysModel->findParameter("dt_sig_tau_Bd") );
 	if ( lifetimePar ) {
 		lifetimePar->fixed( settings.fixLifetime );
 	}
-	auto deltaMPar = dynamic_cast<LauParameter*>( dtPhysModel->findParameter("dt_deltaM_Bd") );
+	auto deltaMPar = dynamic_cast<LauParameter*>( dtPhysModel->findParameter("dt_sig_deltaM_Bd") );
 	if ( deltaMPar ) {
 		deltaMPar->fixed( settings.fixDeltaM );
 	}
 
-	const std::vector<bool> scale {
-		false, false, false
-	//	settings.perEventTimeErr && true,
-	//	settings.perEventTimeErr && true,
-	};
-	const std::size_t nGauss{scale.size()};
-
-	LauParameter * mean0  = new LauParameter("dt_mean_2",  runNo == 1 ? -0.00174035 : -0.00218543 , -0.01, 0.01, kTRUE );
-	LauParameter * mean1  = new LauParameter("dt_mean_1",  runNo == 1 ?  0.043133   :  0.046862   , -0.10, 0.10, kTRUE );
-	LauParameter * mean2  = new LauParameter("dt_mean_0",  runNo == 1 ? -0.00108159 : -0.00159348 , -0.01, 0.01, kTRUE );
-	LauParameter * sigma0 = new LauParameter("dt_sigma_2", runNo == 1 ?  0.046272   :  0.046064   ,  0.0,  2.0,  kTRUE );
-	LauParameter * sigma1 = new LauParameter("dt_sigma_1", runNo == 1 ?  0.11815    :  0.11870    ,  0.0,  2.5,  kTRUE );
-	LauParameter * sigma2 = new LauParameter("dt_sigma_0", runNo == 1 ?  0.025127   :  0.025077   ,  0.0,  2.5,  kTRUE );
-	LauParameter * frac0  = new LauParameter("dt_frac_2",  runNo == 1 ?  0.44359    :  0.43327    ,  0.0,  1.0,  kTRUE );
-	LauParameter * frac1  = new LauParameter("dt_frac_1",  runNo == 1 ?  0.045545   :  0.045310   ,  0.0,  1.0,  kTRUE );
-
-	std::vector<LauAbsRValue*> dtResoPars {
-		mean0,
-		mean1,
-		mean2,
-		sigma0,
-		sigma1,
-		sigma2,
-		frac0,
-		frac1
-	};
+	// resolution models
+	auto dtResolutionModels { LauDecayTimeResolution::readFromJson( "Bd2D0pipi_DecayTime.json", Form("ResolutionModels_Run%d", runNo) ) };
 
-	auto dtResoModel = std::make_unique<LauDecayTimeResolution>( nGauss, dtResoPars, scale );
+	auto& dtResoModel { dtResolutionModels.at("signal") };
 
 	// Decay time error histogram
 	// (always set this so that it gets generated properly, whether we're using it in the PDF or not)
 	TFile* dtErrFile = TFile::Open(Form("%s/%s/dte-hist.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str()));
 	TH1* dtErrHist = dynamic_cast<TH1*>(dtErrFile->Get("dte_hist"));
 
 	LauDecayTimePdf * dtPdf{nullptr};
 	if ( settings.timeResolution ) {
 		if ( settings.perEventTimeErr ) {
+			std::cerr << "Warning : the resolution model has not been constructed with per-event error scaling, but you have requested it" << std::endl;
 			dtPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, "decayTimeErr", minDtErr, maxDtErr, std::move(dtPhysModel), std::move(dtResoModel), dtErrHist );
 		} else {
 			dtPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(dtPhysModel), std::move(dtResoModel) );
 		}
 	} else {
 		dtPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(dtPhysModel) );
 	}
 
 	// Decay time acceptance histogram
 	TFile* dtaFile = TFile::Open(Form("%s/%s/DTAhists.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str()));
 	TH1* dtaHist{nullptr};
 
 	if (settings.directory == "B2Dhh_Kpi"){
 		dtaHist = dynamic_cast<TH1*>(dtaFile->Get(Form("signalMC_Kpi_Run%u_DTAhist",runNo)));
 	} else if (settings.directory == "B2Dhh_KK"){
 		dtaHist = dynamic_cast<TH1*>(dtaFile->Get(Form("signalMC_KK_Run%u_DTAhist",runNo)));
 	} else if (settings.directory == "B2Dhh_pipi"){
 		dtaHist = dynamic_cast<TH1*>(dtaFile->Get(Form("signalMC_pipi_Run%u_DTAhist",runNo)));
 	} else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;}
 
 	if( not std::filesystem::exists("run1DTAsplines.json") and runNo == 1 )
 	{
 		std::cerr << "Warning : couldn't find the json file of splines in Run 1; attempting to fetch from eos ..." << std::endl;
 		int r = system( Form("xrdcp %s/%s/run1DTAsplines.json .", settings.eosRoot.c_str(), settings.eosConfigDir.c_str()) );
 		if(r != 0){ std::cerr << "FATAL : couldn't get the json file from remote!" << std::endl; return 1;}
 	}
 	if( not std::filesystem::exists("run2DTAsplines.json") and runNo == 2 )
 	{
 		std::cerr << "Warning : couldn't find the json file of splines in Run 2; attempting to fetch from eos ..." << std::endl;
 		int r = system( Form("xrdcp %s/%s/run2DTAsplines.json .", settings.eosRoot.c_str(), settings.eosConfigDir.c_str()) );
 		if(r != 0){ std::cerr << "FATAL : couldn't get the json file from remote!" << std::endl; return 1;}
 	}
 
 	//The below is replaced with the new Lau1DCubicSpline::readFromJson
 	// Create the spline knot positions and
 	// starting Y values, to be fit to dtaHist
 	//const std::vector<Double_t> dtvals  {  0.0,   0.5,   1.0,   2.0,   3.0,   4.0,   7.0,  10.0,  15.0};
 	//const std::vector<Double_t> effvals {0.000, 0.010, 0.022, 0.035, 0.042, 0.050, 0.051, 0.052, 0.055};
 	//const std::vector<Double_t> corvals { 1.10,  1.08,  1.06,  1.04,  1.02,  1.00,  0.98,  0.96,  0.94};
 
 	std::unique_ptr<Lau1DCubicSpline> dtEffSpline = nullptr;
 	Lau1DCubicSpline* dtEffSplinePtr = nullptr;
 	std::unique_ptr<Lau1DCubicSpline> dtCorrSpline = nullptr;
 	Lau1DCubicSpline* dtCorrSplinePtr = nullptr;
 
 	LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>* sigDTA = nullptr;
 
 	std::unique_ptr<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>> dtaModel;
 
 	switch(settings.timeEffModel)
 	{
 		case LauDecayTime::EfficiencyMethod::Spline:
 			{
 				fitModel->setASqMaxValue(0.01);
 				if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) {
 					fitModel->setASqMaxValue(0.009);
 				}
 
 				//auto dtEffSpline = std::make_unique<Lau1DCubicSpline>( dtvals, effvals, Lau1DCubicSpline::SplineType::AkimaSpline );
 
 				TString jsonFilename = Form("run%uDTAsplines.json",runNo);
 				dtEffSpline = Lau1DCubicSpline::readFromJson(jsonFilename, Form("signalMC_Kpi_Run%u_DTAhist",runNo));
 				dtEffSplinePtr = dtEffSpline.get();
 
 				dtaModel = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>>( Form("dteff_QFS_Run%u",runNo), std::move(dtEffSpline) );
 				dynamic_cast<LauParameter*>(dtaModel->getParameters().front())->maxValue(1e-1);
 				//dtaModel->fitToTH1(dtaHist);
 
 				// Set which knots to float and which to fix (at least 1 knot must be fixed, not the first one)
 				// Knots should only be floating if requested AND the B lifetime is fixed!
 				if ( settings.fixSplineKnots or not settings.fixLifetime ) {
 					dtaModel->fixKnots();
 				} else {
 					dtaModel->floatKnots();
 					dtaModel->fixKnot( 0, true );
 					dtaModel->fixKnot( 3, true );
 				}
 				sigDTA = dtaModel.get();
 				if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) {
 					// For the QFS mode we just use the cubic model as it is
 					sigDTA = dtaModel.get();
 					dtPdf->setSplineEfficiency( std::move(dtaModel) );
 				} else {
 					// For the CP modes we modify it using a corrective spline
 					// Json available?
 					//auto dtCorrSpline = std::make_unique<Lau1DCubicSpline>( dtvals, corvals, Lau1DCubicSpline::SplineType::AkimaSpline );
 					std::unique_ptr<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>> dtaCPModel;
 					if (settings.directory == "B2Dhh_KK"){
 						dtCorrSpline = Lau1DCubicSpline::readFromJson(jsonFilename, Form("signalMC_KK_Run%u_DTAhist",runNo));
 						dtCorrSplinePtr = dtCorrSpline.get();
 						dtaCPModel = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( Form("dteff_CPKK_Run%u",runNo), *dtaModel, std::move( dtCorrSpline ) );
 					} else if (settings.directory == "B2Dhh_pipi"){
 						dtCorrSpline = Lau1DCubicSpline::readFromJson(jsonFilename, Form("signalMC_pipi_Run%u_DTAhist",runNo));
 						dtCorrSplinePtr = dtCorrSpline.get();
 						dtaCPModel = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( Form("dteff_CPpipi_Run%u",runNo), *dtaModel, std::move( dtCorrSpline ) );
 					}
 					dtPdf->setSplineEfficiency( std::move(dtaCPModel) );
 				}
 
 				break;
 			}
 
 		case LauDecayTime::EfficiencyMethod::Binned:
 			{
 				fitModel->setASqMaxValue(0.06);
 
 				auto dtaBinnedModel = std::make_unique<LauBinnedDecayTimeEfficiency>( *dtaHist );
 				dtPdf->setBinnedEfficiency( std::move(dtaBinnedModel) );
 
 				break;
 			}
 
 		case LauDecayTime::EfficiencyMethod::Uniform:
 			{
 				fitModel->setASqMaxValue(4.45);
 				break;
 			}
 	}
 
 	fitModel->setSignalDtPdf( dtPdf );
 
 	std::map<TString, LauDecayTime::FuncType> bkgndDecayTimeTypes = {
 		{"comb",    LauDecayTime::FuncType::Hist},
 		{"Bd2DKpi", LauDecayTime::FuncType::ExpTrig},
 		{"Bs2DKpi", LauDecayTime::FuncType::ExpHypTrig},
 		{"Lb2Dppi", LauDecayTime::FuncType::Exp}
 	};
 
 	TH1* bgDTA{nullptr};
 	std::map<TString, std::pair<Double_t, Double_t > > decayFlvYields = {
 		{"comb"   , runNo == 1 ? std::make_pair( 1111, 1044) : std::make_pair( 3011, 2795)},
 		{"Bd2DKpi", runNo == 1 ? std::make_pair( 9552, 8933) : std::make_pair(13842,12036)},
 		{"Bs2DKpi", runNo == 1 ? std::make_pair(16989,15634) : std::make_pair(18294,17005)},
 		{"Lb2Dppi", runNo == 1 ? std::make_pair( 3695, 3431) : std::make_pair( 2389, 2501)}
 	};
 	TString DTAnameSuffix{""};
 
 	if (settings.directory == "B2Dhh_Kpi"){
 		DTAnameSuffix = Form("_Kpi_Run%u_DTAhist",runNo);
 	} else if (settings.directory == "B2Dhh_KK"){
 		DTAnameSuffix = Form("_KK_Run%u_DTAhist",runNo);
 	} else if (settings.directory == "B2Dhh_pipi"){
 		DTAnameSuffix = Form("_pipi_Run%u_DTAhist",runNo);
 	}
 
 	for (auto& [bg, _] : BkgndTypes) {
 
 		if ( bkgndDecayTimeTypes[bg] == LauDecayTime::FuncType::Hist ) {
 
 			TString histname = bg + DTAnameSuffix;
 			if(histname.Contains("comb"))
 			{
 				histname = "Sideband" + DTAnameSuffix;
 				histname.ReplaceAll("DTAhist","DThist");
 			}
 			TH1* dt_hist = dynamic_cast<TH1*>(dtaFile->Get(histname.Data()));
 			TH1* dt_err_hist = nullptr; // TODO get these!
 
 			if ( not dt_hist ) { std::cerr << "PROBLEM FINDING THE LIFETIME HIST:" << histname << std::endl; return EXIT_FAILURE; }
 			if ( settings.perEventTimeErr and not dt_err_hist ) { std::cerr << "PROBLEM FINDING THE LIFETIME ERROR HIST:" << bg << std::endl; return EXIT_FAILURE; }
 
 			LauDecayTimePdf* backgroundPdf {nullptr};
 			if ( settings.perEventTimeErr ) {
 				backgroundPdf = new LauDecayTimePdf("decayTime", minDt, maxDt, "decayTimeErr", minDtErr, maxDtErr, dt_hist, dt_err_hist );
 			} else {
 				backgroundPdf = new LauDecayTimePdf("decayTime", minDt, maxDt, dt_hist );
 			}
 
 			fitModel->setBkgndDtPdf(bg,backgroundPdf);
 
 		} else {
 
 			auto& bgPhysModel { dtPhysicsModels.at(bg) };
-
-			std::vector<LauAbsRValue*> bgResoParams {
-				mean0->createClone(),
-				mean1->createClone(),
-				mean2->createClone(),
-				sigma0->createClone(),
-				sigma1->createClone(),
-				sigma2->createClone(),
-				frac0->createClone(),
-				frac1->createClone()
-			};
-			auto bgResoModel = std::make_unique<LauDecayTimeResolution>( nGauss, bgResoParams, scale );
+			auto& bgResoModel { dtResolutionModels.at(bg) };
 
 			TH1* dt_err_hist = nullptr; // TODO get these!
 			if ( settings.perEventTimeErr and not dt_err_hist ) { std::cerr << "PROBLEM FINDING THE LIFETIME ERROR HIST:" << bg << std::endl; return EXIT_FAILURE; }
 
 			LauDecayTimePdf* backgroundPdf {nullptr};
 			if ( settings.timeResolution ) {
 				if ( settings.perEventTimeErr ) {
+					std::cerr << "Warning : the resolution model has not been constructed with per-event error scaling, but you have requested it" << std::endl;
 					backgroundPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, "decayTimeErr", minDtErr, maxDtErr, std::move(bgPhysModel), std::move(bgResoModel), dt_err_hist );
 				} else {
 					backgroundPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(bgPhysModel), std::move(bgResoModel) );
 				}
 			} else {
 				backgroundPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(bgPhysModel) );
 			}
 
 			// background decay time acceptance histogram
 			TString histname = bg + DTAnameSuffix;
 			bgDTA = dynamic_cast<TH1*>(dtaFile->Get(histname));
 			if ( not bgDTA and bg!="comb" ) {
 				std::cerr << "Couldn't find DTA histogram for background: " << bg << std::endl;
 				return EXIT_FAILURE;
 			}
 
 			//XXX if we decide to rebin, do it here
 			bgDTA->Divide(dtaHist);
 			for (Int_t bin = 1; bin <= bgDTA->GetNbinsX(); ++bin)
 			{
 				if( bgDTA->GetBinContent(bin)==0. ){bgDTA->SetBinContent(bin,1.);}
 			}
 
 			TString jsonFilename = Form("run%iDTAsplines.json",runNo);
 			auto dtBkgndCorrSpline = Lau1DCubicSpline::readFromJson(jsonFilename,histname);
 
 			std::unique_ptr<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>> bgDTAModel;
 			if (sigDTA==nullptr){
 				std::cout << "Null pointer" << std::endl;
 			}
 			if (settings.directory == "B2Dhh_Kpi"){
 				bgDTAModel = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( Form("%s_DTA_ratio_Run%u_Kpi",bg.Data(), runNo), *sigDTA, std::move( dtBkgndCorrSpline ) );
 			} else if (settings.directory == "B2Dhh_KK"){
 				bgDTAModel = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( Form("%s_DTA_ratio_Run%u_KK",bg.Data(), runNo), *sigDTA, std::move( dtBkgndCorrSpline ) );
 			} else if (settings.directory == "B2Dhh_pipi"){
 				bgDTAModel = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( Form("%s_DTA_ratio_Run%u_pipi",bg.Data(),runNo), *sigDTA, std::move( dtBkgndCorrSpline ) );
 			}
 
 			//bgDTAModel->fitToTH1(bgDTA, LauOutputLevel::Quiet, true);
 			backgroundPdf->setSplineEfficiency( std::move(bgDTAModel) );
 
 			fitModel->setBkgndDtPdf(bg,backgroundPdf);
 		}
 	}
 
 	// Signal and bkgnd yields
 	Double_t nSigEvents{0};
 	Double_t nCombEvents{0};
 	std::map<TString, Double_t> bkgndYields;
 	if (settings.directory == "B2Dhh_Kpi") {
 		nSigEvents = (runNo == 1) ? 28362 : 115032;
 		nCombEvents = (runNo == 1) ? 2415 : 7637;
 		if ( ! settings.scaleYields ) {
 			nCombEvents += (runNo == 1) ? (303 + 111 + 123 + 143 + 2) : (974 + 390 + 256 + 389 + 4);
 		}
 		bkgndYields = {{"comb", nCombEvents },
 			  {"Bs2DKpi", runNo == 1 ? 796 : 1713  },
 			  {"Bd2DKpi", runNo == 1 ? 254 : 639   },
 			  {"Lb2Dppi", runNo == 1 ? 245 : 304  }};
 	} else if (settings.directory == "B2Dhh_KK"){
 		if (settings.setNegToZero and runNo == 1) {
 			// Mass fit with mis-ID yields fixed to 0
 			nSigEvents = 3098;
 			nCombEvents = 800;
 			if( ! settings.scaleYields){
 				nCombEvents += 20 + 22 + 53 + 39 + 0;
 			}
 			bkgndYields = {{"comb", nCombEvents },
 				{"Bs2DKpi", 0  },
 				{"Bd2DKpi", 0  },
 				{"Lb2Dppi", 0  }};
 		}
 		else if (settings.constrainedYields and runNo == 1)
 		{
 
 			// Mass fit with mis-ID yields constrained based on
 			// Run1/Run2 in Kpi and pipi and Run2 yield of KK
 			nSigEvents = 3093;
 			nCombEvents = 777;
 			if( ! settings.scaleYields){
 				nCombEvents += 19 + 21 + 52 + 41 + 0;
 			}
 			bkgndYields = {{"comb", nCombEvents },
 				{"Bs2DKpi", 19  },
 				{"Bd2DKpi", 6  },
 				{"Lb2Dppi", 7  }};
 		} else {
 			// Original mass fit
 			nSigEvents = runNo == 1 ? 3140 : 12099;
 			nCombEvents = (runNo == 1) ? 1016 : 1501;
 			if ( ! settings.scaleYields ) {
 				nCombEvents += (runNo == 1) ? (63 + 23 + 26 + 13 + 1) : (107 + 43 + 28 + 5 + 1);
 			}
 			bkgndYields = {{"comb", nCombEvents },
 				{"Bs2DKpi", runNo == 1 ? -172 : 167  },
 				{"Bd2DKpi", runNo == 1 ? -64  : 73   },
 				{"Lb2Dppi", runNo == 1 ? -61  : 33   }};
 		}
 
 	} else if (settings.directory == "B2Dhh_pipi"){
 		nSigEvents = runNo == 1 ? 1175 : 4566;
 		nCombEvents = (runNo == 1) ? 115 : 301;
 		if ( ! settings.scaleYields ) {
 			nCombEvents += (runNo == 1) ? (7 + 3 + 3 + 79 + 0) : (36 + 15 + 10 + 138 + 1);
 		}
 		bkgndYields = {{"comb", nCombEvents },
 			  {"Bs2DKpi", runNo == 1 ? 57 : 99  },
 			  {"Bd2DKpi", runNo == 1 ? 18 : 41   },
 			  {"Lb2Dppi", runNo == 1 ? 21 : 19   }};
 	} else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;}
 
 	if(settings.command == Command::Generate)
 	{
 		for(auto& [mode, yield] : bkgndYields)
 		{
 			if(yield < 0.)
 			{
 				std::cerr << "WARNING : setting negative yield for mode " << mode << " to 0 for generation" << std::endl;
 				yield = 0.;
 			}
 		}
 	}
 
 	const Bool_t fixYields { settings.bkgndList.empty() or not settings.floatYields };
 
 	TString yieldName;
 	for (auto& [bg, _] : BkgndTypes) {
 		const Double_t nBkg = bkgndYields[bg];
 		std::cout<<bg<<" Events = "<<nBkg<<std::endl;
 		yieldName = Form("%sEvents_%s_Run%u", bg.Data(), settings.directory.c_str(), runNo);
 		LauParameter* nBkgndEvents = new LauParameter(yieldName,nBkg,-2.0*TMath::Abs(nBkg),2.0*TMath::Abs(nBkg),fixYields);
 		fitModel->setNBkgndEvents( nBkgndEvents );
 	}
 
 	std::cout<<"nSigEvents = "<<nSigEvents<<std::endl;
 	yieldName = Form("signalEvents_%s_Run%u", settings.directory.c_str(), runNo);
 	LauParameter* nSigPar = new LauParameter(yieldName, nSigEvents, -2.0*nSigEvents, 2.0*nSigEvents, fixYields);
 	fitModel->setNSigEvents(nSigPar);
 
 	if ( settings.command == Command::Generate or not fixYields ) {
 		// mB PDFs
 		// PDFs for B+ and B- are set to be identical here but could be different in general
 		const TString mBName { "mB" };
 		const Double_t mBMin { 5.200 };
 		const Double_t mBMax { 5.700 };
 
 		// Signal PDF is a double Crystal Ball function
 		const std::map<TString,Double_t> sigParVals {
 			{"mB_sig_mean",    5.2802},
 				{"mB_sig_sigma",   0.01125},
 				{"mB_sig_alphaL",  1.226},
 				{"mB_sig_alphaR", -1.722},
 				{"mB_sig_orderL",  1.817},
 				{"mB_sig_orderR",  2.861},
 				{"mB_sig_frac",    0.480}
 		};
 		fitModel->setSignalPdfs( makeDCBPDF( mBName, mBMin, mBMax, "sig", sigParVals ) );
 
 		// Combinatorial PDF is an exponential
 		if ( BkgndTypes.find("comb") != BkgndTypes.end() ) {
 			LauParameter* mB_comb_slope = new LauParameter("mB_comb_slope", -0.00419);
 			LauAbsPdf* mB_comb_pdf = new LauExponentialPdf(mBName, {mB_comb_slope}, mBMin, mBMax);
 			fitModel->setBkgndPdf("comb", mB_comb_pdf);
 		}
 
 		// Bd2DKpi PDF is double Crystal Ball function
 		if ( BkgndTypes.find("Bd2DKpi") != BkgndTypes.end() ) {
 			const std::map<TString,Double_t> Bd2DKpiParVals {
 				{"mB_Bd2DKpi_mean",    5.2401},
 					{"mB_Bd2DKpi_sigma",   0.01580},
 					{"mB_Bd2DKpi_alphaL",  0.560},
 					{"mB_Bd2DKpi_alphaR", -2.212},
 					{"mB_Bd2DKpi_orderL",  1.076},
 					{"mB_Bd2DKpi_orderR",  2.453},
 					{"mB_Bd2DKpi_frac",    0.544}
 			};
 			fitModel->setBkgndPdf("Bd2DKpi", makeDCBPDF( mBName, mBMin, mBMax, "Bd2DKpi", Bd2DKpiParVals ) );
 		}
 
 		// Bs2DKpi PDF is double Crystal Ball function
 		if ( BkgndTypes.find("Bs2DKpi") != BkgndTypes.end() ) {
 			const std::map<TString,Double_t> Bs2DKpiParVals {
 				{"mB_Bs2DKpi_mean",    5.3244},
 					{"mB_Bs2DKpi_sigma",   0.01827},
 					{"mB_Bs2DKpi_alphaL",  0.865},
 					{"mB_Bs2DKpi_alphaR", -0.001},
 					{"mB_Bs2DKpi_orderL", 10.000},
 					{"mB_Bs2DKpi_orderR", 10.000},
 					{"mB_Bs2DKpi_frac",    0.954}
 			};
 			fitModel->setBkgndPdf("Bs2DKpi", makeDCBPDF( mBName, mBMin, mBMax, "Bs2DKpi", Bs2DKpiParVals ) );
 		}
 
 		// Lb2Dppi PDF is double Crystal Ball function
 		if ( BkgndTypes.find("Lb2Dppi") != BkgndTypes.end() ) {
 			const std::map<TString,Double_t> Lb2DppiParVals {
 				{"mB_Lb2Dppi_mean",     5.4884},
 					{"mB_Lb2Dppi_sigma",    0.02561},
 					{"mB_Lb2Dppi_alphaL",   0.091},
 					{"mB_Lb2Dppi_alphaR", -18.053},
 					{"mB_Lb2Dppi_orderL",   2.906},
 					{"mB_Lb2Dppi_orderR",   6.087},
 					{"mB_Lb2Dppi_frac",     0.968}
 			};
 			fitModel->setBkgndPdf("Lb2Dppi", makeDCBPDF( mBName, mBMin, mBMax, "Lb2Dppi", Lb2DppiParVals ) );
 		}
 	}
 
 	// set the number of experiments
 	if (settings.command == Command::Generate) {
 		fitModel->setNExpts(settings.nExptGen, settings.firstExptGen);
 	} else {
 		fitModel->setNExpts(settings.nExptFit, settings.firstExptFit);
 	}
 
 	fitModel->useAsymmFitErrors(kFALSE);
 	fitModel->useRandomInitFitPars(kFALSE);
 	fitModel->writeLatexTable(kFALSE);
 	const Bool_t haveBkgnds = ( fitModel->nBkgndClasses() > 0 );
 	fitModel->doPoissonSmearing(haveBkgnds);
 	fitModel->doEMLFit(haveBkgnds);
 
 	TString dTypeStr;
 	switch (settings.dType) {
 		case LauTimeDepFitModel::CPEigenvalue::CPEven :
 			dTypeStr = "CPEven";
 			break;
 		case LauTimeDepFitModel::CPEigenvalue::CPOdd :
 			dTypeStr = "CPOdd";
 			break;
 		case LauTimeDepFitModel::CPEigenvalue::QFS :
 			dTypeStr = "QFS";
 			break;
 	}
 
 	TString dataFile("");
 	TString treeName("fitTree");
 	TString rootFileName("");
 	TString tableFileName("");
 	TString fitToyFileName("");
 	TString splotFileName("");
 
 	dataFile = "TEST-Dpipi";
 
 	dataFile += "_"+dTypeStr+"_"+settings.directory;
 
 	switch(settings.timeEffModel)
 	{
 		case LauDecayTime::EfficiencyMethod::Spline:
 			dataFile += "_Spline";
 			break;
 		case LauDecayTime::EfficiencyMethod::Binned:
 			dataFile += "_Binned";
 			break;
 		case LauDecayTime::EfficiencyMethod::Uniform:
 			dataFile += "_Uniform";
 			break;
 	}
 
 	if (settings.timeResolution) {
 		if (settings.perEventTimeErr) {
 			dataFile += "_DTRperevt";
 		} else {
 			dataFile += "_DTRavg";
 		}
 	} else {
 		dataFile += "_DTRoff";
 	}
 
 	dataFile += "_expts";
 	dataFile += settings.firstExptGen;
 	dataFile += "-";
 	dataFile += settings.firstExptGen+settings.nExptGen-1;
 	dataFile += Form("-Run%u",runNo);
 
 	dataFile += ".root";
 
 	if (settings.dataFit)
 	{
 		if(runNo == 2)
 		{
 			dataFile = Form("%s/%s/Run2/B2Dhh-CollisionCombined-MagCombined-Stripping24r2-withDNN-withPIDcorr-withMVA-withPIDMVA-WithDpiMatching-withSwappedMassHypotheses-withKKvetoes-withSelection-LauraPrepped.root",settings.eosRoot.c_str(),settings.eosDataDir.c_str());
 		} else {
 			dataFile = Form("%s/%s/Run1/B2Dhh-CollisionCombined-MagCombined-Stripping21r1p2-withDNN-withPIDcorr-withMVA-withPIDMVA-WithDpiMatching-withSwappedMassHypotheses-withKKvetoes-withSelection-LauraPrepped.root",settings.eosRoot.c_str(),settings.eosDataDir.c_str());
 		}
 		treeName = Form("%s/B2DhhReco",settings.directory.data());
 	}
 
 	if (settings.command == Command::Generate) {
 		rootFileName = "dummy.root";
 		tableFileName = "genResults";
 	} else {
 		rootFileName = "fit"; rootFileName += settings.iFit;
 		rootFileName += "_Results_"; rootFileName += dTypeStr;
 		rootFileName += "_"; rootFileName += settings.directory;
 		rootFileName += "_expts"; rootFileName += settings.firstExptFit; rootFileName += "-"; rootFileName += settings.firstExptFit+settings.nExptFit-1;
 		rootFileName += Form("_Run%u",runNo);
 		rootFileName += ".root";
 
 		tableFileName = "fit"; tableFileName += settings.iFit;
 		tableFileName += "_Results_"; tableFileName += dTypeStr;
 		tableFileName += "_"; tableFileName += settings.directory;
 		tableFileName += "_expts"; tableFileName += settings.firstExptFit; tableFileName += "-"; tableFileName += settings.firstExptFit+settings.nExptFit-1;
 		tableFileName += Form("_Run%u",runNo);
 
 		fitToyFileName = "fit"; fitToyFileName += settings.iFit;
 		fitToyFileName += "_ToyMC_"; fitToyFileName += dTypeStr;
 		fitToyFileName += "_"; fitToyFileName += settings.directory;
 		fitToyFileName += "_expts"; fitToyFileName += settings.firstExptFit; fitToyFileName += "-"; fitToyFileName += settings.firstExptFit+settings.nExptFit-1;
 		fitToyFileName += Form("_Run%u",runNo);
 		fitToyFileName += ".root";
 
 		splotFileName = "fit"; splotFileName += settings.iFit;
 		splotFileName += "_sPlot_"; splotFileName += dTypeStr;
 		splotFileName += "_"; splotFileName += settings.directory;
 		splotFileName += "_expts"; splotFileName += settings.firstExptFit; splotFileName += "-"; splotFileName += settings.firstExptFit+settings.nExptFit-1;
 		splotFileName += Form("_Run%u",runNo);
 		splotFileName += ".root";
 	}
 
 	if( settings.genToy and not settings.blindFit )
 	{
 		fitModel->compareFitData(50, fitToyFileName);
 	}
 
 	// Generate toy from the fitted parameters
 	//fitModel->compareFitData(1, fitToyFileName);
 
 	// Write out per-event likelihoods and sWeights
 	//fitModel->writeSPlotData(splotFileName, "splot", kFALSE);
 
 	// Execute the generation/fit
 	switch (settings.command) {
 		case Command::Generate :
 			fitModel->run( "gen", dataFile, treeName, rootFileName, tableFileName );
 			break;
 		case Command::Fit :
 			fitModel->run( "fit", dataFile, treeName, rootFileName, tableFileName );
 			break;
 		case Command::SimFit :
 			fitModel->runTask( dataFile, treeName, rootFileName, tableFileName, "localhost", settings.port );
 			break;
 	}
 
 	return EXIT_SUCCESS;
 }
 
 LauAbsPdf* makeDCBPDF( const TString& varName, const Double_t varMin, const Double_t varMax, const TString& componentName, const std::map<TString,Double_t> paramVals )
 {
 	const TString meanName   { Form( "%s_%s_mean",   varName.Data(), componentName.Data() ) };
 	const TString sigmaName  { Form( "%s_%s_sigma",  varName.Data(), componentName.Data() ) };
 	const TString alphaLName { Form( "%s_%s_alphaL", varName.Data(), componentName.Data() ) };
 	const TString alphaRName { Form( "%s_%s_alphaR", varName.Data(), componentName.Data() ) };
 	const TString orderLName { Form( "%s_%s_orderL", varName.Data(), componentName.Data() ) };
 	const TString orderRName { Form( "%s_%s_orderR", varName.Data(), componentName.Data() ) };
 	const TString fracName   { Form( "%s_%s_frac",   varName.Data(), componentName.Data() ) };
 
 	LauParameter* mean   = new LauParameter( meanName,   paramVals.at( meanName )   );
         LauParameter* sigma  = new LauParameter( sigmaName,  paramVals.at( sigmaName )  );
 	LauParameter* alphaL = new LauParameter( alphaLName, paramVals.at( alphaLName ) );
         LauParameter* alphaR = new LauParameter( alphaRName, paramVals.at( alphaRName ) );
 	LauParameter* orderL = new LauParameter( orderLName, paramVals.at( orderLName ) );
         LauParameter* orderR = new LauParameter( orderRName, paramVals.at( orderRName ) );
 	LauParameter* frac   = new LauParameter( fracName,   paramVals.at( fracName )   );
 
 	std::vector<LauAbsRValue*> pars;
 	pars.reserve(4);
 
 	pars.clear();
         pars.push_back(mean);
         pars.push_back(sigma);
         pars.push_back(alphaL);
         pars.push_back(orderL);
 	LauAbsPdf* pdf_L = new LauCrystalBallPdf(varName, pars, varMin, varMax);
 
 	pars.clear();
         pars.push_back(mean);
         pars.push_back(sigma);
         pars.push_back(alphaR);
         pars.push_back(orderR);
 	LauAbsPdf* pdf_R = new LauCrystalBallPdf(varName, pars, varMin, varMax);
 
 	LauAbsPdf* pdf = new LauSumPdf(pdf_L, pdf_R, frac);
 
 	return pdf;
 }
diff --git a/inc/LauDecayTimeResolution.hh b/inc/LauDecayTimeResolution.hh
index b21d39e..efafa63 100644
--- a/inc/LauDecayTimeResolution.hh
+++ b/inc/LauDecayTimeResolution.hh
@@ -1,194 +1,332 @@
 
 /*
 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 LauDecayTimeResolution.hh
     \brief File containing declaration of LauDecayTimeResolution class.
 */
 
-/*! \class LauDecayTimeResolution
-    \brief Class for defining the model for decay time resolution
-*/
-
 #ifndef LAU_DECAYTIME_RESOLUTION
 #define LAU_DECAYTIME_RESOLUTION
 
-#include <vector>
+#include "LauAbsRValue.hh"
+#include "LauJsonTools.hh"
 
 #include "Rtypes.h"
+#include "TString.h"
 
-class LauAbsRValue;
+#include <nlohmann/json.hpp>
 
+#include <map>
+#include <memory>
+#include <iostream>
+#include <vector>
+
+/*! \class LauDecayTimeResolution
+    \brief Class for defining the model for decay time resolution
+*/
 
 class LauDecayTimeResolution final {
 
 	public:
+		/*! \enum ScalingType
+		    \brief Define the types of per-event error scaling that can be performed
+		*/
+		enum class ScalingType {
+			None,		//!< No scaling is performed
+			Global,		//!< All means and widths are scaled
+			PerGaussian,	//!< Scaling is enabled/disabled for each Gaussian
+			PerParameter	//!< Scaling is enabled/disabled for individual parameters
+		};
+
 		//! Constructor
 		/*!
 		    \param nGauss the number of Gaussians in the resolution model
 		    \param resolutionParams the parameters of the resolution model
 		    \param scale if true all means and widths are scaled by per-event decay time error, if false none are scaled
 		*/
-		LauDecayTimeResolution( const std::size_t nGauss, const std::vector<LauAbsRValue*>& resolutionParams, const bool scale = false );
+		LauDecayTimeResolution( const std::size_t nGauss, std::vector<std::unique_ptr<LauAbsRValue>> resolutionParams, const bool scale = false );
 
 		//! Constructor
 		/*!
 		    \param nGauss the number of Gaussians in the resolution model
 		    \param resolutionParams the parameters of the resolution model
 		    \param scale for each Gaussian, if true its mean and width are both scaled by per-event decay time error, if false neither are scaled
 		*/
-		LauDecayTimeResolution( const std::size_t nGauss, const std::vector<LauAbsRValue*>& resolutionParams, const std::vector<bool>& scale );
+		LauDecayTimeResolution( const std::size_t nGauss, std::vector<std::unique_ptr<LauAbsRValue>> resolutionParams, const std::vector<bool>& scale );
 
 		//! Constructor
 		/*!
 		    \param nGauss the number of Gaussians in the resolution model
 		    \param resolutionParams the parameters of the resolution model
 		    \param scaleMeans for each Gaussian, if true its mean is scaled by per-event decay time error, if false its mean is not scaled
 		    \param scaleWidths for each Gaussian, if true its width is scaled by per-event decay time error, if false its width is not scaled
 		*/
-		LauDecayTimeResolution( const std::size_t nGauss, const std::vector<LauAbsRValue*>& resolutionParams, const std::vector<bool>& scaleMeans, const std::vector<bool>& scaleWidths );
+		LauDecayTimeResolution( const std::size_t nGauss, std::vector<std::unique_ptr<LauAbsRValue>> resolutionParams, const std::vector<bool>& scaleMeans, const std::vector<bool>& scaleWidths );
 
 		//! Retrieve the parameters of the resolution model so that they can be loaded into a fit
 		/*!
 		    \return the parameters of the resolution model
 		*/
 		std::vector<LauAbsRValue*> getParameters() { return params_; }
 
 		//! Retrieve the number of Gaussians in the model
 		/*!
 		    \return the number of Gaussians
 		*/
 		const std::size_t& nGauss() const { return nGauss_; }
 
+		//! Retrieve the scaling type
+		/*!
+		    \return the scaling type
+		*/
+		ScalingType getScalingType() const { return scalingType_; }
+
 		//! Retrieve whether any of the parameters of the resolution function scaled by the per-event error
 		/*!
 		    \return whether scaling or not
 		*/
 		bool scaleWithPerEventError() const { return scaleWithPerEventError_; }
 
 		//! Retrieve whether the mean of each Gaussian is scaled by the per-event decay time error
 		/*!
 		    \return the mean scaling flags
 		*/
 		const std::vector<bool>& scaleMeans() const { return scaleMeans_; }
 
 		//! Retrieve whether the width of each Gaussian is scaled by the per-event decay time error
 		/*!
 		    \return the width scaling flags
 		*/
 		const std::vector<bool>& scaleWidths() const { return scaleWidths_; }
 
 		//! Retrieve the up-to-date values of the means
 		/*!
 		    \return the mean values
 		*/
 		const std::vector<Double_t>& getMeanValues() const { return meanVals_; }
 
 		//! Retrieve the up-to-date values of the widths
 		/*!
 		    \return the width values
 		*/
 		const std::vector<Double_t>& getWidthValues() const { return widthVals_; }
 
 		//! Retrieve the up-to-date values of the fractions
 		/*!
 		    \return the fraction values
 		*/
 		const std::vector<Double_t>& getFractionValues() const { return fractionVals_; }
 
 		//! Retrieve whether any of the parameter values have changed in the last fit iteration
 		/*!
 		    \return the any param changed flag
 		*/
 		const bool& anythingChanged() const { return anythingChanged_; }
 
 		//! Initialise the parameter cache
 		/*!
 			Must be called prior to starting fitting or generation
 		*/
 		void initialise();
 
 		//! Propagate any updates to parameters and recalculate information as neeeded
 		/*!
 			Should be called at each fit iteration
 		*/
 		void propagateParUpdates();
 
-	private:
+		//! Construct a collection of resolution model objects based on values read from a JSON value
+		/*!
+			The JSON value should be an Object, which must contain the following elements:
+			- "parameters", an Array, which contains an Object for each parameter of the models to be constructed
+			- "signalModel", an Object, which contains:
+			  + "nGauss", a Number_Unsigned to specify the number of Gaussian functions in the model
+			  + "scalingType", a String to specify the scaling type (must match a state of ScalingType)
+			  + "parameters", an Array of Strings that specifies the parameters from the top-level "parameters" Array to be used to build the signal model
+			- "backgroundModels", an Array, which contains Objects for each background component, each of which contain:
+			  + "name", a String specifying the background component name
+			  + "model", an Object to define the model (definition as per the "signalModel" entry above)
+
+			\param [in] j the JSON value to deserialise
+			\return the collection of newly constructed models
+		*/
+		static std::map<TString,std::unique_ptr<LauDecayTimeResolution>> readFromJson( const nlohmann::json& j );
+
+		//! Construct a collection of resolution model objects based on values read from a JSON file
+		/*!
+		  	The JSON value structure is as defined in readFromJson(nlohmann::json)
+
+			\param [in] fileName the name of the file from which the JSON should be read
+			\param [in] elementName the optional name of the JSON element that contains the definitions (defaults to using the root record)
+			\return the collection of newly constructed models
+		*/
+		static std::map<TString,std::unique_ptr<LauDecayTimeResolution>> readFromJson( const TString& fileName, const TString& elementName = "" );
+
+		//! Write a collection of physics model objects to a JSON file
+		/*!
+			\tparam ResModelPtr a pointer-like type that points to LauDecayTimeResolution objects
+			\param [in] models the collection of models to be written out
+			\param [in] fileName the name of the file to which the JSON should be written
+			\param [in] elementName the (optional) name of the JSON element to which the coefficients should be written
+			\param [in] append if true, append the spline to the existing JSON within the file - if using this option it is then mandatory to provide elementName
+			\param [in] indent the indentation level to use in the output in number of spaces (defaults to 4)
+		*/
+		template <typename ResModelPtr>
+		static void writeToJson( const std::map<TString,ResModelPtr>& models, const TString& fileName, const TString& elementName = "", const bool append = false, const int indent = 4 );
+
 		//! Retrieve the specified parameter
 		/*!
 			\param [in] parName the parameter to retrieve
 		*/
 		const LauAbsRValue* findParameter(const TString& parName) const;
 
 		//! Retrieve the specified parameter
 		/*!
 			\param [in] parName the parameter to retrieve
 		*/
 		LauAbsRValue* findParameter(const TString& parName);
 
+	private:
 		//! Resize the internal vectors to match nGauss_
 		void resizeVectors();
 
 		//! Do an initial sanity check of our setup
 		void checkSetup();
 
 		//! Update the cached parameter values
 		void updateParameterCache();
 
+		//! Read the resolution model part of the JSON serialisation
+		/*!
+			\param [in] j the JSON value from which to read the model
+			\param [in] parameters the collection of parameters, some of which should be used to construct the model
+			\return the newly constructed signal model
+		*/
+		static std::unique_ptr<LauDecayTimeResolution> readResolutionModelFromJson( const nlohmann::json& j, std::vector<std::unique_ptr<LauAbsRValue>>& parameters );
+
+		//! Read the background models part of the JSON serialisation
+		/*!
+			\param [in] j the JSON value from which to read the model
+			\param [in] parameters the collection of parameters, which should be used to construct the models
+			\return the newly constructed background models
+		*/
+		static std::map<TString,std::unique_ptr<LauDecayTimeResolution>> readBackgroundModelsFromJson( const nlohmann::json& j, std::vector<std::unique_ptr<LauAbsRValue>>& parameters );
+
 		//! The number of Gaussians in the resolution function
 		const std::size_t nGauss_;
 
 		//! Are any of the parameters of the resolution function scaled by the per-event error?
 		const bool scaleWithPerEventError_{false};
 		//! Scale the mean of each Gaussian by the per-event decay time error?
 		const std::vector<bool> scaleMeans_;
 		//! Scale the width of each Gaussian by the per-event decay time error?
 		const std::vector<bool> scaleWidths_;
+		//! Scaling type
+		const ScalingType scalingType_;
 
 		//! Store of all parameters of the resolution function
+		std::vector<std::unique_ptr<LauAbsRValue>> paramsOwned_;
+		//! Store of all parameters of the resolution function
 		std::vector<LauAbsRValue*> params_;
 
 		//! Fraction parameter for each Gaussian in the resolution function
 		std::vector<LauAbsRValue*> fractions_;
 		//! Mean parameter for each Gaussian in the resolution function
 		std::vector<LauAbsRValue*> means_;
 		//! Width parameter for each Gaussian in the resolution function
 		std::vector<LauAbsRValue*> widths_;
 
 		//! Fraction of each Gaussian in the resolution function
 		std::vector<Double_t> fractionVals_;
 		//! Mean of each Gaussian in the resolution function
 		std::vector<Double_t> meanVals_;
 		//! Width of each Gaussian in the resolution function
 		std::vector<Double_t> widthVals_;
 
 		//! Is any parameter floating in the fit?
 		bool anythingFloating_{false};
 		//! Has any floating parameter changed in the last fit iteration?
 		bool anythingChanged_{false};
 
 		ClassDef(LauDecayTimeResolution, 0)
 };
 
+//! \cond DOXYGEN_IGNORE
+// map LauDecayTimeResolution::ScalingType values to JSON as strings
+NLOHMANN_JSON_SERIALIZE_ENUM( LauDecayTimeResolution::ScalingType, {
+		{LauDecayTimeResolution::ScalingType::None, "None"},
+		{LauDecayTimeResolution::ScalingType::Global, "Global"},
+		{LauDecayTimeResolution::ScalingType::PerGaussian, "PerGaussian"},
+		{LauDecayTimeResolution::ScalingType::PerParameter, "PerParameter"},
+		})
+//! \endcond DOXYGEN_IGNORE
+
+template<typename ResModelPtr>
+void LauDecayTimeResolution::writeToJson( const std::map<TString,ResModelPtr>& models, const TString& fileName, const TString& elementName, const bool append, const int indent )
+{
+	using nlohmann::json;
+
+	// Check for a signal model
+	auto iter { models.find("signal") };
+	if ( iter == models.end() ) {
+		std::cerr << "ERROR in LauDecayTimeResolution::writeToJson : cannot locate the signal model in the collection, aborting" << std::endl;
+		return;
+	}
+
+	json j;
+
+	json& parameters = j["parameters"] = json::array();
+	json& signalModel = j["signalModel"] = json::object();
+	json& backgroundModels = j["backgroundModels"] = json::array();
+
+	for ( const auto& [ name, model ] : models ) {
+		if ( name == "signal" ) {
+			signalModel["scalingType"] = model->getScalingType();
+			signalModel["nGauss"] = model->nGauss();
+			json& sigPars = signalModel["parameters"] = json::array();
+			for ( auto par : model->getParameters() ) {
+				sigPars.push_back( par->name() );
+				parameters.push_back( *par );
+			}
+		} else {
+			json modelDef = json::object();
+			modelDef["name"] = name;
+			json& bkgndModel = modelDef["model"] = json::object();
+			bkgndModel["scalingType"] = model->getScalingType();
+			bkgndModel["nGauss"] = model->nGauss();
+			json& bkgndPars = bkgndModel["parameters"] = json::array();
+			for ( auto par : model->getParameters() ) {
+				bkgndPars.push_back( par->name() );
+				parameters.push_back( *par );
+			}
+			backgroundModels.push_back( modelDef );
+		}
+	}
+
+	const bool writeOK { LauJsonTools::writeJsonFile( j, fileName.Data(), elementName.Data(), append, indent ) };
+	if ( ! writeOK ) {
+		std::cerr << "ERROR in LauDecayTimePhysicsModel::writeToJson : could not successfully write to file \"" << fileName << std::endl;
+	}
+}
 #endif
diff --git a/src/LauDecayTimeResolution.cc b/src/LauDecayTimeResolution.cc
index d0525d7..4ca371d 100644
--- a/src/LauDecayTimeResolution.cc
+++ b/src/LauDecayTimeResolution.cc
@@ -1,209 +1,349 @@
 
 /*
 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 LauDecayTimeResolution.cc
     \brief File containing implementation of LauDecayTimeResolution class.
 */
 
 #include <algorithm>
 #include <iostream>
 #include <numeric>
 #include <vector>
 
 #include "TString.h"
 #include "TSystem.h"
 
 #include "LauAbsRValue.hh"
 #include "LauDecayTimeResolution.hh"
 
 ClassImp(LauDecayTimeResolution);
 
 
-LauDecayTimeResolution::LauDecayTimeResolution( const std::size_t nGauss, const std::vector<LauAbsRValue*>& resolutionParams, const bool scale ) :
+LauDecayTimeResolution::LauDecayTimeResolution( const std::size_t nGauss, std::vector<std::unique_ptr<LauAbsRValue>> resolutionParams, const bool scale ) :
 	nGauss_{nGauss},
 	scaleWithPerEventError_{scale},
 	scaleMeans_(nGauss_,scale),
 	scaleWidths_(nGauss_,scale),
-	params_{resolutionParams}
+	scalingType_{scaleWithPerEventError_ ? ScalingType::Global : ScalingType::None},
+	paramsOwned_{std::move(resolutionParams)}
 {
 	this->checkSetup();
 }
 
-LauDecayTimeResolution::LauDecayTimeResolution( const std::size_t nGauss, const std::vector<LauAbsRValue*>& resolutionParams, const std::vector<bool>& scale ) :
+LauDecayTimeResolution::LauDecayTimeResolution( const std::size_t nGauss, std::vector<std::unique_ptr<LauAbsRValue>> resolutionParams, const std::vector<bool>& scale ) :
 	nGauss_{nGauss},
 	scaleWithPerEventError_( std::accumulate( scale.begin(), scale.end(), false, std::logical_or<bool>() ) ),
 	scaleMeans_{scale},
 	scaleWidths_{scale},
-	params_{resolutionParams}
+	scalingType_{scaleWithPerEventError_ ? ScalingType::PerGaussian : ScalingType::None},
+	paramsOwned_{std::move(resolutionParams)}
 {
 	this->checkSetup();
 }
 
-LauDecayTimeResolution::LauDecayTimeResolution( const std::size_t nGauss, const std::vector<LauAbsRValue*>& resolutionParams, const std::vector<bool>& scaleMeans, const std::vector<bool>& scaleWidths ) :
+LauDecayTimeResolution::LauDecayTimeResolution( const std::size_t nGauss, std::vector<std::unique_ptr<LauAbsRValue>> resolutionParams, const std::vector<bool>& scaleMeans, const std::vector<bool>& scaleWidths ) :
 	nGauss_{nGauss},
-	scaleWithPerEventError_( std::accumulate( scaleMeans.begin(), scaleMeans.end(), false, std::logical_or<bool>() ) || std::accumulate( scaleWidths.begin(), scaleWidths.end(), kFALSE, std::logical_or<bool>() ) ),
+	scaleWithPerEventError_( std::accumulate( scaleMeans.begin(), scaleMeans.end(), false, std::logical_or<bool>() ) || std::accumulate( scaleWidths.begin(), scaleWidths.end(), false, std::logical_or<bool>() ) ),
 	scaleMeans_{scaleMeans},
 	scaleWidths_{scaleWidths},
-	params_{resolutionParams}
+	scalingType_{scaleWithPerEventError_ ? ScalingType::PerParameter : ScalingType::None},
+	paramsOwned_{std::move(resolutionParams)}
 {
 	this->checkSetup();
 }
 
 void LauDecayTimeResolution::resizeVectors()
 {
+	params_.resize( paramsOwned_.size() );
+	std::transform( paramsOwned_.begin(), paramsOwned_.end(), params_.begin(), [](std::unique_ptr<LauAbsRValue>& par){return par.get();} );
+
 	fractions_.assign( nGauss_-1, nullptr );
 	means_.assign( nGauss_, nullptr );
 	widths_.assign( nGauss_, nullptr );
 
 	fractionVals_.assign( nGauss_, 0.0 );
 	meanVals_.assign( nGauss_, 0.0 );
 	widthVals_.assign( nGauss_, 0.0 );
 }
 
 void LauDecayTimeResolution::checkSetup()
 {
 	if ( nGauss_ == 0 ) {
 		std::cerr << "ERROR in LauDecayTimeResolution::checkSetup : number of Gaussians is zero!" << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	if ( scaleWithPerEventError_ ) {
 		// if we're scaling by the per-event error, check that the scale vectors are the right length
 		if ( scaleMeans_.size() != nGauss_ || scaleWidths_.size() != nGauss_ ) {
 			std::cerr << "ERROR in LauDecayTimeResolution::checkSetup : number of Gaussians = " << nGauss_ << " but scale vectors are of length " << scaleMeans_.size() << " and " << scaleWidths_.size() << ", cannot continue!" << std::endl;
 			gSystem->Exit(EXIT_FAILURE);
 		}
 	}
 
 	this->resizeVectors();
 
 	const TString meanNameBase{"mean_"};
 	const TString sigmaNameBase{"sigma_"};
 	const TString fracNameBase{"frac_"};
 
 	bool foundParams{true};
 
 	for (std::size_t i{0}; i < nGauss_; ++i) {
 		TString tempName{meanNameBase}; tempName += i;
 		TString tempName2{sigmaNameBase}; tempName2 += i;
 		TString tempName3{fracNameBase}; tempName3 += i;
 
 		means_[i] = this->findParameter(tempName);
 		foundParams &= (means_[i] != nullptr);
 
 		widths_[i] = this->findParameter(tempName2);
 		foundParams &= (widths_[i] != nullptr);
 
 		if (i!=0) {
 			fractions_[i-1] = this->findParameter(tempName3);
 			foundParams &= (fractions_[i-1] != nullptr);
 		}
 	}
 
 	if ( ! foundParams ) {
 		std::cerr << "ERROR in LauDecayTimeResolution::checkSetup : decay time resolution function with " << nGauss_ << " Gaussians requires:\n";
 		std::cerr << "                                              - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\"\n";
 		std::cerr << "                                              - " << nGauss_-1 << " fractions: \"frac_i\", where i!=0" << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 }
 
 LauAbsRValue* LauDecayTimeResolution::findParameter(const TString& parName)
 {
 	for ( LauAbsRValue* par : params_ ) {
 		if ( par && par->name().Contains(parName) ) {
 			return par;
 		}
 	}
 	std::cerr << "ERROR in LauDecayTimeResolution::findParameter : Parameter \"" << parName << "\" not found." << std::endl;
 	return nullptr;
 }
 
 const LauAbsRValue* LauDecayTimeResolution::findParameter(const TString& parName) const
 {
 	for ( const LauAbsRValue* par : params_ ) {
 		if ( par && par->name().Contains(parName) ) {
 			return par;
 		}
 	}
 	std::cerr << "ERROR in LauDecayTimeResolution::findParameter : Parameter \"" << parName << "\" not found." << std::endl;
 	return nullptr;
 }
 
 void LauDecayTimeResolution::initialise()
 {
 	anythingFloating_ = false;
 	for ( std::size_t i{0}; i < nGauss_; ++i )
 	{
 		const bool meanFloating { not means_[i]->fixed() };
 		const bool widthFloating { not widths_[i]->fixed() };
 
 		anythingFloating_ |= (meanFloating or widthFloating);
 
 		std::cout << "INFO in LauDecayTimeResolution::initialise : mean[" << i << "] floating set to:  " << (meanFloating  ? "True" : "False") << std::endl;
 		std::cout << "INFO in LauDecayTimeResolution::initialise : width[" << i << "] floating set to: " << (widthFloating ? "True" : "False") << std::endl;
 
 		if ( i < (nGauss_ - 1) ) {
 			const bool fracFloating { not fractions_[i]->fixed() };
 
 			anythingFloating_ |= fracFloating;
 
 			std::cout << "INFO in LauDecayTimeResolution::initialise : fraction[" << i << "] floating set to:  " << (fracFloating ? "True" : "False") << std::endl;
 		}
 	}
 
 	this->updateParameterCache();
 }
 
 void LauDecayTimeResolution::updateParameterCache()
 {
 	static auto assignValue = [](const LauAbsRValue* par){return par->unblindValue();};
 
 	std::transform( means_.begin(), means_.end(), meanVals_.begin(), assignValue );
 	std::transform( widths_.begin(), widths_.end(), widthVals_.begin(), assignValue );
 	std::transform( fractions_.begin(), fractions_.end(), fractionVals_.begin()+1, assignValue );
 	fractionVals_[0] = std::accumulate( fractionVals_.begin()+1, fractionVals_.end(), 1.0, std::minus<Double_t>{} );
 }
 
 void LauDecayTimeResolution::propagateParUpdates()
 {
 	// If none of the parameters are floating there's nothing to do
 	if ( not anythingFloating_ ) {
 		return;
 	}
 
 	// Otherwise, determine whether any of the floating parameters have changed and act accordingly
 
 	static auto checkEquality = [](const LauAbsRValue* par, const Double_t cacheVal){return par->unblindValue() == cacheVal;};
 
 	anythingChanged_ = false;
 	anythingChanged_ |= not std::equal( means_.begin(), means_.end(), meanVals_.begin(), checkEquality );
 	anythingChanged_ |= not std::equal( widths_.begin(), widths_.end(), widthVals_.begin(), checkEquality );
 	anythingChanged_ |= not std::equal( fractions_.begin(), fractions_.end(), fractionVals_.begin()+1, checkEquality );
 
 	if ( anythingChanged_ ) {
 		this->updateParameterCache();
 	}
 }
 
+std::map<TString,std::unique_ptr<LauDecayTimeResolution>> LauDecayTimeResolution::readFromJson( const TString& fileName, const TString& elementName)
+{
+	// 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(), elementName.Data(), LauJsonTools::JsonType::Object );
+
+	if ( j.is_null() ) {
+		if ( elementName != "" ) {
+			std::cerr << "ERROR in LauDecayTimeResolution::readFromJson : unable to retrieve JSON object from element \"" << elementName << "\" in file \"" << fileName << "\"" << std::endl;
+		} else {
+			std::cerr << "ERROR in LauDecayTimeResolution::readFromJson : unable to retrieve JSON object from root element of file \"" << fileName << "\"" << std::endl;
+		}
+		return {};
+	}
+
+	return readFromJson( j );
+}
+
+std::map<TString,std::unique_ptr<LauDecayTimeResolution>> LauDecayTimeResolution::readFromJson( const nlohmann::json& j )
+{
+	using JsonType = LauJsonTools::JsonType;
+
+	// the object should have the following elements:
+	// - an array of parameter objects
+	// - an object for the signal model settings
+	// - an array of background model settings
+	const std::vector<LauJsonTools::ElementNameType> mandatoryElements {
+		std::make_pair("parameters", JsonType::Array),
+		std::make_pair("signalModel", JsonType::Object),
+		std::make_pair("backgroundModels", JsonType::Array)
+	};
+	if ( ! LauJsonTools::checkObjectElements( j, mandatoryElements ) ) {
+		std::cerr << "ERROR in LauDecayTimeResolution::readFromJson : JSON object does not contain required elements: \"parameters\" (array), \"signalModel\" (object) and \"backgroundModels\" (array)" << std::endl;
+		return {};
+	}
+
+	auto parameters { LauAbsRValue::readFromJson( j.at("parameters") ) };
+
+	auto sigModel { readResolutionModelFromJson( j.at("signalModel"), parameters ) };
+	if ( ! sigModel ) {
+		return {};
+	}
+
+	auto backgroundModels { readBackgroundModelsFromJson( j.at("backgroundModels"), parameters ) };
+	if ( backgroundModels.size() != j.at("backgroundModels").size() ) {
+		return {};
+	}
+
+	backgroundModels.insert( std::make_pair( "signal", std::move(sigModel) ) );
+
+	return backgroundModels;
+}
+
+std::unique_ptr<LauDecayTimeResolution> LauDecayTimeResolution::readResolutionModelFromJson( const nlohmann::json& j, std::vector<std::unique_ptr<LauAbsRValue>>& parameters )
+{
+	using JsonType = LauJsonTools::JsonType;
+
+	// the object should have the following elements:
+	// - an unsigned integer specifying the number of Gaussians in the model
+	// - a string representing the LauDecayTimeResolution::ScalingType function type
+	// - an array of strings specifying the parameters we need
+	const std::vector<LauJsonTools::ElementNameType> mandatoryElements {
+		std::make_pair("nGauss", JsonType::Number_Unsigned),
+		std::make_pair("scalingType", JsonType::String),
+		std::make_pair("parameters", JsonType::Array)
+	};
+	if ( ! LauJsonTools::checkObjectElements( j, mandatoryElements ) ) {
+		std::cerr << "ERROR in LauDecayTimeResolution::readResolutionModelFromJson : JSON object from \"signalModel\" element does not contain required elements: \"nGauss\" (unsigned), \"scalingType\" (string) and \"parameters\" (array)" << std::endl;
+		return {};
+	}
+
+	auto nGauss { LauJsonTools::getValue<std::size_t>( j, "nGauss" ) };
+	auto scalingType { LauJsonTools::getValue<ScalingType>( j, "scalingType" ) };
+	auto parNames { LauJsonTools::getValue<std::vector<TString>>( j, "parameters" ) };
+
+	std::vector<std::unique_ptr<LauAbsRValue>> params;
+	params.reserve( parNames.size() );
+
+	for ( const auto& name : parNames ) {
+
+		auto iter { std::find_if( parameters.begin(), parameters.end(), [&name]( std::unique_ptr<LauAbsRValue>& par ){ return par->name() == name; } ) };
+		if ( iter == parameters.end() ) {
+			const TString errMsg {"ERROR in LauDecayTimeResolution::readResolutionModelFromJson : Cannot locate parameter \"" + name + "\""};
+			throw LauJsonTools::InvalidJson{errMsg.Data()};
+		}
+
+		params.emplace_back( std::move(*iter) );
+		parameters.erase(iter);
+	}
+
+	switch ( scalingType ) {
+		case ScalingType::None :
+			return std::make_unique<LauDecayTimeResolution>( nGauss, std::move(params), false );
+		case ScalingType::Global :
+			return std::make_unique<LauDecayTimeResolution>( nGauss, std::move(params), true );
+		case ScalingType::PerGaussian :
+			return std::make_unique<LauDecayTimeResolution>( nGauss, std::move(params), LauJsonTools::getValue<std::vector<bool>>( j, "scale" ) );
+		case ScalingType::PerParameter :
+			return std::make_unique<LauDecayTimeResolution>( nGauss, std::move(params), LauJsonTools::getValue<std::vector<bool>>( j, "scaleMeans" ), LauJsonTools::getValue<std::vector<bool>>( j, "scaleWidths" ) );
+	}
+
+	// we shouldn't ever reach here
+	return {};
+}
+
+std::map<TString,std::unique_ptr<LauDecayTimeResolution>> LauDecayTimeResolution::readBackgroundModelsFromJson( const nlohmann::json& j, std::vector<std::unique_ptr<LauAbsRValue>>& parameters )
+{
+	using LauJsonTools::JsonType;
+	using LauJsonTools::checkObjectElements;
+	using LauJsonTools::getValue;
+
+	// each object in the array should have the following elements:
+	// - a string for the name of the background component
+	// - an object defining the model
+	const std::vector<LauJsonTools::ElementNameType> mandatoryElements {
+		std::make_pair("name", JsonType::String),
+		std::make_pair("model", JsonType::Object)
+	};
+	for ( const auto& comp : j ) {
+		if ( ! checkObjectElements( comp, mandatoryElements ) ) {
+			std::cerr << "ERROR in LauDecayTimeResolution::readBackgroundModelsFromJson : JSON object in \"backgroundModels\" array does not contain required elements: \"name\" (string) and \"model\" (object)" << std::endl;
+			return {};
+		}
+	}
+
+	std::map<TString,std::unique_ptr<LauDecayTimeResolution>> models;
+
+	for ( const auto& comp : j ) {
+		auto name { getValue<TString>( comp, "name" ) };
+		models[name] = readResolutionModelFromJson( comp.at("model"), parameters );
+	}
+
+	return models;
+}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index b78b3fb..0dd7395 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,29 +1,30 @@
 
 list(APPEND TEST_SOURCES
     TestCovariant
     TestCovariant2
     TestFitDoubleSplineToTH1
     TestFitSplineToTH1
     TestNewKinematicsMethods
     TestReadCoeffSetFromJson
     TestReadDPModelFromJson
     TestReadDecayTimePhysicsModelFromJson
+    TestReadDecayTimeResolutionModelFromJson
     TestReadParametersFromJson
     TestReadSplineFromJson
     TestReadVetoesFromJson
     TestSplineDTAdivision
     TestSplineFindMax
     TestWriteCoeffSetToJson
     TestWriteParametersToJson
     TestWriteSplineToJson
     TestWriteVetoesToJson
     )
 
 foreach( _test ${TEST_SOURCES})
     add_executable(${_test} ${_test}.cc)
     target_link_libraries(${_test} PRIVATE Laura++)
     install(TARGETS ${_test} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
 endforeach()
 
 install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/test-model.json ${CMAKE_CURRENT_SOURCE_DIR}/test-coeffset.json DESTINATION ${CMAKE_INSTALL_BINDIR})
 
diff --git a/test/TestReadDecayTimeResolutionModelFromJson.cc b/test/TestReadDecayTimeResolutionModelFromJson.cc
new file mode 100644
index 0000000..66ab8f2
--- /dev/null
+++ b/test/TestReadDecayTimeResolutionModelFromJson.cc
@@ -0,0 +1,100 @@
+#include "LauDecayTimeResolution.hh"
+#include "LauParameter.hh"
+
+#include "TString.h"
+
+#include <nlohmann/json.hpp>
+
+#include <fstream>
+#include <iostream>
+#include <memory>
+
+int main()
+{
+	using nlohmann::json;
+
+	for ( UInt_t runNo{1}; runNo <=2; ++runNo ) {
+
+		// signal parameters
+		constexpr std::size_t nGauss{3};
+		constexpr LauDecayTimeResolution::ScalingType scalingType{ LauDecayTimeResolution::ScalingType::None };
+		auto mean0  = std::make_unique<LauParameter>("dt_sig_mean_2",  runNo == 1 ? -0.00174035 : -0.00218543 , -0.01, 0.01, kTRUE );
+		auto mean1  = std::make_unique<LauParameter>("dt_sig_mean_1",  runNo == 1 ?  0.043133   :  0.046862   , -0.10, 0.10, kTRUE );
+		auto mean2  = std::make_unique<LauParameter>("dt_sig_mean_0",  runNo == 1 ? -0.00108159 : -0.00159348 , -0.01, 0.01, kTRUE );
+		auto sigma0 = std::make_unique<LauParameter>("dt_sig_sigma_2", runNo == 1 ?  0.046272   :  0.046064   ,  0.0,  2.0,  kTRUE );
+		auto sigma1 = std::make_unique<LauParameter>("dt_sig_sigma_1", runNo == 1 ?  0.11815    :  0.11870    ,  0.0,  2.5,  kTRUE );
+		auto sigma2 = std::make_unique<LauParameter>("dt_sig_sigma_0", runNo == 1 ?  0.025127   :  0.025077   ,  0.0,  2.5,  kTRUE );
+		auto frac0  = std::make_unique<LauParameter>("dt_sig_frac_2",  runNo == 1 ?  0.44359    :  0.43327    ,  0.0,  1.0,  kTRUE );
+		auto frac1  = std::make_unique<LauParameter>("dt_sig_frac_1",  runNo == 1 ?  0.045545   :  0.045310   ,  0.0,  1.0,  kTRUE );
+
+		std::vector<std::unique_ptr<LauParameter>> sigParams;
+		sigParams.emplace_back( std::move(mean0) );
+		sigParams.emplace_back( std::move(mean1) );
+		sigParams.emplace_back( std::move(mean2) );
+		sigParams.emplace_back( std::move(sigma0) );
+		sigParams.emplace_back( std::move(sigma1) );
+		sigParams.emplace_back( std::move(sigma2) );
+		sigParams.emplace_back( std::move(frac0) );
+		sigParams.emplace_back( std::move(frac1) );
+
+		constexpr std::size_t nBkgnd{3};
+
+		const std::array<TString,nBkgnd> bkgndNames { "Bd2DKpi", "Bs2DKpi", "Lb2Dppi" };
+		std::array<std::vector<std::unique_ptr<LauAbsRValue>>,nBkgnd+1> allParams;
+
+		for ( std::size_t i{0}; i < nBkgnd+1; ++i ) {
+			allParams[i].reserve( sigParams.size() );
+			for ( auto& par : sigParams ) {
+				if ( i < nBkgnd ) {
+					TString parName{ par->name() };
+					parName.ReplaceAll( "_sig_", "_"+bkgndNames[i]+"_" );
+					allParams[i].emplace_back( par->createClone( parName ) );
+				} else {
+					allParams[i].emplace_back( std::move(par) );
+				}
+			}
+		}
+
+		json j = json::object();
+
+		json& parameters = j["parameters"] = json::array();
+		json& signalModel = j["signalModel"] = json::object();
+		json& backgroundModels = j["backgroundModels"] = json::array();
+
+		for ( std::size_t i{0}; i < nBkgnd+1; ++i ) {
+			for ( auto& par : allParams[i] ) {
+				parameters.push_back( *par );
+			}
+		}
+
+		signalModel["nGauss"] = nGauss;
+		signalModel["scalingType"] = scalingType;
+		signalModel["parameters"] = json::array();
+		for ( auto& par : allParams[nBkgnd] ) {
+			signalModel.at("parameters").push_back( par->name() );
+		}
+
+		for ( std::size_t i{0}; i < nBkgnd; ++i ) {
+			backgroundModels.push_back( json::object( { {"name", bkgndNames[i]}, {"model", json::object()} } ) );
+			json& model = backgroundModels.back().at("model");
+			model["nGauss"] = nGauss;
+			model["scalingType"] = scalingType;
+			json& pars = model["parameters"] = json::array();
+			for ( auto& par : allParams[i] ) {
+				pars.push_back( par->name() );
+			}
+		}
+
+		{
+			std::ofstream fout(Form("dt-resolution-Run%d-test1.json",runNo));
+			fout << j.dump(4) << std::endl;
+		}
+
+		auto models = LauDecayTimeResolution::readFromJson( j );
+		for ( auto& [ name, model ] : models ) {
+			std::cout << "Loaded LauDecayTimeResolution for component: " << name << " with nGauss = " << model->nGauss() << std::endl;
+		}
+
+		LauDecayTimeResolution::writeToJson( models, Form("dt-resolution-Run%d-test2.json",runNo) );
+	}
+}