diff --git a/examples/GenFitTimeDep_Bd2Dpipi.cc b/examples/GenFitTimeDep_Bd2Dpipi.cc
index 3494eac..2528429 100644
--- a/examples/GenFitTimeDep_Bd2Dpipi.cc
+++ b/examples/GenFitTimeDep_Bd2Dpipi.cc
@@ -1,1082 +1,1054 @@
 
 #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.200103);
-	//const Double_t maxDt(14.9998);
-	const Double_t minDt(0.2);
-	const Double_t maxDt(15.0);
-	const Double_t minDtErr(0.0);
-	const Double_t maxDtErr(0.15);
-
-	LauParameter * tauB0    = new LauParameter("dt_tau",    1.519,  0.5, 5.0, settings.fixLifetime);
-	LauParameter * deltaMd  = new LauParameter("dt_deltaM", 0.5065, 0.0, 1.0, settings.fixDeltaM);
-
-	std::vector<LauAbsRValue*> dtPhysPars {
-		tauB0,
-		deltaMd
-	};
+	const Double_t minDt{0.2};
+	const Double_t maxDt{15.0};
+	const Double_t minDtErr{0.0};
+	const Double_t maxDtErr{0.15};
+
 
-	auto dtPhysModel = std::make_unique<LauDecayTimePhysicsModel>( LauDecayTime::FuncType::ExpTrig, dtPhysPars );
+	// 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") );
+	if ( lifetimePar ) {
+		lifetimePar->fixed( settings.fixLifetime );
+	}
+	auto deltaMPar = dynamic_cast<LauParameter*>( dtPhysModel->findParameter("dt_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
 	};
 
 	auto dtResoModel = std::make_unique<LauDecayTimeResolution>( nGauss, dtResoPars, scale );
 
-	//Bs backgrounds
-	LauParameter * tauBs    = new LauParameter("dt_tau_Bs", 1.527, 1.526, 1.528, kTRUE);
-	LauParameter * deltaMs  = new LauParameter("dt_deltaM_Bs", 17.765, 17.60, 17.80, kTRUE);
-	LauParameter * deltaGs  = new LauParameter("dt_deltaGamma_Bs", 0.084, 0.080, 0.090, kTRUE);
-
-	//Lb backgrounds
-	LauParameter * tauLb    = new LauParameter("dt_tau_Lb", 1.464, 1.450, 1.470, kTRUE);
-
 	// 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 ) {
 			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 {
 
-			std::vector<LauAbsRValue*> bgPhysParams;
-			if (bg == "Bd2DKpi") {
-				bgPhysParams =
-				{
-					tauB0->createClone(),
-					deltaMd->createClone()
-				};
-			}
-			else if (bg == "Bs2DKpi") {
-				bgPhysParams =
-				{
-					tauBs,
-					deltaGs,
-					deltaMs
-				};
-			}
-			else if (bg == "Lb2Dppi") {
-				bgPhysParams =
-				{
-					tauLb
-				};
-			}
-			auto bgPhysModel = std::make_unique<LauDecayTimePhysicsModel>( bkgndDecayTimeTypes[bg], bgPhysParams );
+			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 );
 
 			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 ) {
 					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/examples/runTimeDepTest.sh b/examples/runTimeDepTest.sh
index 383690e..6749f0c 100755
--- a/examples/runTimeDepTest.sh
+++ b/examples/runTimeDepTest.sh
@@ -1,248 +1,254 @@
 #!/bin/bash
 #SBATCH --partition=epp
 #SBATCH --ntasks=1
 #SBATCH --cpus-per-task=6
 #SBATCH --mem-per-cpu=3997
 #SBATCH --time=02:00:00
 #Run with sbatch runTimeDepTest.sh
 
 # Copyright 2020 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
 
 processes=()
 
 # trap ^C and call ctrl_c()
 trap ctrl_c INT
 
 function ctrl_c()
 {
 	echo " ** Caught ^C event!"
 	for process in ${processes[@]}
 	do
 		if ps -u $USER | grep -q "$process"
 		then
 			echo "Killing : $process"
 			kill $process
 		fi
 	done
 	exit 1
 }
 
 # Check that we have a kerberos token, since nothing will work without one!
 klist -s > /dev/null 2>&1
 if [ $? -ne 0 ]
 then
     echo "You don't have a valid Kerberos ticket!"
     exit 1
 fi
 
 iFit=0
 dta_model="spline"
 dtr=1
 dtr_perevent=0
 fixTau=1
 fixDeltaM=1
 seed=140279 # $(pwd | xargs basename)
 dataFit=1
 runNo=12
 
 # Modify these if you want to run more experiments
 nExpt=1
 firstExpt=0
 
 logName_genQFS_run1=Logs/gen-QFS-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
 logName_genQFS_run2=Logs/gen-QFS-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
 logName_genCPE_KK_run1=Logs/gen-CPEvenKK-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
 logName_genCPE_KK_run2=Logs/gen-CPEvenKK-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
 logName_genCPE_pipi_run1=Logs/gen-CPEvenpipi-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
 logName_genCPE_pipi_run2=Logs/gen-CPEvenpipi-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
 
 logName_coord=Logs/coordinator-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
 logName_taskQFS_run1=Logs/task-QFS-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
 logName_taskQFS_run2=Logs/task-QFS-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
 logName_taskCPE_KK_run1=Logs/task-CPEvenKK-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
 logName_taskCPE_KK_run2=Logs/task-CPEvenKK-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
 logName_taskCPE_pipi_run1=Logs/task-CPEvenpipi-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
 logName_taskCPE_pipi_run2=Logs/task-CPEvenpipi-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
 
 mkdir -p Logs
 
 # make sure we have the most up-to-date json files
 eosArea="root://eoslhcb.cern.ch//eos/lhcb/wg/b2oc/B2D0pipi/timedep_fit_histos"
-rm -f run1DTAsplines.json run2DTAsplines.json Bd2D0pipi_DP_Model_Coeffs.json Bd2D0pipi_DP_Vetoes.json
+rm -f run1DTAsplines.json run2DTAsplines.json Bd2D0pipi_DP_Model_Coeffs.json Bd2D0pipi_DP_Vetoes.json Bd2D0pipi_DecayTime.json
 xrdcp $eosArea/run1DTAsplines.json .
 if [ $? -ne 0 ]
 then
     echo "Problem downloading json files from EOS, exiting"
     exit 1
 fi
 xrdcp $eosArea/run2DTAsplines.json .
 if [ $? -ne 0 ]
 then
     echo "Problem downloading json files from EOS, exiting"
     exit 1
 fi
 xrdcp $eosArea/Bd2D0pipi_DP_Model_Coeffs.json .
 if [ $? -ne 0 ]
 then
     echo "Problem downloading json files from EOS, exiting"
     exit 1
 fi
 xrdcp $eosArea/Bd2D0pipi_DP_Vetoes.json .
 if [ $? -ne 0 ]
 then
     echo "Problem downloading json files from EOS, exiting"
     exit 1
 fi
+xrdcp $eosArea/Bd2D0pipi_DecayTime.json .
+if [ $? -ne 0 ]
+then
+    echo "Problem downloading json files from EOS, exiting"
+    exit 1
+fi
 
 if [ $iFit == 0 -a $dataFit == 0 ]
 then
     #run generation of the modes in parallel (put the tasks in the bg)
     echo "Generating samples..."
     if [ $runNo -eq 1 -o $runNo -eq 12 ]
     then
         GenFitTimeDep_Bd2Dpipi gen --dtype QFS    --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_Kpi > $logName_genQFS_run1 2>&1 &
         task=$!
         processes+=( $task )
 
         GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_KK > $logName_genCPE_KK_run1 2>&1 &
         task=$!
         processes+=( $task )
 
         GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_pipi > $logName_genCPE_pipi_run1 2>&1 &
         task=$!
         processes+=( $task )
     fi
 
     if [ $runNo -eq 2 -o $runNo -eq 12 ]
     then
         GenFitTimeDep_Bd2Dpipi gen --dtype QFS    --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_Kpi > $logName_genQFS_run2 2>&1 &
         task=$!
         processes+=( $task )
 
         GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_KK > $logName_genCPE_KK_run2 2>&1 &
         task=$!
         processes+=( $task )
 
         GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_pipi > $logName_genCPE_pipi_run2 2>&1 &
         task=$!
         processes+=( $task )
     fi
 
     wait #wait for the generation to complete
 fi
 
 echo "Running fit $iFit"
 if [ $runNo -eq 12 ]
 then
     nTasks=6
 else
     nTasks=3
 fi
 
 SimFitCoordinator $iFit $nExpt $firstExpt $nTasks > $logName_coord 2>&1 &
 task=$!
 echo "SimFitCoordinator process : $task"
 processes+=( $task )
 
 echo "Initialised coordinator process for fit $iFit"
 
 sleep 5
 
 port=""
 while [ "x$port" == "x" ]
 do
     sleep 2
     port=`tail $logName_coord | grep "Waiting for connection" | awk '{print $NF}'`
 done
 
 echo "Coordinator is listening on port $port"
 echo "Initialising tasks..."
 
 if [ $runNo -eq 1 -o $runNo -eq 12 ]
 then
     GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1  --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run1 2>&1 &
     task=$!
     echo "QFS Run 1 process : $task"
     processes+=( $task )
     sleep 5
     GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1  --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run1 2>&1 &
     task=$!
     echo "CP KK Run 1 process : $task"
     processes+=( $task )
     sleep 5
     GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1  --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run1 2>&1 &
     task=$!
     echo "CP pipi Run 1 process : $task"
     processes+=( $task )
 fi
 if [ $runNo -eq 2 -o $runNo -eq 12 ]
 then
     GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2  --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run2 2>&1 &
     task=$!
     echo "QFS Run 2 process : $task"
     processes+=( $task )
     sleep 5
     GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2  --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run2 2>&1 &
     task=$!
     echo "CP KK Run 2 process : $task"
     processes+=( $task )
     sleep 5
     GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2  --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run2 2>&1 &
     task=$!
     echo "CP pipi Run 2 process : $task"
     processes+=( $task )
 fi
 
 wait #wait for the fits and coordinator to terminate
 
 echo "Fit completed"
 
 echo "Warnings/Errors from coordinator:"
 grep -e ERROR -e WARNING -e Error -e Warning $logName_coord
 
 echo "Warnings/Errors from CPeven KK tasks:"
 if [ $runNo -eq 1 -o $runNo -eq 12 ]
 then
     grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_KK_run1
 fi
 if [ $runNo -eq 2 -o $runNo -eq 12 ]
 then
     grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_KK_run2
 fi
 
 echo "Warnings/Errors from CPeven pipi tasks:"
 if [ $runNo -eq 1 -o $runNo -eq 12 ]
 then
     grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_pipi_run1
 fi
 if [ $runNo -eq 2 -o $runNo -eq 12 ]
 then
     grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_pipi_run2
 fi
 
 echo "Warnings/Errors from QFS tasks:"
 if [ $runNo -eq 1 -o $runNo -eq 12 ]
 then
     grep -e ERROR -e WARNING -e Error -e Warning $logName_taskQFS_run1
 fi
 if [ $runNo -eq 2 -o $runNo -eq 12 ]
 then
     grep -e ERROR -e WARNING -e Error -e Warning $logName_taskQFS_run2
 fi
diff --git a/inc/LauDecayTime.hh b/inc/LauDecayTime.hh
index ef1d1c2..7fe4ac4 100644
--- a/inc/LauDecayTime.hh
+++ b/inc/LauDecayTime.hh
@@ -1,63 +1,88 @@
 
 /*
 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 LauDecayTime.hh
     \brief File containing declaration of utilities for decay time modelling
 */
 
 #ifndef LAU_DECAYTIME_UTILS
 #define LAU_DECAYTIME_UTILS
 
+#include <nlohmann/json.hpp>
+
 /*! \namespace LauDecayTime
     \brief Namespace to contain various definitions that are common to the description of the decay time dependence of the amplitudes
 */
 
 namespace LauDecayTime {
 
 	// TODO - can we think of better names?
 	//! The functional form of the decay time PDF
 	enum class FuncType {
 		Hist,		//< Hist PDF for fixed background
 		Delta,		//< Delta function - for prompt background
 		Exp,		//< Exponential function - for non-prompt background or charged B's
 		DeltaExp,	//< Delta + Exponential function - for background with prompt and non-prompt parts
 		ExpTrig,	//< Exponential function with Delta m driven mixing - for neutral B_d's
 		ExpHypTrig	//< Exponential function with both Delta m and Delta Gamma driven mixing - for neutral B_s's
 	};
 
 	//! How is the decay time measured - absolute or difference?
 	enum class TimeMeasurementMethod {
 		DecayTime,	//< Absolute measurement of decay time, e.g. LHCb scenario
 		DecayTimeDiff	//< Measurement of the difference of two decay times, e.g. BaBar/Belle(II) scenario
 	};
 
 	//! How is the TD efficiency information going to be given?
 	enum class EfficiencyMethod {
 		Uniform,	//< As a uniform distribution (constant)
 		Binned,		//< As a histogram (TH1D/TH1F)
 		Spline		//< As a cubic spline (or products thereof)
 	};
+
+	//! \cond DOXYGEN_IGNORE
+	// map FuncType values to JSON as strings
+	NLOHMANN_JSON_SERIALIZE_ENUM( FuncType, {
+			{FuncType::Hist, "Hist"},
+			{FuncType::Delta, "Delta"},
+			{FuncType::Exp, "Exp"},
+			{FuncType::DeltaExp, "DeltaExp"},
+			{FuncType::ExpTrig, "ExpTrig"},
+			{FuncType::ExpHypTrig, "ExpHypTrig"},
+			})
+	// map TimeMeasurementMethod values to JSON as strings
+	NLOHMANN_JSON_SERIALIZE_ENUM( TimeMeasurementMethod, {
+			{TimeMeasurementMethod::DecayTime, "DecayTime"},
+			{TimeMeasurementMethod::DecayTimeDiff, "DecayTimeDiff"},
+			})
+	// map EfficiencyMethod values to JSON as strings
+	NLOHMANN_JSON_SERIALIZE_ENUM( EfficiencyMethod, {
+			{EfficiencyMethod::Uniform, "Uniform"},
+			{EfficiencyMethod::Binned, "Binned"},
+			{EfficiencyMethod::Spline, "Spline"},
+			})
+	//! \endcond DOXYGEN_IGNORE
 }
 
 #endif
diff --git a/inc/LauDecayTimePhysicsModel.hh b/inc/LauDecayTimePhysicsModel.hh
index a496316..b6354fd 100644
--- a/inc/LauDecayTimePhysicsModel.hh
+++ b/inc/LauDecayTimePhysicsModel.hh
@@ -1,220 +1,270 @@
 
 /*
 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 LauDecayTimePhysicsModel.hh
     \brief File containing declaration of LauDecayTimePhysicsModel class.
 */
 
 #ifndef LAU_DECAYTIME_PHYSICSMODEL
 #define LAU_DECAYTIME_PHYSICSMODEL
 
-#include <vector>
+#include "LauDecayTime.hh"
 
 #include "Rtypes.h"
+#include "TString.h"
 
-#include "LauDecayTime.hh"
+#include <nlohmann/json_fwd.hpp>
+
+#include <map>
+#include <memory>
+#include <vector>
 
 class LauAbsRValue;
 
 /*! \class LauDecayTimePhysicsModel
     \brief Class for defining the physics and associated parameters for the decay time model
 */
 
 class LauDecayTimePhysicsModel final {
 
 	public:
 		//! Constructor
 		/*!
 		    \param type the type of the physics model
 		    \param parameters the parameters of the physics model
 		*/
-		LauDecayTimePhysicsModel( const LauDecayTime::FuncType type, std::vector<LauAbsRValue*> parameters );
+		LauDecayTimePhysicsModel( const LauDecayTime::FuncType type, std::vector<std::unique_ptr<LauAbsRValue>> parameters );
 
 		//! Retrieve the type of the physics function
 		/*!
 		    \return the function type
 		*/
 		LauDecayTime::FuncType getFunctionType() const { return type_; }
 
 		//! Retrieve the parameters of the physics model so that they can be loaded into a fit
 		/*!
 		    \return the parameters of the physics model
 		*/
 		std::vector<LauAbsRValue*> getParameters() { return params_; }
 
 		//! Retrieve the up-to-date value of the lifetime
 		/*!
 		    \return the lifetime value
 		*/
 		const Double_t& getLifetimeValue() const { return tauVal_; }
 
 		//! Retrieve the up-to-date value of the inverse lifetime
 		/*!
 		    \return the inverse lifetime value
 		*/
 		const Double_t& getGammaValue() const { return gammaVal_; }
 
 		//! Retrieve the up-to-date value of the mass difference
 		/*!
 		    \return the mass difference value
 		*/
 		const Double_t& getDeltaMValue() const { return deltaMVal_; }
 
 		//! Retrieve the up-to-date value of the width difference
 		/*!
 		    \return the width difference value
 		*/
 		const Double_t& getDeltaGammaValue() const { return deltaGammaVal_; }
 
 		//! Retrieve the up-to-date value of the prompt fraction
 		/*!
 		    \return the prompt fraction value
 		*/
 		const Double_t& getFracPromptValue() const { return fracPromptVal_; }
 
 		//! 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_; }
 
 		//! Retrieve whether the lifetime value has changed in the last fit iteration
 		/*!
 		    \return the lifetime changed flag
 		*/
 		const bool& lifetimeChanged() const { return tauChanged_; }
 
 		//! Retrieve the up-to-date value of the mass difference
 		/*!
 		    \return the mass difference value
 		*/
 		const bool& deltaMChanged() const { return deltaMChanged_; }
 
 		//! Retrieve the up-to-date value of the width difference
 		/*!
 		    \return the width difference value
 		*/
 		const bool& deltaGammaChanged() const { return deltaGammaChanged_; }
 
 		//! Retrieve the up-to-date value of the prompt fraction
 		/*!
 		    \return the prompt fraction value
 		*/
 		const bool& fracPromptChanged() const { return fracPromptChanged_; }
 
 		//! 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 physics 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:
+			  + "functionType", a String to specify the function type (must match a state of LauDecayTime::FuncType)
+			  + "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<LauDecayTimePhysicsModel>> readFromJson( const nlohmann::json& j );
+
+		//! Construct a collection of physics 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<LauDecayTimePhysicsModel>> readFromJson( const TString& fileName, const TString& elementName = "" );
+
 		//! 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:
 		//! Update the cached parameter values
 		void updateParameterCache();
 
+		//! Read the signal 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<LauDecayTimePhysicsModel> readSignalModelFromJson( 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<LauDecayTimePhysicsModel>> readBackgroundModelsFromJson( const nlohmann::json& j, std::vector<std::unique_ptr<LauAbsRValue>>& parameters );
+
 		//! Which type of decay time function is this?
 		const LauDecayTime::FuncType type_;
 
 		//! Store of all parameters of the physics model
+		std::vector<std::unique_ptr<LauAbsRValue>> paramsOwned_;
+
+		//! Store of all parameters of the physics model
 		std::vector<LauAbsRValue*> params_;
 
 		//! Lifetime parameter
 		LauAbsRValue* tau_{nullptr};
 
 		//! Mass difference parameter
 		LauAbsRValue* deltaM_{nullptr};
 
 		//! Width difference parameter
 		LauAbsRValue* deltaGamma_{nullptr};
 
 		//! Parameter for the fraction of prompt events in DeltaExp
 		LauAbsRValue* fracPrompt_{nullptr};
 
 		//! Cached value of lifetime
 		Double_t tauVal_{0.0};
 
 		//! Cached value of 1/lifetime
 		Double_t gammaVal_{0.0};
 
 		//! Cached value of mass difference
 		Double_t deltaMVal_{0.0};
 
 		//! Cached value of width difference
 		Double_t deltaGammaVal_{0.0};
 
 		//! Cached value of prompt fraction
 		Double_t fracPromptVal_{0.0};
 
 		//! Is the lifetime floating in the fit?
 		bool tauFloating_{false};
 
 		//! Is the mass difference floating in the fit?
 		bool deltaMFloating_{false};
 
 		//! Is the width difference floating in the fit?
 		bool deltaGammaFloating_{false};
 
 		//! Is the prompt fraction floating in the fit?
 		bool fracPromptFloating_{false};
 
 		//! Are any of the physics parameters floating in the fit?
 		bool anythingFloating_{false};
 
 		//! Has the lifetime parameter changed in the last fit iteration?
 		bool tauChanged_{false};
 
 		//! Has the mass difference parameter changed in the last fit iteration?
 		bool deltaMChanged_{false};
 
 		//! Has the width difference parameter changed in the last fit iteration?
 		bool deltaGammaChanged_{false};
 
 		//! Has the prompt fraction parameter changed in the last fit iteration?
 		bool fracPromptChanged_{false};
 
 		//! Have any of the physics parameters changed in the last fit iteration?
 		bool anythingChanged_{false};
 
 		ClassDef(LauDecayTimePhysicsModel, 0)
 };
 
 #endif
diff --git a/src/LauDecayTimePhysicsModel.cc b/src/LauDecayTimePhysicsModel.cc
index c6db2fe..6c4c664 100644
--- a/src/LauDecayTimePhysicsModel.cc
+++ b/src/LauDecayTimePhysicsModel.cc
@@ -1,196 +1,352 @@
 
 /*
 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 LauDecayTimePhysicsModel.cc
     \brief File containing implementation of LauDecayTimePhysicsModel class.
 */
 
+#include <algorithm>
 #include <iostream>
 #include <vector>
 
 #include "TSystem.h"
 
 #include "LauAbsRValue.hh"
 #include "LauDecayTime.hh"
 #include "LauDecayTimePhysicsModel.hh"
 
 ClassImp(LauDecayTimePhysicsModel);
 
 
-LauDecayTimePhysicsModel::LauDecayTimePhysicsModel( const LauDecayTime::FuncType type, std::vector<LauAbsRValue*> parameters ) :
+LauDecayTimePhysicsModel::LauDecayTimePhysicsModel( const LauDecayTime::FuncType type, std::vector<std::unique_ptr<LauAbsRValue>> parameters ) :
 	type_{type},
-	params_{std::move(parameters)}
+	paramsOwned_{std::move(parameters)}
 {
+	params_.resize( paramsOwned_.size() );
+	std::transform( paramsOwned_.begin(), paramsOwned_.end(), params_.begin(), [](std::unique_ptr<LauAbsRValue>& par){return par.get();} );
+
 	bool foundParams{true};
 	switch ( type_ ) {
 		case LauDecayTime::FuncType::Hist :
 			std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : Hist type needs no physics model" << std::endl;
 			gSystem->Exit(EXIT_FAILURE);
 			break;
 		case LauDecayTime::FuncType::Delta :
 			if ( ! params_.empty() ) {
 				std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : Delta type model requires no parameters" << std::endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 			break;
 		case LauDecayTime::FuncType::Exp :
 			tau_ = this->findParameter("tau");
 			foundParams &= (tau_ != nullptr);
 			if ( params_.size() != 1 || (!foundParams)) {
 				std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : Exp type model requires the following parameters:\n";
 				std::cerr << "                                                              - the lifetime of the exponential decay: \"tau\"" << std::endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 			break;
 		case LauDecayTime::FuncType::DeltaExp :
 			tau_ = this->findParameter("tau");
 			fracPrompt_ = this->findParameter("frac_prompt");
 			foundParams &= (tau_ != nullptr);
 			foundParams &= (fracPrompt_ != nullptr);
 			if ( params_.size() != 2 || (!foundParams)) {
 				std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : DeltaExp type model requires the following parameters:\n";
 				std::cerr << "                                                              - the lifetime of the exponential decay: \"tau\"\n";
 				std::cerr << "                                                              - the fraction of the prompt part: \"frac_prompt\"" << std::endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 			break;
 		case LauDecayTime::FuncType::ExpTrig :
 			tau_ = this->findParameter("tau");
 			deltaM_ = this->findParameter("deltaM");
 			foundParams &= (tau_ != nullptr);
 			foundParams &= (deltaM_ != nullptr);
 			if ( params_.size() != 2 || (!foundParams)) {
 				std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : ExpTrig type model requires the following parameters:\n";
 				std::cerr << "                                                              - the lifetime of the exponential decay: \"tau\"\n";
 				std::cerr << "                                                              - the mass difference: \"deltaM\"" << std::endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 			break;
 		case LauDecayTime::FuncType::ExpHypTrig :
 			tau_        = this->findParameter("tau");
 			deltaM_     = this->findParameter("deltaM");
 			deltaGamma_ = this->findParameter("deltaGamma");
 			foundParams &= (tau_ != nullptr);
 			foundParams &= (deltaM_ != nullptr);
 			foundParams &= (deltaGamma_ != nullptr);
 			if ( params_.size() != 3 || (!foundParams)) {
 				std::cerr << "ERROR in LauDecayTimePhysicsModel::LauDecayTimePhysicsModel : ExpHypTrig type model requires the following parameters:\n";
 				std::cerr << "                                                              - the lifetime of the exponential decay: \"tau\"\n";
 				std::cerr << "                                                              - the mass difference: \"deltaM\"\n";
 				std::cerr << "                                                              - the width difference: \"deltaGamma\"" << std::endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 			break;
 	}
 }
 
 LauAbsRValue* LauDecayTimePhysicsModel::findParameter(const TString& parName)
 {
 	for ( LauAbsRValue* par : params_ ) {
 		if ( par && par->name().Contains(parName) ) {
 			return par;
 		}
 	}
 	std::cerr << "ERROR in LauDecayTimePhysicsModel::findParameter : Parameter \"" << parName << "\" not found." << std::endl;
 	return nullptr;
 }
 
 const LauAbsRValue* LauDecayTimePhysicsModel::findParameter(const TString& parName) const
 {
 	for ( const LauAbsRValue* par : params_ ) {
 		if ( par && par->name().Contains(parName) ) {
 			return par;
 		}
 	}
 	std::cerr << "ERROR in LauDecayTimePhysicsModel::findParameter : Parameter \"" << parName << "\" not found." << std::endl;
 	return nullptr;
 }
 
 void LauDecayTimePhysicsModel::initialise()
 {
 	tauFloating_ = tauChanged_ = ( tau_ != nullptr );
 	deltaMFloating_ = deltaMChanged_ = ( deltaM_ != nullptr );
 	deltaGammaFloating_ = deltaGammaChanged_ = ( deltaGamma_ != nullptr );
 	fracPromptFloating_ = fracPromptChanged_ = ( fracPrompt_ != nullptr );
 
 	this->updateParameterCache();
 
 	tauFloating_ = ( tau_ != nullptr ) and not tau_->fixed();
 	deltaMFloating_ = ( deltaM_ != nullptr ) and not deltaM_->fixed();
 	deltaGammaFloating_ = ( deltaGamma_ != nullptr ) and not deltaGamma_->fixed();
 	fracPromptFloating_ = ( fracPrompt_ != nullptr ) and not fracPrompt_->fixed();
 
 	anythingFloating_ = ( tauFloating_ or deltaMFloating_ or deltaGammaFloating_ or fracPromptFloating_ );
 
 	std::cout << "INFO in LauDecayTimePhysicsModel::initialise : tau floating set to:        " << (tauFloating_        ? "True" : "False") << std::endl;
 	std::cout << "INFO in LauDecayTimePhysicsModel::initialise : deltaM floating set to:     " << (deltaMFloating_     ? "True" : "False") << std::endl;
 	std::cout << "INFO in LauDecayTimePhysicsModel::initialise : deltaGamma floating set to: " << (deltaGammaFloating_ ? "True" : "False") << std::endl;
 	std::cout << "INFO in LauDecayTimePhysicsModel::initialise : fracPrompt floating set to: " << (fracPromptFloating_ ? "True" : "False") << std::endl;
 }
 
 void LauDecayTimePhysicsModel::updateParameterCache()
 {
 	// Get the updated values of all parameters
 
 	if ( tauChanged_ ) {
 		tauVal_ = tau_->unblindValue();
 		gammaVal_ = 1.0 / tauVal_;
 	}
 
 	if ( deltaMChanged_ ) {
 		deltaMVal_ = deltaM_->unblindValue();
 	}
 
 	if ( deltaGammaChanged_ ) {
 		deltaGammaVal_ = deltaGamma_->unblindValue();
 	}
 
 	if ( fracPromptChanged_ ) {
 		fracPromptVal_ = fracPrompt_->unblindValue();
 	}
 }
 
 void LauDecayTimePhysicsModel::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;};
 
 	tauChanged_ = tauFloating_ and not checkEquality(tau_, tauVal_);
 	deltaMChanged_ = deltaMFloating_ and not checkEquality(deltaM_, deltaMVal_);
 	deltaGammaChanged_ = deltaGammaFloating_ and not checkEquality(deltaGamma_, deltaGammaVal_);
 	fracPromptChanged_ = fracPromptFloating_ and not checkEquality(fracPrompt_, fracPromptVal_);
 
 	anythingChanged_ = ( tauChanged_ or deltaMChanged_ or deltaGammaChanged_ or fracPromptChanged_ );
 
 	if ( anythingChanged_ ) {
 		this->updateParameterCache();
 	}
 }
 
+std::map<TString,std::unique_ptr<LauDecayTimePhysicsModel>> LauDecayTimePhysicsModel::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 LauDecayTimePhysicsModel::readFromJson : unable to retrieve JSON object from element \"" << elementName << "\" in file \"" << fileName << "\"" << std::endl;
+		} else {
+			std::cerr << "ERROR in LauDecayTimePhysicsModel::readFromJson : unable to retrieve JSON object from root element of file \"" << fileName << "\"" << std::endl;
+		}
+		return {};
+	}
+
+	return readFromJson( j );
+}
+
+std::map<TString,std::unique_ptr<LauDecayTimePhysicsModel>> LauDecayTimePhysicsModel::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 LauDecayTimePhysicsModel::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 { readSignalModelFromJson( 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<LauDecayTimePhysicsModel> LauDecayTimePhysicsModel::readSignalModelFromJson( const nlohmann::json& j, std::vector<std::unique_ptr<LauAbsRValue>>& parameters )
+{
+	using JsonType = LauJsonTools::JsonType;
+
+	// the object should have the following elements:
+	// - a string representing the LauDecayTime::FuncType function type
+	// - an array of strings specifying the parameters we need
+	const std::vector<LauJsonTools::ElementNameType> mandatoryElements {
+		std::make_pair("functionType", JsonType::String),
+		std::make_pair("parameters", JsonType::Array)
+	};
+	if ( ! LauJsonTools::checkObjectElements( j, mandatoryElements ) ) {
+		std::cerr << "ERROR in LauDecayTimePhysicsModel::readSignalModelFromJson : JSON object from \"signalModel\" element does not contain required elements: \"functionType\" (string) and \"parameters\" (array)" << std::endl;
+		return {};
+	}
+
+	auto funcType { LauJsonTools::getValue<LauDecayTime::FuncType>( j, "functionType" ) };
+	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 LauDecayTimePhysicsModel::readSignalModelFromJson : Cannot locate parameter \"" + name + "\""};
+			throw LauJsonTools::InvalidJson{errMsg.Data()};
+		}
+
+		params.emplace_back( std::move(*iter) );
+		parameters.erase(iter);
+	}
+
+	return std::make_unique<LauDecayTimePhysicsModel>( funcType, std::move(params) );
+}
+
+std::map<TString,std::unique_ptr<LauDecayTimePhysicsModel>> LauDecayTimePhysicsModel::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, which itself has elements:
+	//   = a string representing the LauDecayTime::FuncType function type
+	//   = an array of strings specifying the parameters we need
+	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 LauDecayTimePhysicsModel::readBackgroundModelsFromJson : JSON object in \"backgroundModels\" array does not contain required elements: \"name\" (string) and \"model\" (object)" << std::endl;
+			return {};
+		}
+	}
+
+	mandatoryElements = {
+		std::make_pair("functionType", JsonType::String),
+		std::make_pair("parameters", JsonType::Array)
+	};
+
+	std::map<TString,std::unique_ptr<LauDecayTimePhysicsModel>> models;
+
+	for ( const auto& comp : j ) {
+		auto name { getValue<TString>( comp, "name" ) };
+
+		const auto& model { comp.at("model") };
+		if ( ! checkObjectElements( model, mandatoryElements ) ) {
+			std::cerr << "ERROR in LauDecayTimePhysicsModel::readBackgroundModelsFromJson : JSON object defining model in \"backgroundModels\" array does not contain required elements: \"functionType\" (string) and \"parameters\" (array)" << std::endl;
+			return {};
+		}
+
+		auto funcType { getValue<LauDecayTime::FuncType>( model, "functionType" ) };
+		auto parNames { getValue<std::vector<TString>>( model, "parameters" ) };
+
+		std::vector<std::unique_ptr<LauAbsRValue>> params;
+		params.reserve( parNames.size() );
+
+		for ( const auto& parName : parNames ) {
+
+			auto iter { std::find_if( parameters.begin(), parameters.end(), [&parName]( std::unique_ptr<LauAbsRValue>& par ){ return par->name() == parName; } ) };
+			if ( iter == parameters.end() ) {
+				const TString errMsg {"ERROR in LauDecayTimePhysicsModel::readBackgroundModelsFromJson : Cannot locate parameter \"" + parName + "\""};
+				throw LauJsonTools::InvalidJson{errMsg.Data()};
+			}
+
+			params.emplace_back( std::move(*iter) );
+			parameters.erase(iter);
+		}
+
+		models[name] = std::make_unique<LauDecayTimePhysicsModel>( funcType, std::move(params) );
+	}
+
+	return models;
+}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f7a6dc4..b78b3fb 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,28 +1,29 @@
 
 list(APPEND TEST_SOURCES
     TestCovariant
     TestCovariant2
-    TestNewKinematicsMethods
-    TestFitSplineToTH1
     TestFitDoubleSplineToTH1
-    TestSplineDTAdivision
-    TestWriteCoeffSetToJson
-    TestWriteParametersToJson
-    TestWriteSplineToJson
-    TestWriteVetoesToJson
+    TestFitSplineToTH1
+    TestNewKinematicsMethods
     TestReadCoeffSetFromJson
     TestReadDPModelFromJson
+    TestReadDecayTimePhysicsModelFromJson
     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/TestReadDecayTimePhysicsModelFromJson.cc b/test/TestReadDecayTimePhysicsModelFromJson.cc
new file mode 100644
index 0000000..81563e1
--- /dev/null
+++ b/test/TestReadDecayTimePhysicsModelFromJson.cc
@@ -0,0 +1,63 @@
+#include "LauDecayTimePhysicsModel.hh"
+#include "LauParameter.hh"
+
+#include <nlohmann/json.hpp>
+
+#include <fstream>
+#include <iostream>
+#include <memory>
+
+int main()
+{
+	using nlohmann::json;
+
+	json j = json::object();
+
+	json& parameters = j["parameters"] = json::array();
+	json& signalModel = j["signalModel"] = json::object();
+	json& backgroundModels = j["backgroundModels"] = json::array();
+
+	// signal parameters
+	auto tauBd   = std::make_unique<LauParameter>("dt_sig_tau_Bd",    1.519,  0.5, 5.0, kTRUE);
+	auto deltaMd = std::make_unique<LauParameter>("dt_sig_deltaM_Bd", 0.5065, 0.0, 1.0, kFALSE);
+
+	// Bd2DKpi background parameters
+	auto tauBd_Bkgnd   = std::unique_ptr<LauParameter>(tauBd->createClone("dt_bkgnd_tau_Bd"));
+	auto deltaMd_Bkgnd = std::unique_ptr<LauParameter>(deltaMd->createClone("dt_bkgnd_deltaM_Bd"));
+
+	// Bs2DKpi background parameters
+	auto tauBs   = std::make_unique<LauParameter>("dt_tau_Bs", 1.527, 1.526, 1.528, kTRUE);
+	auto deltaMs = std::make_unique<LauParameter>("dt_deltaM_Bs", 17.765, 17.60, 17.80, kTRUE);
+	auto deltaGs = std::make_unique<LauParameter>("dt_deltaGamma_Bs", 0.084, 0.080, 0.090, kTRUE);
+
+	// Lb2Dppi background parameters
+	auto tauLb = std::make_unique<LauParameter>("dt_tau_Lb", 1.464, 1.450, 1.470, kTRUE);
+
+	std::vector<std::unique_ptr<LauAbsRValue>> params;
+	params.emplace_back( std::move(tauBd) );
+	params.emplace_back( std::move(deltaMd) );
+	params.emplace_back( std::move(tauBd_Bkgnd) );
+	params.emplace_back( std::move(deltaMd_Bkgnd) );
+	params.emplace_back( std::move(tauBs) );
+	params.emplace_back( std::move(deltaMs) );
+	params.emplace_back( std::move(deltaGs) );
+	params.emplace_back( std::move(tauLb) );
+
+	for ( auto& par : params ) {
+		parameters.push_back( *par );
+	}
+
+	signalModel["functionType"] = LauDecayTime::FuncType::ExpTrig;
+	signalModel["parameters"] = { "dt_sig_tau_Bd", "dt_sig_deltaM_Bd" };
+
+	backgroundModels.push_back( json::object( { {"name", "Bd2DKpi"}, {"model", json::object( { {"functionType", LauDecayTime::FuncType::ExpTrig}, {"parameters", json::array({"dt_bkgnd_tau_Bd","dt_bkgnd_deltaM_Bd"})} } )} } ) );
+	backgroundModels.push_back( json::object( { {"name", "Bs2DKpi"}, {"model", json::object( { {"functionType", LauDecayTime::FuncType::ExpHypTrig}, {"parameters", json::array({"dt_tau_Bs","dt_deltaM_Bs","dt_deltaGamma_Bs"})} } )} } ) );
+	backgroundModels.push_back( json::object( { {"name", "Lb2Dppi"}, {"model", json::object( { {"functionType", LauDecayTime::FuncType::Exp}, {"parameters", json::array({"dt_tau_Lb"})} } )} } ) );
+
+	std::cout << j.dump(4) << std::endl;
+
+	auto models = LauDecayTimePhysicsModel::readFromJson( j );
+	for ( auto& [ name, model ] : models ) {
+		std::cout << "Loaded LauDecayTimePhysicsModel for component: " << name << " with type " << static_cast<int>(model->getFunctionType()) << std::endl;
+	}
+}