diff --git a/examples/GenFitTimeDep_Bd2Dpipi.cc b/examples/GenFitTimeDep_Bd2Dpipi.cc
index d3da8e1..7aab72e 100644
--- a/examples/GenFitTimeDep_Bd2Dpipi.cc
+++ b/examples/GenFitTimeDep_Bd2Dpipi.cc
@@ -1,934 +1,1005 @@
 
 #include "TimeDep_Dpipi_ProgOpts.hh"
 
 #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 "TFile.h"
 #include "TH2.h"
 #include "TRandom.h"
 #include "TString.h"
 #include "TSystem.h"
 #include "TF1.h"
 #include "TCanvas.h"
 
 #include "nlohmann/json.hpp"
 
 #include <algorithm>
 #include <filesystem>
 #include <iostream>
 #include <map>
 #include <tuple>
 #include <vector>
 
 // 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
+// Helper functions for creating sets of DTA splines
 using LauCubicSplineDecayTimeEfficiencyPtr = std::unique_ptr<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>>;
 using LauSixthSplineDecayTimeEfficiencyPtr = std::unique_ptr<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>;
+
 auto createSplineDTASet( const TString& jsonFileName, const TString& dModeStr, const TString& runStr, const TString& suffix, const TString& signalName, const std::vector<TString>& bkgndNames )
 {
+	std::cout << "INFO - creating QFS signal DTA model" << std::endl;
 	const TString controlSplineName { Form("%s_Kpi_%s%s",signalName.Data(),runStr.Data(),suffix.Data()) };
 	const TString controlEffName { Form("dteff_QFS_%s",runStr.Data()) };
 	auto controlSpline = Lau1DCubicSpline::readFromJson( jsonFileName, controlSplineName );
 	LauCubicSplineDecayTimeEfficiencyPtr controlEff { std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>>( controlEffName, std::move( controlSpline ) ) };
 
 	LauSixthSplineDecayTimeEfficiencyPtr signalEff;
 	if ( dModeStr != "Kpi" ) {
+		std::cout << "INFO - creating " << dModeStr << " signal DTA model" << std::endl;
 		const TString signalSplineName { Form("%s_%s_%s%s",signalName.Data(),dModeStr.Data(),runStr.Data(),suffix.Data())  };
 		const TString signalEffName { Form("dteff_CP%s_%s",dModeStr.Data(),runStr.Data()) };
 		auto signalSpline = Lau1DCubicSpline::readFromJson( jsonFileName, signalSplineName );
 		signalEff = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( signalEffName, *controlEff, std::move( signalSpline ) );
 	}
 
 	std::map<TString,LauSixthSplineDecayTimeEfficiencyPtr> bkgndEffs;
 	for ( auto& name : bkgndNames ) {
+		std::cout << "INFO - creating " << dModeStr << " " << name << " background DTA model" << std::endl;
 		const TString bkgndSplineName { Form("%s_%s_%s%s",name.Data(),dModeStr.Data(),runStr.Data(),suffix.Data()) };
 		const TString bkgndEffName { Form("%s_DTA_ratio_%s_%s",name.Data(),runStr.Data(),dModeStr.Data()) };
 		auto bkgndSpline = Lau1DCubicSpline::readFromJson( jsonFileName, bkgndSplineName );
 		auto bkgndEff = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( bkgndEffName, *controlEff, std::move( bkgndSpline ) );
 		bkgndEffs.emplace( std::make_pair( name, std::move(bkgndEff) ) );
 	}
 
 	return std::make_tuple( std::move(controlEff), std::move(signalEff), std::move(bkgndEffs) );
 }
 
