diff --git a/examples/GenFitTimeDep_Bd2Dpipi.cc b/examples/GenFitTimeDep_Bd2Dpipi.cc
index 649ad21..38a40ff 100644
--- a/examples/GenFitTimeDep_Bd2Dpipi.cc
+++ b/examples/GenFitTimeDep_Bd2Dpipi.cc
@@ -1,1082 +1,1084 @@
 
 #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
 	LauVetoes* vetoes = new LauVetoes();
+	vetoes->addMassVeto( 1, 2.0, 2.1 ); //D* veto
+	vetoes->addMassVeto( 2, 2.0, 2.1 ); //D* veto
 	LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes);
 	LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes);
 
 	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);
 		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);
 			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
 	};
 
 	auto dtPhysModel = std::make_unique<LauDecayTimePhysicsModel>( LauDecayTime::FuncType::ExpTrig, dtPhysPars );
 
 	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 );
 
 			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;
 }