+auto createSplineDTASetPCBkgnds( const TString& jsonFileName, const TString& dModeStr, const TString& runStr, const TString& suffix, const TString& signalName, const std::vector<TString>& bkgndNames )
+{
+	std::cout << "INFO - creating QFS Dpi DTA model" << std::endl;
+
+	const TString splineNameDpi { Form("Dpi_%s",runStr.Data()) };
+	auto splineDpi = Lau1DCubicSpline::readFromJson( jsonFileName, splineNameDpi );
+
+	// used in case we need a ratio that is just unity across the range
+	const std::vector<Double_t>& xVals { splineDpi->getXValues() };
+	std::vector<Double_t> yVals;
+	yVals.assign( xVals.size(), 1.0 );
+
+	const TString effNameDpi { Form("dteff_QFS_Dpi_%s",runStr.Data()) };
+	LauCubicSplineDecayTimeEfficiencyPtr effDpi { std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>>( effNameDpi, std::move( splineDpi ) ) };
+
+	std::map<TString,LauSixthSplineDecayTimeEfficiencyPtr> bkgndEffs;
+	for ( auto& name : bkgndNames ) {
+		std::cout << "INFO - creating " << dModeStr << " " << name << " background DTA model" << std::endl;
+
+		std::unique_ptr<Lau1DCubicSpline> signalRatioSpline;
+		if ( dModeStr == "Kpi" ) {
+			signalRatioSpline = std::make_unique<Lau1DCubicSpline>( xVals, yVals );
+		} else {
+			const TString signalRatioSplineName { Form("%s_%s_%s%s",signalName.Data(),dModeStr.Data(),runStr.Data(),suffix.Data())  };
+			signalRatioSpline = Lau1DCubicSpline::readFromJson( jsonFileName, signalRatioSplineName );
+		}
+
+		const TString bkgndEffName { Form("%s_DTA_ratio_%s_%s",name.Data(),runStr.Data(),dModeStr.Data()) };
+		auto bkgndEff = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( bkgndEffName, *effDpi, std::move( signalRatioSpline ) );
+		bkgndEffs.emplace( std::make_pair( name, std::move(bkgndEff) ) );
+	}
+
+	return std::make_tuple( std::move(effDpi), 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}
 		{}
 
 		BgTagDecayFlvAsym( const nlohmann::json& j ) : BgTagDecayFlvAsym( j[0], j[1], j[2], j[3] )
 		{}
 
 		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;
 	}
 
 	// check that all required data files are present
 	std::string bkgndsJsonFile { "Bd2D0pipi_Bkgnds" };
 	std::string signalYieldsJsonFile { "Bd2D0pipi_SignalYields" };
 	if ( settings.setNegToZero ) {
 		bkgndsJsonFile       += "_ZeroNegYields";
 		signalYieldsJsonFile += "_ZeroNegYields";
 	} else if ( ! settings.constrainedYields ) {
 		bkgndsJsonFile       += "_WithNegYields";
 		signalYieldsJsonFile += "_WithNegYields";
 	}
 	if ( ! settings.scaleYields ) {
 		bkgndsJsonFile += "_AddNeglectedToComb";
 	}
 	bkgndsJsonFile       += ".json";
 	signalYieldsJsonFile += ".json";
 	const std::string dpModelJsonFile { "Bd2D0pipi_DP_Model_Coeffs_Vetoes.json" };
 	const std::string dtModelJsonFile { "Bd2D0pipi_DecayTime.json" };
 	const std::string dtaModelJsonFile { "Bd2D0pipi_DecayTime_Acceptance.json" };
 
 	const std::array jsonFiles {
 		bkgndsJsonFile,
 		dpModelJsonFile,
 		dtModelJsonFile,
 		dtaModelJsonFile,
 		signalYieldsJsonFile
 	};
 
 	for ( auto& jsonFile : jsonFiles ) {
 		if ( std::filesystem::exists( jsonFile ) ) {
 			continue;
 		}
 
 		std::cerr << "WARNING : couldn't find the json file \"" << jsonFile << "\" attempting to fetch from eos ..." << std::endl;
 		if ( system( Form("xrdcp %s/%s/%s .", settings.eosRoot.c_str(), settings.eosConfigDir.c_str(), jsonFile.c_str()) ) != 0 ) {
 			std::cerr << "FATAL : couldn't get the json file from remote!" << std::endl;
 			return EXIT_FAILURE;
 		}
 	}
 
 
 	const UInt_t runNo = settings.run;
 	const std::string runStr = Form("Run%u", runNo);
 
 	const TString dModeStr { TString{settings.directory}.ReplaceAll("B2Dhh_","") };
 
 	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( dpModelJsonFile, "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
 
 	// NB deliberately not using uniform initialisation here because of this issue:
 	// https://json.nlohmann.me/home/faq/#brace-initialization-yields-arrays
 	const nlohmann::json bkgndsJson = LauJsonTools::readJsonFile( bkgndsJsonFile, "", LauJsonTools::JsonType::Object );
 	if ( bkgndsJson.is_null() ) {
 		std::cerr << "ERROR : unable to retrieve root JSON object from file \"" << bkgndsJsonFile << "\"" << std::endl;
 		return EXIT_FAILURE;
 	}
 
 	std::map<TString,LauFlavTag::BkgndType> BkgndTypes;
 	std::map<TString, LauDecayTime::FuncType> bkgndDecayTimeTypes;
 
 	for ( auto& bkgndName : settings.bkgndList ) {
 		if ( ! bkgndsJson.contains( bkgndName ) ) {
 			std::cerr << "ERROR : unable to retrieve \"" << bkgndName << "\" JSON object from file \"" << bkgndsJsonFile << "\"" << std::endl;
 			return EXIT_FAILURE;
 		}
 
 		auto bkgndType { LauJsonTools::getValue<LauFlavTag::BkgndType>( bkgndsJson.at(bkgndName), "type" ) };
 		auto bkgndDtFunc { LauJsonTools::getValue<LauDecayTime::FuncType>( bkgndsJson.at(bkgndName), "dt_func" ) };
 
 		BkgndTypes.insert( std::make_pair( bkgndName, bkgndType ) );
 		bkgndDecayTimeTypes.insert( std::make_pair( bkgndName, bkgndDtFunc ) );
 	}
 
 	// setup flavour tagging
 	const Bool_t useAveDeltaCalibVals { settings.useAveDeltaCalibVals };
 	const Bool_t useEtaPrime { kFALSE };
 	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();
+		TString decayName = name;
+		if ( name == "comb" ) {
+			decayName = "sideband";
+		} else if ( name.Contains("Dst") ) {
+			decayName = "signalMC";
+		}
+
 		if (settings.directory == "B2Dhh_Kpi"){
-			bkgndEtaHist1      = dynamic_cast<TH1*>(etaFile->Get( Form("%s_Kpi_Run%u_etaHist_dcyB" ,decayName,runNo) ));
+			bkgndEtaHist1      = dynamic_cast<TH1*>(etaFile->Get( Form("%s_Kpi_Run%u_etaHist_dcyB" ,decayName.Data(),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) ));
+			bkgndEtaHistMinus1 = dynamic_cast<TH1*>(etaFile->Get( Form("%s_Kpi_Run%u_etaHist_dcyBbar",decayName.Data(),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) ));
+			bkgndEtaHist1      = dynamic_cast<TH1*>(etaFile->Get( Form("%s_KK_Run%u_etaHist" ,decayName.Data(),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) ));
+			bkgndEtaHist1      = dynamic_cast<TH1*>(etaFile->Get( Form("%s_pipi_Run%u_etaHist" ,decayName.Data(),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 ( name == "comb" ) {
 			if ( ! bkgndsJson.contains("comb") || ! bkgndsJson.at("comb").contains("tag_asym") ) {
 				std::cerr << "ERROR : unable to retrieve \"tag_asym\" JSON object from \"comb\" object of file \"" << bkgndsJsonFile << "\"" << std::endl;
 				return EXIT_FAILURE;
 			}
 			const nlohmann::json& combTagAsymJson { bkgndsJson.at("comb").at("tag_asym") };
 
 			if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) {
 				const BgTagDecayFlvAsym bgt { combTagAsymJson.at(runStr).at(settings.directory) };
 				const std::pair<Double_t,Double_t> effsB0    { useAveDeltaCalibVals ? bgt.aveDeltaB0()    : bgt.decayB0()    };
 				const std::pair<Double_t,Double_t> effsB0bar { useAveDeltaCalibVals ? bgt.aveDeltaB0bar() : bgt.decayB0bar() };
 				flavTag->setBkgndParams( name, Form("IFT_Run%u",runNo),
 							 hists.first,  effsB0,
 							 hists.second, effsB0bar );
 			} else if (settings.dType == LauTimeDepFitModel::CPEigenvalue::CPEven){
 				const BgTagDecayFlvAsym bgt { combTagAsymJson.at(runStr).at(settings.directory) };
 				const std::pair<Double_t,Double_t> effs { useAveDeltaCalibVals ? bgt.aveDeltaB0() : bgt.decayB0() };
 				flavTag->setBkgndParams( name, Form("IFT_Run%u",runNo),
 							 hists.first, effs);
 			}
 		} else {
 			if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) {
 				flavTag->setBkgndParams( name, Form("IFT_Run%u",runNo),
 							 hists.first,  tagEff,
 							 hists.second, tagEff );
 			} else if (settings.dType == LauTimeDepFitModel::CPEigenvalue::CPEven){
 				flavTag->setBkgndParams( name, Form("IFT_Run%u",runNo),
 							 hists.first, tagEff );
 			}
 		}
 	}
 
 	// signal dynamics
 	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( dpModelJsonFile, "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( dpModelJsonFile, "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( dpModelJsonFile, "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* BkgndHistFile1 = 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()));
+	TFile* BkgndHistFile3 = TFile::Open(Form("%s/%s/BG_SDPs_Dst.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 if ( name.Contains("Dst") ) {
+			TString histNameB = Form("%s_Kpi_Run%u_B_SDP",name.Data(),runNo);
+			TString histNameBbar = Form("%s_Kpi_Run%u_Bbar_SDP",name.Data(),runNo);
+			BkgndHistB0 = dynamic_cast<TH2*>(BkgndHistFile3->Get(histNameB));
+			BkgndHistB0bar = dynamic_cast<TH2*>(BkgndHistFile3->Get(histNameBbar));
+			if ( !BkgndHistB0 || !BkgndHistB0bar ) {
+				std::cerr << "ERROR - unable to retrieve histograms " << histNameB << " and " << histNameBbar << " from file " << BkgndHistFile3->GetName() << std::endl;
+				return 1;
+			}
 		} 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"));
+			BkgndHistB0 = dynamic_cast<TH2*>(BkgndHistFile1->Get(histName+"1"));
+			BkgndHistB0bar = dynamic_cast<TH2*>(BkgndHistFile1->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( dtModelJsonFile, "PhysicsModels" ) };
 
 	auto& dtPhysModel { dtPhysicsModels.at("signal") };
 	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_sig_deltaM_Bd") );
 	if ( deltaMPar ) {
 		deltaMPar->fixed( settings.fixDeltaM );
 	}
 
 	// resolution models
 	auto dtResolutionModels { LauDecayTimeResolution::readFromJson( dtModelJsonFile, Form("ResolutionModels_Run%d", runNo) ) };
 
 	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;}
 
 	const TString signalName { "signalMC" };
 	std::vector<TString> bkgndNames;
+	std::vector<TString> bkgndNamesDst;
 	bkgndNames.reserve( bkgndDecayTimeTypes.size() );
+	bkgndNamesDst.reserve( bkgndDecayTimeTypes.size() );
 	for ( auto& [name, type] : bkgndDecayTimeTypes ) {
-		if ( type != LauDecayTime::FuncType::Hist ) {
+		if ( type == LauDecayTime::FuncType::Hist ) {
+			continue;
+		}
+
+		if ( name.Contains("Dst") ) {
+			bkgndNamesDst.push_back( name );
+		} else {
 			bkgndNames.push_back( name );
 		}
 	}
 	const TString suffix { "_DTAhist" };
 
 	auto [ qfsDTAModel, cpDTAModel, bkgndDTAModels ] = createSplineDTASet( dtaModelJsonFile, dModeStr, runStr, suffix, signalName, bkgndNames );
+	auto [ dtaDpiModel, pcBkgndDTAModels ] = createSplineDTASetPCBkgnds( dtaModelJsonFile, dModeStr, runStr, suffix, signalName, bkgndNamesDst );
 
 	switch(settings.timeEffModel)
 	{
 		case LauDecayTime::EfficiencyMethod::Spline:
 			{
 				fitModel->setASqMaxValue(0.01);
 				if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) {
 					fitModel->setASqMaxValue(0.009);
 				}
 
 				dynamic_cast<LauParameter*>(qfsDTAModel->getParameters().front())->maxValue(1e-1);
 				//qfsDTAModel->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 ) {
 					qfsDTAModel->fixKnots();
 				} else {
 					qfsDTAModel->floatKnots();
 					qfsDTAModel->fixKnot( 0, true );
 					qfsDTAModel->fixKnot( 3, true );
 				}
 
 				if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) {
 					dtPdf->setSplineEfficiency( std::move(qfsDTAModel) );
 				} else {
 					dtPdf->setSplineEfficiency( std::move(cpDTAModel) );
 				}
 
 				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 );
 
 	const TString dtaHistSuffix { Form("_%s_%s_DTAhist",dModeStr.Data(),runStr.c_str()) };
 
 	for (auto& [bg, _] : BkgndTypes) {
 
 		if ( bkgndDecayTimeTypes[bg] == LauDecayTime::FuncType::Hist ) {
 
 			TString histname = bg + dtaHistSuffix;
 			if(histname.Contains("comb"))
 			{
 				histname = "Sideband" + dtaHistSuffix;
 				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 {
 
+			try { dtPhysicsModels.at(bg); } catch ( std::out_of_range& ) { std::cerr << "ERROR - cannot locate physics model for " << bg << std::endl; return 1; }
+			try { dtResolutionModels.at(bg); } catch ( std::out_of_range& ) { std::cerr << "ERROR - cannot locate resolution model for " << bg << std::endl; return 1; }
+			try { ( bg.Contains("Dst") ) ? pcBkgndDTAModels.at(bg) : bkgndDTAModels.at(bg); } catch ( std::out_of_range& ) { std::cerr << "ERROR - cannot locate DTA model for " << bg << std::endl; return 1; }
+
 			auto& bgPhysModel { dtPhysicsModels.at(bg) };
 			auto& bgResoModel { dtResolutionModels.at(bg) };
-			auto& bgDTAModel  { bkgndDTAModels.at(bg) };
+			auto& bgDTAModel  { ( bg.Contains("Dst") ) ? pcBkgndDTAModels.at(bg) : bkgndDTAModels.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 + dtaHistSuffix;
 			auto 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.);}
 			}
 
-			//bgDTAModel->fitToTH1(bgDTA, LauOutputLevel::Quiet, true);
+			bgDTAModel->fitToTH1(bgDTA, LauOutputLevel::Quiet, true);
+			*/
+
 			backgroundPdf->setSplineEfficiency( std::move(bgDTAModel) );
 
 			fitModel->setBkgndDtPdf(bg,backgroundPdf);
 		}
 	}
 
 	// Signal and bkgnd yields
 	// NB deliberately not using uniform initialisation here because of this issue:
 	// https://json.nlohmann.me/home/faq/#brace-initialization-yields-arrays
 	const nlohmann::json signalYieldsJson = LauJsonTools::readJsonFile( signalYieldsJsonFile, "", LauJsonTools::JsonType::Object );
 	if ( signalYieldsJson.is_null() ) {
 		std::cerr << "ERROR : unable to retrieve root JSON object from file \"" << signalYieldsJsonFile << "\"" << std::endl;
 		return EXIT_FAILURE;
 	}
 
 	const Double_t nSigEvents{ LauJsonTools::getValue<Double_t>( signalYieldsJson.at(runStr), settings.directory ) };
 
 	std::map<TString, Double_t> bkgndYields;
 	for ( auto& bkgndName : settings.bkgndList ) {
 		const auto& yieldsJson { bkgndsJson.at(bkgndName).at("yields") };
 		const Double_t nEvents{ LauJsonTools::getValue<Double_t>( yieldsJson.at(runStr), settings.directory ) };
 		bkgndYields.insert( std::make_pair( bkgndName, nEvents ) );
 	}
 
 	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}
+			{"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}
+				{"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}
+				{"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}
+				{"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, kTRUE);
 	} else {
 		fitModel->setNExpts(settings.nExptFit, settings.firstExptFit, !settings.dataFit);
 	}
 
 	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/ProgOpts/TimeDep_Dpipi_ProgOpts.cc b/examples/ProgOpts/TimeDep_Dpipi_ProgOpts.cc
index 9ac5d3f..4979246 100644
--- a/examples/ProgOpts/TimeDep_Dpipi_ProgOpts.cc
+++ b/examples/ProgOpts/TimeDep_Dpipi_ProgOpts.cc
@@ -1,174 +1,174 @@
 
 #include <iostream>
 #include <string>
 
 #include "boost/program_options.hpp"
 #include "boost/tokenizer.hpp"
 
 #include "TimeDep_Dpipi_ProgOpts.hh"
 
 namespace po = boost::program_options;
 
 void validate( boost::any& v, const std::vector<std::string>& values, Command* /*target_type*/, int)
 {
 	// Make sure no previous assignment to 'v' has been made
 	po::validators::check_first_occurrence(v);
 
 	// Extract the first string from 'values'. If there is more than
 	// one string, it's an error, and exception will be thrown.
 	const std::string& s { po::validators::get_single_string(values) };
 
 	if ( s == "gen" ) {
 		v = boost::any( Command::Generate );
 	} else if ( s == "fit" ) {
 		v = boost::any( Command::Fit );
 	} else if ( s == "simfit" ) {
 		v = boost::any( Command::SimFit );
 	} else {
 		throw po::validation_error(po::validation_error::invalid_option_value, "command", s, 3);
 	}
 }
 
 void validate( boost::any& v, const std::vector<std::string>& values, LauTimeDepFitModel::CPEigenvalue* /*target_type*/, int)
 {
 	// Make sure no previous assignment to 'v' has been made
 	po::validators::check_first_occurrence(v);
 
 	// Extract the first string from 'values'. If there is more than
 	// one string, it's an error, and exception will be thrown.
 	const std::string& s { po::validators::get_single_string(values) };
 
 	if ( s == "QFS" ) {
 		v = boost::any( LauTimeDepFitModel::CPEigenvalue::QFS );
 	} else if ( s == "CPEven" ) {
 		v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPEven );
 	} else if ( s == "CPOdd" ) {
 		v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPOdd );
 	} else {
 		throw po::validation_error(po::validation_error::invalid_option_value);
 	}
 }
 
 namespace LauDecayTime {
 void validate( boost::any& v, const std::vector<std::string>& values, EfficiencyMethod* /*target_type*/, int)
 {
 	// Make sure no previous assignment to 'v' has been made
 	po::validators::check_first_occurrence(v);
 
 	// Extract the first string from 'values'. If there is more than
 	// one string, it's an error, and exception will be thrown.
 	const std::string& s { po::validators::get_single_string(values) };
 
 	if ( s == "uniform" || s == "flat" ) {
 		v = boost::any( EfficiencyMethod::Uniform );
 	} else if ( s == "binned" || s == "hist" ) {
 		v = boost::any( EfficiencyMethod::Binned );
 	} else if ( s == "spline" ) {
 		v = boost::any( EfficiencyMethod::Spline );
 	} else {
 		throw po::validation_error(po::validation_error::invalid_option_value);
 	}
 }
 };
 
 TestDpipi_ProgramSettings::TestDpipi_ProgramSettings(const int argc, const  char ** argv)
 {
 	po::options_description main_desc{"Main options"};
 	main_desc.add_options()
 		("command", po::value<Command>(&command)->required(), "main command: gen, fit, or simfit")
 		;
 	po::positional_options_description p;
 	p.add("command", 1);
 
 	po::options_description common_desc{"Common options"};
 	common_desc.add_options()
 		("help", "produce help message")
 		("dtype", po::value<LauTimeDepFitModel::CPEigenvalue>(&dType)->default_value(LauTimeDepFitModel::CPEigenvalue::QFS,"QFS"), "type of D decay: QFS, CPOdd, or CPEven")
 		("dta-model", po::value<LauDecayTime::EfficiencyMethod>(&timeEffModel)->default_value(LauDecayTime::EfficiencyMethod::Uniform,"uniform"), "decay-time acceptance model: uniform/flat, binned/hist, spline")
 		("dtr", po::value<Bool_t>(&timeResolution)->default_value(kTRUE), "enable/disable decay-time resolution")
 		("dtr-perevent", po::value<Bool_t>(&perEventTimeErr)->default_value(kFALSE), "enable/disable use of per-event decay-time error (requires decay-time resolution to be enabled to take effect)")
 		("scaleYields", po::value<Bool_t>(&scaleYields)->default_value(kTRUE), "toggle whether (kTRUE) the yields are scaled to account for the neglected components (the default strategy) or (kFALSE) the neglected components' yields are added to the combinatorial yield")
 		("setNegToZero", po::value<Bool_t>(&setNegToZero)->default_value(kFALSE), "toggle whether (kTRUE) the negative yields from the mass fits for Run 1 KK are set to 0")
 		("constrainedYields", po::value<Bool_t>(&constrainedYields)->default_value(kTRUE), "toggle whether (kTRUE) the yields in the mass fits for Run 1 KK are constrained based on the other samples")
 		("seed", po::value<UInt_t>(&RNGseed)->default_value(0), "set the seed for the RNG; if not set, the time is used to generate a seed")
 		("run", po::value<UInt_t>(&run)->default_value(2), "set the Run number (1 or 2)")
 		("dir", po::value<std::string>(&directory)->default_value("<none>"), "set the directory used to find the nTuples within the root file. Defaults to no directory")
-		("bkgndList", po::value<std::string>(&bkgndListStr)->default_value("comb,Bd2DKpi,Bs2DKpi,Lb2Dppi"), "the comma-separated list of background components to include")
+		("bkgndList", po::value<std::string>(&bkgndListStr)->default_value("comb,Bd2DKpi,Bs2DKpi,Lb2Dppi,Bd2Dstpi,Bu2Dstpi,Bu2Dstpi_Dpi0"), "the comma-separated list of background components to include")
 		("fitBackInput", po::value<std::string>(&fitBackInput)->default_value("fit0_ToyMC_QFS_expts0-0_expt0.root"), "set the file used as input to the fit: only used if fitBack is also selected")
 		("eosRoot", po::value<std::string>(&eosRoot)->default_value("root://eoslhcb.cern.ch//eos/lhcb/wg/b2oc/B2D0pipi"), "set the base EOS area from which all the input data/histograms/json files/etc. should be read")
 		("eosConfigDir", po::value<std::string>(&eosConfigDir)->default_value("timedep_fit_histos"), "set the directory name from which the configuration (json files / histogram files) should be read")
 		("eosDataDir", po::value<std::string>(&eosDataDir)->default_value("Tuples22/Data"), "set the directory name from which the data should be read")
 		;
 
 	po::options_description gen_desc{"Generation options"};
 	gen_desc.add_options()
 		("firstExptGen", po::value<UInt_t>(&firstExptGen)->default_value(0), "first experiment to generate")
 		("nExptGen", po::value<UInt_t>(&nExptGen)->default_value(1), "number of experiments to generate")
 		;
 
 	po::options_description fit_desc{"Fitting options"};
 	fit_desc.add_options()
 		("firstExpt", po::value<UInt_t>(&firstExptFit)->default_value(0), "first experiment to fit")
 		("nExpt", po::value<UInt_t>(&nExptFit)->default_value(1), "number of experiments to fit")
 		("iFit", po::value<UInt_t>(&iFit)->default_value(0), "fit ID - used to distinguish multiple fits to a given dataset (used to disentangle multiple minima)")
 		("fixTau", po::value<Bool_t>(&fixLifetime)->default_value(kTRUE), "enable/disable floating of B lifetime parameter")
 		("fixDm", po::value<Bool_t>(&fixDeltaM)->default_value(kTRUE), "enable/disable floating of B mixing frequency parameter")
 		("fixPhiMix", po::value<Bool_t>(&fixPhiMix)->default_value(kFALSE), "enable/disable floating of B mixing phase parameter(s)")
 		("fixSplineKnots", po::value<Bool_t>(&fixSplineKnots)->default_value(kFALSE), "enable/disable floating of the decay-time-acceptance spline")
 		("useSinCos", po::value<Bool_t>(&useSinCos)->default_value(kTRUE), "enable/disable using sine and cosine of B mixing phase as the floating parameters rather than the phase itself")
 		("useAveDeltaCalibVals", po::value<Bool_t>(&useAveDeltaCalibVals)->default_value(kTRUE), "enable/disable fitting calib values as their average +/- delta values rather than having separate values for B and Bbar")
 		("addTaggers", po::value<Bool_t>(&addTaggers)->default_value(kTRUE), "add/omit taggers")
 		("floatCalibPars", po::value<Bool_t>(&floatCalibPars)->default_value(kTRUE), "enable/disable floating of the FT calibration parameters")
 		("floatYields", po::value<Bool_t>(&floatYields)->default_value(kFALSE), "enable/disable floating of the yields of each component - if enabled/disabled implies inclusion/exclusion of non-DP PDFs in/from the likelihood")
 		("fitBack", po::value<Bool_t>(&fitBack)->default_value(kFALSE), "enable/disable refitting some generated data from a model")
 		("blindFit", po::value<Bool_t>(&blindFit)->default_value(kFALSE), "enable/disable blinding of the physics parameters")
 		("dataFit", po::value<Bool_t>(&dataFit)->default_value(kFALSE), "fit to data (not MC or toys)")
 		("genToy", po::value<Bool_t>(&genToy)->default_value(kFALSE), "enable/disable generation of a toy to a fit result")
 		;
 
 	po::options_description simfit_desc{"Simultaneous fitting options"};
 	simfit_desc.add_options()
 		("port", po::value<UInt_t>(&port)->default_value(0), "port number on which sim fit coordinator is listening")
 		;
 
 	po::options_description all_desc;
 	all_desc.add(main_desc).add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc);
 
 	po::variables_map vm;
 	try {
 		po::store(po::command_line_parser(argc, argv).
 				options(all_desc).
 				positional(p).
 				run(),
 				vm);
 		po::notify(vm);
 	}
 	catch ( boost::wrapexcept<po::required_option>& e ) {
 		std::cout << argv[0] << " requires a command, one of 'gen', 'fit', or 'simfit'\n"
 		          << "Append --help to those commands to see help on the related options" << std::endl;
 		parsedOK = kFALSE;
 		return;
 	}
 	catch ( po::validation_error& e ) {
 		std::cerr << e.what() << std::endl;
 		parsedOK = kFALSE;
 		return;
 	}
 
 	if ( vm.count("help") ) {
 		po::options_description help_desc;
 		help_desc.add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc);
 
 		std::cout << help_desc << std::endl;
 		helpRequested = kTRUE;
 		return;
 	}
 
 	if ( ! bkgndListStr.empty() ) {
 		const boost::char_separator<char> sep{","};
 		const boost::tokenizer<boost::char_separator<char>> tokens{bkgndListStr, sep};
 		for (const auto& t : tokens) {
 			bkgndList.insert( t );
 		}
 	}
 }