No OneTemporary
Actions

Size

87 KB

Subscribers

None

View Options

	diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
	index 521aac2..66531ac 100644
	--- a/examples/CMakeLists.txt
	+++ b/examples/CMakeLists.txt
	@@ -1,49 +1,50 @@

	add_subdirectory(ProgOpts)

	list(APPEND EXAMPLE_SOURCES
	B3piKMatrixMassProj
	B3piKMatrixPlots
	CalcChiSq
	GenFit3K
	GenFit3KS
	GenFit3pi
	GenFitBelleCPKpipi
	GenFitDpipi
	GenFitDs2KKpi
	GenFitEFKLLM
	GenFitIncoherent_Bs2KSKpi
	GenFitKpipi
	GenFitNoDP
	GenFitNoDPMultiDim
	GenFitTimeDep
	GenFitTimeDep_Bs2KSKpi
	+ GenFitTimeDep_Bd2Dpipi
	KMatrixDto3pi
	KMatrixExample
	MergeDataFiles
	mixedSampleTest
	PlotKMatrixTAmp
	PlotResults
	point2PointTestSample
	QuasiFlatSqDalitz
	QuasiFlatSqDalitz-CustomMasses
	ResultsExtractor
	SimFitCoordinator
	SimFitTask
	SimFitTaskRooFit
	)

	if(NOT LAURA_BUILD_ROOFIT_TASK)
	list(REMOVE_ITEM EXAMPLE_SOURCES SimFitTaskRooFit)
	endif()

	foreach( _example ${EXAMPLE_SOURCES})
	add_executable(${_example} ${_example}.cc)
	target_link_libraries(${_example} PRIVATE Laura++ ProgramOptions)
	install(TARGETS ${_example} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
	endforeach()

	# Also install the python script version of GenFit3pi
	configure_file(GenFit3pi.py.in GenFit3pi.py @ONLY)
	install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/GenFit3pi.py DESTINATION ${CMAKE_INSTALL_BINDIR})

	diff --git a/examples/GenFitTimeDep_Bd2Dpipi.cc b/examples/GenFitTimeDep_Bd2Dpipi.cc
	new file mode 100644
	index 0000000..649ad21
	--- /dev/null
	+++ b/examples/GenFitTimeDep_Bd2Dpipi.cc
	@@ -0,0 +1,1082 @@
	+
	+#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();
	+ 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.0TMath::Abs(nBkg),2.0TMath::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.0nSigEvents, 2.0nSigEvents, fixYields);
	+ fitModel->setNSigEvents(nSigPar);
	+
	+ if ( settings.command == Command::Generate or not fixYields ) {
	+ // mB PDFs
	+ // PDFs for B+ and B- are set to be identical here but could be different in general
	+ const TString mBName { "mB" };
	+ const Double_t mBMin { 5.200 };
	+ const Double_t mBMax { 5.700 };
	+
	+ // Signal PDF is a double Crystal Ball function
	+ const std::map<TString,Double_t> sigParVals {
	+ {"mB_sig_mean", 5.2802},
	+ {"mB_sig_sigma", 0.01125},
	+ {"mB_sig_alphaL", 1.226},
	+ {"mB_sig_alphaR", -1.722},
	+ {"mB_sig_orderL", 1.817},
	+ {"mB_sig_orderR", 2.861},
	+ {"mB_sig_frac", 0.480}
	+ };
	+ fitModel->setSignalPdfs( makeDCBPDF( mBName, mBMin, mBMax, "sig", sigParVals ) );
	+
	+ // Combinatorial PDF is an exponential
	+ if ( BkgndTypes.find("comb") != BkgndTypes.end() ) {
	+ LauParameter* mB_comb_slope = new LauParameter("mB_comb_slope", -0.00419);
	+ LauAbsPdf* mB_comb_pdf = new LauExponentialPdf(mBName, {mB_comb_slope}, mBMin, mBMax);
	+ fitModel->setBkgndPdf("comb", mB_comb_pdf);
	+ }
	+
	+ // Bd2DKpi PDF is double Crystal Ball function
	+ if ( BkgndTypes.find("Bd2DKpi") != BkgndTypes.end() ) {
	+ const std::map<TString,Double_t> Bd2DKpiParVals {
	+ {"mB_Bd2DKpi_mean", 5.2401},
	+ {"mB_Bd2DKpi_sigma", 0.01580},
	+ {"mB_Bd2DKpi_alphaL", 0.560},
	+ {"mB_Bd2DKpi_alphaR", -2.212},
	+ {"mB_Bd2DKpi_orderL", 1.076},
	+ {"mB_Bd2DKpi_orderR", 2.453},
	+ {"mB_Bd2DKpi_frac", 0.544}
	+ };
	+ fitModel->setBkgndPdf("Bd2DKpi", makeDCBPDF( mBName, mBMin, mBMax, "Bd2DKpi", Bd2DKpiParVals ) );
	+ }
	+
	+ // Bs2DKpi PDF is double Crystal Ball function
	+ if ( BkgndTypes.find("Bs2DKpi") != BkgndTypes.end() ) {
	+ const std::map<TString,Double_t> Bs2DKpiParVals {
	+ {"mB_Bs2DKpi_mean", 5.3244},
	+ {"mB_Bs2DKpi_sigma", 0.01827},
	+ {"mB_Bs2DKpi_alphaL", 0.865},
	+ {"mB_Bs2DKpi_alphaR", -0.001},
	+ {"mB_Bs2DKpi_orderL", 10.000},
	+ {"mB_Bs2DKpi_orderR", 10.000},
	+ {"mB_Bs2DKpi_frac", 0.954}
	+ };
	+ fitModel->setBkgndPdf("Bs2DKpi", makeDCBPDF( mBName, mBMin, mBMax, "Bs2DKpi", Bs2DKpiParVals ) );
	+ }
	+
	+ // Lb2Dppi PDF is double Crystal Ball function
	+ if ( BkgndTypes.find("Lb2Dppi") != BkgndTypes.end() ) {
	+ const std::map<TString,Double_t> Lb2DppiParVals {
	+ {"mB_Lb2Dppi_mean", 5.4884},
	+ {"mB_Lb2Dppi_sigma", 0.02561},
	+ {"mB_Lb2Dppi_alphaL", 0.091},
	+ {"mB_Lb2Dppi_alphaR", -18.053},
	+ {"mB_Lb2Dppi_orderL", 2.906},
	+ {"mB_Lb2Dppi_orderR", 6.087},
	+ {"mB_Lb2Dppi_frac", 0.968}
	+ };
	+ fitModel->setBkgndPdf("Lb2Dppi", makeDCBPDF( mBName, mBMin, mBMax, "Lb2Dppi", Lb2DppiParVals ) );
	+ }
	+ }
	+
	+ // set the number of experiments
	+ if (settings.command == Command::Generate) {
	+ fitModel->setNExpts(settings.nExptGen, settings.firstExptGen);
	+ } else {
	+ fitModel->setNExpts(settings.nExptFit, settings.firstExptFit);
	+ }
	+
	+ fitModel->useAsymmFitErrors(kFALSE);
	+ fitModel->useRandomInitFitPars(kFALSE);
	+ fitModel->writeLatexTable(kFALSE);
	+ const Bool_t haveBkgnds = ( fitModel->nBkgndClasses() > 0 );
	+ fitModel->doPoissonSmearing(haveBkgnds);
	+ fitModel->doEMLFit(haveBkgnds);
	+
	+ TString dTypeStr;
	+ switch (settings.dType) {
	+ case LauTimeDepFitModel::CPEigenvalue::CPEven :
	+ dTypeStr = "CPEven";
	+ break;
	+ case LauTimeDepFitModel::CPEigenvalue::CPOdd :
	+ dTypeStr = "CPOdd";
	+ break;
	+ case LauTimeDepFitModel::CPEigenvalue::QFS :
	+ dTypeStr = "QFS";
	+ break;
	+ }
	+
	+ TString dataFile("");
	+ TString treeName("fitTree");
	+ TString rootFileName("");
	+ TString tableFileName("");
	+ TString fitToyFileName("");
	+ TString splotFileName("");
	+
	+ dataFile = "TEST-Dpipi";
	+
	+ dataFile += "_"+dTypeStr+"_"+settings.directory;
	+
	+ switch(settings.timeEffModel)
	+ {
	+ case LauDecayTime::EfficiencyMethod::Spline:
	+ dataFile += "_Spline";
	+ break;
	+ case LauDecayTime::EfficiencyMethod::Binned:
	+ dataFile += "_Binned";
	+ break;
	+ case LauDecayTime::EfficiencyMethod::Uniform:
	+ dataFile += "_Uniform";
	+ break;
	+ }
	+
	+ if (settings.timeResolution) {
	+ if (settings.perEventTimeErr) {
	+ dataFile += "_DTRperevt";
	+ } else {
	+ dataFile += "_DTRavg";
	+ }
	+ } else {
	+ dataFile += "_DTRoff";
	+ }
	+
	+ dataFile += "_expts";
	+ dataFile += settings.firstExptGen;
	+ dataFile += "-";
	+ dataFile += settings.firstExptGen+settings.nExptGen-1;
	+ dataFile += Form("-Run%u",runNo);
	+
	+ dataFile += ".root";
	+
	+ if (settings.dataFit)
	+ {
	+ if(runNo == 2)
	+ {
	+ dataFile = Form("%s/%s/Run2/B2Dhh-CollisionCombined-MagCombined-Stripping24r2-withDNN-withPIDcorr-withMVA-withPIDMVA-WithDpiMatching-withSwappedMassHypotheses-withKKvetoes-withSelection-LauraPrepped.root",settings.eosRoot.c_str(),settings.eosDataDir.c_str());
	+ } else {
	+ dataFile = Form("%s/%s/Run1/B2Dhh-CollisionCombined-MagCombined-Stripping21r1p2-withDNN-withPIDcorr-withMVA-withPIDMVA-WithDpiMatching-withSwappedMassHypotheses-withKKvetoes-withSelection-LauraPrepped.root",settings.eosRoot.c_str(),settings.eosDataDir.c_str());
	+ }
	+ treeName = Form("%s/B2DhhReco",settings.directory.data());
	+ }
	+
	+ if (settings.command == Command::Generate) {
	+ rootFileName = "dummy.root";
	+ tableFileName = "genResults";
	+ } else {
	+ rootFileName = "fit"; rootFileName += settings.iFit;
	+ rootFileName += "_Results_"; rootFileName += dTypeStr;
	+ rootFileName += "_"; rootFileName += settings.directory;
	+ rootFileName += "_expts"; rootFileName += settings.firstExptFit; rootFileName += "-"; rootFileName += settings.firstExptFit+settings.nExptFit-1;
	+ rootFileName += Form("_Run%u",runNo);
	+ rootFileName += ".root";
	+
	+ tableFileName = "fit"; tableFileName += settings.iFit;
	+ tableFileName += "_Results_"; tableFileName += dTypeStr;
	+ tableFileName += "_"; tableFileName += settings.directory;
	+ tableFileName += "_expts"; tableFileName += settings.firstExptFit; tableFileName += "-"; tableFileName += settings.firstExptFit+settings.nExptFit-1;
	+ tableFileName += Form("_Run%u",runNo);
	+
	+ fitToyFileName = "fit"; fitToyFileName += settings.iFit;
	+ fitToyFileName += "_ToyMC_"; fitToyFileName += dTypeStr;
	+ fitToyFileName += "_"; fitToyFileName += settings.directory;
	+ fitToyFileName += "_expts"; fitToyFileName += settings.firstExptFit; fitToyFileName += "-"; fitToyFileName += settings.firstExptFit+settings.nExptFit-1;
	+ fitToyFileName += Form("_Run%u",runNo);
	+ fitToyFileName += ".root";
	+
	+ splotFileName = "fit"; splotFileName += settings.iFit;
	+ splotFileName += "_sPlot_"; splotFileName += dTypeStr;
	+ splotFileName += "_"; splotFileName += settings.directory;
	+ splotFileName += "_expts"; splotFileName += settings.firstExptFit; splotFileName += "-"; splotFileName += settings.firstExptFit+settings.nExptFit-1;
	+ splotFileName += Form("_Run%u",runNo);
	+ splotFileName += ".root";
	+ }
	+
	+ if( settings.genToy and not settings.blindFit )
	+ {
	+ fitModel->compareFitData(50, fitToyFileName);
	+ }
	+
	+ // Generate toy from the fitted parameters
	+ //fitModel->compareFitData(1, fitToyFileName);
	+
	+ // Write out per-event likelihoods and sWeights
	+ //fitModel->writeSPlotData(splotFileName, "splot", kFALSE);
	+
	+ // Execute the generation/fit
	+ switch (settings.command) {
	+ case Command::Generate :
	+ fitModel->run( "gen", dataFile, treeName, rootFileName, tableFileName );
	+ break;
	+ case Command::Fit :
	+ fitModel->run( "fit", dataFile, treeName, rootFileName, tableFileName );
	+ break;
	+ case Command::SimFit :
	+ fitModel->runTask( dataFile, treeName, rootFileName, tableFileName, "localhost", settings.port );
	+ break;
	+ }
	+
	+ return EXIT_SUCCESS;
	+}
	+
	+LauAbsPdf* makeDCBPDF( const TString& varName, const Double_t varMin, const Double_t varMax, const TString& componentName, const std::map<TString,Double_t> paramVals )
	+{
	+ const TString meanName { Form( "%s_%s_mean", varName.Data(), componentName.Data() ) };
	+ const TString sigmaName { Form( "%s_%s_sigma", varName.Data(), componentName.Data() ) };
	+ const TString alphaLName { Form( "%s_%s_alphaL", varName.Data(), componentName.Data() ) };
	+ const TString alphaRName { Form( "%s_%s_alphaR", varName.Data(), componentName.Data() ) };
	+ const TString orderLName { Form( "%s_%s_orderL", varName.Data(), componentName.Data() ) };
	+ const TString orderRName { Form( "%s_%s_orderR", varName.Data(), componentName.Data() ) };
	+ const TString fracName { Form( "%s_%s_frac", varName.Data(), componentName.Data() ) };
	+
	+ LauParameter* mean = new LauParameter( meanName, paramVals.at( meanName ) );
	+ LauParameter* sigma = new LauParameter( sigmaName, paramVals.at( sigmaName ) );
	+ LauParameter* alphaL = new LauParameter( alphaLName, paramVals.at( alphaLName ) );
	+ LauParameter* alphaR = new LauParameter( alphaRName, paramVals.at( alphaRName ) );
	+ LauParameter* orderL = new LauParameter( orderLName, paramVals.at( orderLName ) );
	+ LauParameter* orderR = new LauParameter( orderRName, paramVals.at( orderRName ) );
	+ LauParameter* frac = new LauParameter( fracName, paramVals.at( fracName ) );
	+
	+ std::vector<LauAbsRValue*> pars;
	+ pars.reserve(4);
	+
	+ pars.clear();
	+ pars.push_back(mean);
	+ pars.push_back(sigma);
	+ pars.push_back(alphaL);
	+ pars.push_back(orderL);
	+ LauAbsPdf* pdf_L = new LauCrystalBallPdf(varName, pars, varMin, varMax);
	+
	+ pars.clear();
	+ pars.push_back(mean);
	+ pars.push_back(sigma);
	+ pars.push_back(alphaR);
	+ pars.push_back(orderR);
	+ LauAbsPdf* pdf_R = new LauCrystalBallPdf(varName, pars, varMin, varMax);
	+
	+ LauAbsPdf* pdf = new LauSumPdf(pdf_L, pdf_R, frac);
	+
	+ return pdf;
	+}
	diff --git a/examples/ProgOpts/CMakeLists.txt b/examples/ProgOpts/CMakeLists.txt
	index 24de35c..30fa487 100644
	--- a/examples/ProgOpts/CMakeLists.txt
	+++ b/examples/ProgOpts/CMakeLists.txt
	@@ -1,7 +1,7 @@

	find_package(Boost REQUIRED COMPONENTS program_options)

	-add_library(ProgramOptions STATIC Test_Dpipi_ProgOpts.cc Test_Dpipi_ProgOpts.hh Command.hh)
	+add_library(ProgramOptions STATIC TimeDep_Dpipi_ProgOpts.cc TimeDep_Dpipi_ProgOpts.hh Command.hh)
	target_include_directories(ProgramOptions PUBLIC ${CMAKE_CURRENT_LIST_DIR})
	target_link_libraries(ProgramOptions PUBLIC Laura++ Boost::boost Boost::program_options)

	diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.cc b/examples/ProgOpts/TimeDep_Dpipi_ProgOpts.cc
	similarity index 99%
	rename from examples/ProgOpts/Test_Dpipi_ProgOpts.cc
	rename to examples/ProgOpts/TimeDep_Dpipi_ProgOpts.cc
	index a30c73c..9ac5d3f 100644
	--- a/examples/ProgOpts/Test_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 "Test_Dpipi_ProgOpts.hh"
	+#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")
	("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 );
	}
	}
	}
	diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.hh b/examples/ProgOpts/TimeDep_Dpipi_ProgOpts.hh
	similarity index 100%
	rename from examples/ProgOpts/Test_Dpipi_ProgOpts.hh
	rename to examples/ProgOpts/TimeDep_Dpipi_ProgOpts.hh
	diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc
	deleted file mode 100644
	index 5a98ab0..0000000
	--- a/examples/Test_Dpipi.cc
	+++ /dev/null
	@@ -1,448 +0,0 @@
	-
	-#include <iostream>
	-
	-#include <vector>
	-#include <map>
	-#include <algorithm>
	-
	-#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 "LauDaughters.hh"
	-#include "LauDecayTime.hh"
	-#include "LauDecayTimePdf.hh"
	-#include "LauDecayTimePhysicsModel.hh"
	-#include "LauDecayTimeResolution.hh"
	-#include "LauEffModel.hh"
	-#include "LauFlavTag.hh"
	-#include "LauIsobarDynamics.hh"
	-#include "LauMagPhaseCoeffSet.hh"
	-#include "LauRandom.hh"
	-#include "LauRealImagCoeffSet.hh"
	-#include "LauSplineDecayTimeEfficiency.hh"
	-#include "LauTimeDepFitModel.hh"
	-#include "LauVetoes.hh"
	-
	-#include "Test_Dpipi_ProgOpts.hh"
	-
	-int main(const int argc, const char ** argv)
	-{
	- const TestDpipi_ProgramSettings settings{argc,argv};
	- if ( settings.helpRequested ) {
	- return EXIT_SUCCESS;
	- }
	- if ( ! settings.parsedOK ) {
	- return EXIT_FAILURE;
	- }
	-
	- LauRandom::setSeed(settings.RNGseed);
	-
	- LauDaughters* daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0");
	- LauDaughters* daughtersB0 = new LauDaughters("B0", "pi+", "pi-", "D0_bar");
	-
	- // efficiency
	- LauVetoes* vetoes = new LauVetoes();
	- LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes);
	- LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes);
	-
	- // background types
	- /*
	- std::map<TString,LauFlavTag::BkgndType> bkgndInfo {
	- {"Bkgnd1",LauFlavTag::BkgndType::Combinatorial},
	- {"Bkgnd2",LauFlavTag::BkgndType::SignalLike}
	- };
	- */
	-
	- // 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");
	- }
	- flavTag->setTrueTagVarName("trueTag");
	-
	- TFile* etaFile = TFile::Open("ft-eta-hist.root");
	- TH1* etaHist = dynamic_cast<TH1*>(etaFile->Get("ft_eta_hist"));
	-
	- // 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
	- const Double_t meanEta { etaHist->GetMean() };
	- const Bool_t interpolateEtaHist { meanEta > etaHist->GetBinWidth(1) };
	- Lau1DHistPdf* etaHistPDF = new Lau1DHistPdf( "eta", etaHist, 0.0, 0.5, interpolateEtaHist, kFALSE );
	-
	- const Double_t tagEffVal { (interpolateEtaHist) ? 0.5 : 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("OSTagger", "tagVal_OS", "mistagVal_OS", etaHistPDF, tagEff, calib0, calib1);
	- if (settings.floatCalibPars) {
	- flavTag->floatAllCalibPars();
	- }
	- }
	-
	- // flavour tagging setup for backgrounds
	- /*
	- std::map<TString, Lau1DHistPdf*> BkgndEtas;
	- TFile* bkgndEtaFile = TFile::Open("ft-bkgnd-eta-hists.root");
	- for(auto& [name, _] : BkgndTypes)
	- {
	- TH1* bkgndEtaHist = dynamic_cast<TH1*>(bkgndEtaFile->Get( name+"_eta_hist" ));
	- Lau1DHistPdf* bkgndEtaHistPDF = new Lau1DHistPdf("eta",bkgndEtaHist,0.0,0.5,kTRUE,kFALSE);
	- BkgndEtas.emplace( std::make_pair(name, bkgndEtaHistPDF) );
	- }
	- for(auto& [name,hist] : BkgndEtas)
	- {flavTag->setBkgndParams(name, "IFT", hist, tagEff );}
	- */
	-
	- // signal dynamics
	- LauIsobarDynamics* sigModelB0bar = new LauIsobarDynamics(daughtersB0bar, effModelB0bar);
	- sigModelB0bar->setIntFileName("integ_B0bar.dat");
	- sigModelB0bar->addResonance("D*+_2", 2, LauAbsResonance::RelBW);
	- sigModelB0bar->addResonance("D*+_0", 2, LauAbsResonance::RelBW);
	- sigModelB0bar->addResonance("rho0(770)", 3, LauAbsResonance::RelBW);
	- sigModelB0bar->addResonance("f_0(980)", 3, LauAbsResonance::RelBW);
	- sigModelB0bar->addResonance("f_2(1270)", 3, LauAbsResonance::RelBW);
	-
	- LauIsobarDynamics* sigModelB0 = new LauIsobarDynamics(daughtersB0, effModelB0);
	- sigModelB0->setIntFileName("integ_B0.dat");
	- sigModelB0->addResonance("D*-_2", 1, LauAbsResonance::RelBW);
	- sigModelB0->addResonance("D*-_0", 1, LauAbsResonance::RelBW);
	- sigModelB0->addResonance("rho0(770)", 3, LauAbsResonance::RelBW);
	- sigModelB0->addResonance("f_0(980)", 3, LauAbsResonance::RelBW);
	- sigModelB0->addResonance("f_2(1270)", 3, LauAbsResonance::RelBW);
	-
	- // fit model
	- LauTimeDepFitModel* fitModel = new LauTimeDepFitModel(sigModelB0bar,sigModelB0,flavTag);
	-
	- std::vector<LauAbsCoeffSet*> coeffset;
	- coeffset.push_back( new LauRealImagCoeffSet("D*+_2", 1.00, 0.00, kTRUE, kTRUE) );
	- coeffset.push_back( new LauRealImagCoeffSet("D+_0", 0.53TMath::Cos( 3.00), 0.53*TMath::Sin( 3.00), kFALSE, kFALSE) );
	- coeffset.push_back( new LauRealImagCoeffSet("rho0(770)", 1.22TMath::Cos( 2.25), 1.22TMath::Sin( 2.25), kFALSE, kFALSE) );
	- coeffset.push_back( new LauRealImagCoeffSet("f_0(980)", 0.19TMath::Cos(-2.48), 0.19TMath::Sin(-2.48), kFALSE, kFALSE) );
	- coeffset.push_back( new LauRealImagCoeffSet("f_2(1270)", 0.75TMath::Cos( 2.97), 0.75TMath::Sin( 2.97), kFALSE, kFALSE) );
	-
	- for (std::vector<LauAbsCoeffSet*>::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) {
	- fitModel->setAmpCoeffSet(*iter);
	- }
	-
	- // background DP models
	- /*
	- LauBkgndDPModel* Bkgnd1Model = new LauBkgndDPModel(daughtersB0, vetoes);
	- LauBkgndDPModel* Bkgnd1Modelbar = new LauBkgndDPModel(daughtersB0bar, vetoes);
	- LauBkgndDPModel* Bkgnd2Model = new LauBkgndDPModel(daughtersB0, vetoes);
	- LauBkgndDPModel* Bkgnd2Modelbar = new LauBkgndDPModel(daughtersB0bar, vetoes);
	- fitModel->setBkgndDPModels( "Bkgnd1", Bkgnd1Model, Bkgnd1Modelbar );
	- fitModel->setBkgndDPModels( "Bkgnd2", Bkgnd2Model, Bkgnd2Modelbar );
	- */
	-
	- // 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 );
	-
	- // 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.0);
	- const Double_t maxDt(15.0);
	- const Double_t minDtErr(0.0);
	- const Double_t maxDtErr(0.15);
	-
	- LauParameter * tau = new LauParameter("dt_tau", 1.519, 0.5, 5.0, settings.fixLifetime);
	- LauParameter * freq = new LauParameter("dt_deltaM", 0.5064, 0.0, 1.0, settings.fixDeltaM);
	-
	- std::vector<LauAbsRValue*> dtPhysPars {
	- tau,
	- freq
	- };
	-
	- auto dtPhysModel = std::make_unique<LauDecayTimePhysicsModel>( LauDecayTime::FuncType::ExpTrig, dtPhysPars );
	-
	- const std::vector<bool> scale {
	- settings.perEventTimeErr && true,
	- settings.perEventTimeErr && true,
	- };
	- const std::size_t nGauss{scale.size()};
	-
	- LauParameter * mean0 = new LauParameter("dt_mean_0", scale[0] ? -2.01290e-03 : -2.25084e-03, -0.01, 0.01, kTRUE );
	- LauParameter * mean1 = new LauParameter("dt_mean_1", scale[1] ? -2.01290e-03 : -5.04275e-03, -0.01, 0.01, kTRUE );
	- LauParameter * sigma0 = new LauParameter("dt_sigma_0", scale[0] ? 9.95145e-01 : 3.03923e-02, 0.0, 2.0, kTRUE );
	- LauParameter * sigma1 = new LauParameter("dt_sigma_1", scale[1] ? 1.81715e+00 : 6.22376e-02, 0.0, 2.5, kTRUE );
	- LauParameter * frac1 = new LauParameter("dt_frac_1", scale[0] && scale[1] ? 1.-9.35940e-01 : 1.-7.69603e-01, 0.0, 1.0, kTRUE);
	-
	- std::vector<LauAbsRValue*> dtResoPars {
	- mean0,
	- mean1,
	- sigma0,
	- sigma1,
	- frac1
	- };
	-
	- auto dtResoModel = std::make_unique<LauDecayTimeResolution>( nGauss, dtResoPars, scale );
	-
	- // 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("dte-hist.root");
	- 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("dta-hist.root");
	- TH1* dtaHist = dynamic_cast<TH1*>(dtaFile->Get("dta_hist"));
	-
	- // 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};
	-
	- switch(settings.timeEffModel)
	- {
	- case LauDecayTime::EfficiencyMethod::Spline:
	- {
	- fitModel->setASqMaxValue(0.06);
	-
	- auto dtEffSpline = std::make_unique<Lau1DCubicSpline>( dtvals, effvals, Lau1DCubicSpline::SplineType::AkimaSpline );
	-
	- auto dtaModel = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Cubic>>( "dteff_QFS", std::move(dtEffSpline) );
	- 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 );
	- }
	-
	- if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) {
	- // For the QFS mode we just use the cubic model as it is
	- dtPdf->setSplineEfficiency( std::move(dtaModel) );
	- } else {
	- // For the CP modes we modify it using a corrective spline
	- auto dtCorrSpline = std::make_unique<Lau1DCubicSpline>( dtvals, corvals, Lau1DCubicSpline::SplineType::AkimaSpline );
	- auto dtaCPModel = std::make_unique<LauSplineDecayTimeEfficiency<LauSplineOrder::Sixth>>( "dteff_CP", *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 );
	-
	- //Background dt part
	- /*
	- TFile* background_dt = TFile::Open("Lifetimes_PV_WL.root");
	- TH1* bkgnd1DtHist = dynamic_cast<TH1*>( background_dt->Get("Bkgnd1") );
	- TH1* bkgnd2DtHist = dynamic_cast<TH1*>( background_dt->Get("Bkgnd2") );
	- TH1* bkgnd1DtErrHist = dynamic_cast<TH1*>( background_dt->Get("Bkgnd1_Err") );
	- TH1* bkgnd2DtErrHist = dynamic_cast<TH1*>( background_dt->Get("Bkgnd2_Err") );
	- LauDecayTimePdf* bkgnd1DtPdf{nullptr};
	- LauDecayTimePdf* bkgnd2DtPdf{nullptr};
	- if ( settings.timeResolution and settings.perEventTimeErr ) {
	- bkgnd1DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd1DtHist, bkgnd1DtErrHist );
	- bkgnd2DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd2DtHist, bkgnd2DtErrHist );
	- } else {
	- bkgnd1DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd1DtHist );
	- bkgnd2DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd2DtHist );
	- }
	- fitModel->setBkgndDtPdf("Bkgnd1",bkgnd1DtPdf);
	- fitModel->setBkgndDtPdf("Bkgnd2",bkgnd2DtPdf);
	- */
	-
	- // set the signal yield
	- Double_t nSigEvents{0};
	- switch (settings.dType) {
	- case LauTimeDepFitModel::CPEigenvalue::CPEven :
	- nSigEvents = 15000;
	- break;
	- case LauTimeDepFitModel::CPEigenvalue::CPOdd :
	- nSigEvents = 5000;
	- break;
	- case LauTimeDepFitModel::CPEigenvalue::QFS :
	- nSigEvents = 50000;
	- break;
	- }
	-
	- // set the background yields
	- /*
	- const Double_t nBkgnd1(100), nBkgnd2(200);
	- LauParameter* Bkgnd1Yield = new LauParameter("Bkgnd1",nBkgnd1,-1.0nBkgnd1,2.0nBkgnd1,kFALSE);
	- LauParameter* Bkgnd2Yield = new LauParameter("Bkgnd2",nBkgnd2,-1.0nBkgnd2,2.0nBkgnd2,kFALSE);
	- fitModel->setNBkgndEvents(Bkgnd1Yield);
	- fitModel->setNBkgndEvents(Bkgnd2Yield);
	- */
	-
	- std::cout<<"nSigEvents = "<<nSigEvents<<std::endl;
	- LauParameter* nSigPar = new LauParameter("signalEvents", nSigEvents, -2.0nSigEvents, 2.0nSigEvents, kTRUE);
	-
	- fitModel->setNSigEvents(nSigPar);
	-
	- // 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;
	-
	- 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 += ".root";
	-
	- if (settings.command == Command::Generate) {
	- rootFileName = "dummy.root";
	- tableFileName = "genResults";
	- } else {
	- rootFileName = "fit"; rootFileName += settings.iFit;
	- rootFileName += "_Results_"; rootFileName += dTypeStr;
	- rootFileName += "_expts"; rootFileName += settings.firstExptFit; rootFileName += "-"; rootFileName += settings.firstExptFit+settings.nExptFit-1;
	- rootFileName += ".root";
	-
	- tableFileName = "fit"; tableFileName += settings.iFit;
	- tableFileName += "_Results_"; tableFileName += dTypeStr;
	- tableFileName += "_expts"; tableFileName += settings.firstExptFit; tableFileName += "-"; tableFileName += settings.firstExptFit+settings.nExptFit-1;
	-
	- fitToyFileName = "fit"; fitToyFileName += settings.iFit;
	- fitToyFileName += "_ToyMC_"; fitToyFileName += dTypeStr;
	- fitToyFileName += "_expts"; fitToyFileName += settings.firstExptFit; fitToyFileName += "-"; fitToyFileName += settings.firstExptFit+settings.nExptFit-1;
	- fitToyFileName += ".root";
	-
	- splotFileName = "fit"; splotFileName += settings.iFit;
	- splotFileName += "_sPlot_"; splotFileName += dTypeStr;
	- splotFileName += "_expts"; splotFileName += settings.firstExptFit; splotFileName += "-"; splotFileName += settings.firstExptFit+settings.nExptFit-1;
	- splotFileName += ".root";
	- }
	-
	- // 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;
	-}
	diff --git a/examples/runTimeDepTest.sh b/examples/runTimeDepTest.sh
	index 508c0ee..99d2fdd 100755
	--- a/examples/runTimeDepTest.sh
	+++ b/examples/runTimeDepTest.sh
	@@ -1,242 +1,242 @@
	#!/bin/bash
	#SBATCH --partition=epp
	#SBATCH --ntasks=1
	#SBATCH --cpus-per-task=6
	#SBATCH --mem-per-cpu=3997
	#SBATCH --time=02:00:00
	#Run with sbatch runTimeDepTest.sh

	# Copyright 2020 University of Warwick
	#
	# Licensed under the Apache License, Version 2.0 (the "License");
	# you may not use this file except in compliance with the License.
	# You may obtain a copy of the License at
	#
	# http://www.apache.org/licenses/LICENSE-2.0
	#
	# Unless required by applicable law or agreed to in writing, software
	# distributed under the License is distributed on an "AS IS" BASIS,
	# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
	# See the License for the specific language governing permissions and
	# limitations under the License.

	# Laura++ package authors:
	# John Back
	# Paul Harrison
	# Thomas Latham

	processes=()

	# trap ^C and call ctrl_c()
	trap ctrl_c INT

	function ctrl_c()
	{
	echo " ** Caught ^C event!"
	for process in ${processes[@]}
	do
	if ps -u $USER \| grep -q "$process"
	then
	echo "Killing : $process"
	kill $process
	fi
	done
	exit 1
	}

	# Check that we have a kerberos token, since nothing will work without one!
	klist -s > /dev/null 2>&1
	if [ $? -ne 0 ]
	then
	echo "You don't have a valid Kerberos ticket!"
	exit 1
	fi

	iFit=0
	dta_model="spline"
	dtr=1
	dtr_perevent=0
	fixTau=1
	fixDeltaM=1
	seed=140279 # $(pwd \| xargs basename)
	dataFit=1
	runNo=12

	# Modify these if you want to run more experiments
	nExpt=1
	firstExpt=0

	logName_genQFS_run1=Logs/gen-QFS-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
	logName_genQFS_run2=Logs/gen-QFS-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
	logName_genCPE_KK_run1=Logs/gen-CPEvenKK-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
	logName_genCPE_KK_run2=Logs/gen-CPEvenKK-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
	logName_genCPE_pipi_run1=Logs/gen-CPEvenpipi-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log
	logName_genCPE_pipi_run2=Logs/gen-CPEvenpipi-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log

	logName_coord=Logs/coordinator-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
	logName_taskQFS_run1=Logs/task-QFS-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
	logName_taskQFS_run2=Logs/task-QFS-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
	logName_taskCPE_KK_run1=Logs/task-CPEvenKK-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
	logName_taskCPE_KK_run2=Logs/task-CPEvenKK-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
	logName_taskCPE_pipi_run1=Logs/task-CPEvenpipi-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log
	logName_taskCPE_pipi_run2=Logs/task-CPEvenpipi-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log

	mkdir -p Logs

	# make sure we have the most up-to-date json files
	eosArea="root://eoslhcb.cern.ch//eos/lhcb/wg/b2oc/B2D0pipi/timedep_fit_histos"
	rm -f run1DTAsplines.json run2DTAsplines.json Bd2D0pipi_DP_Model_Coeffs.json
	xrdcp $eosArea/run1DTAsplines.json .
	if [ $? -ne 0 ]
	then
	echo "Problem downloading json files from EOS, exiting"
	exit 1
	fi
	xrdcp $eosArea/run2DTAsplines.json .
	if [ $? -ne 0 ]
	then
	echo "Problem downloading json files from EOS, exiting"
	exit 1
	fi
	xrdcp $eosArea/Bd2D0pipi_DP_Model_Coeffs.json .
	if [ $? -ne 0 ]
	then
	echo "Problem downloading json files from EOS, exiting"
	exit 1
	fi

	if [ $iFit == 0 -a $dataFit == 0 ]
	then
	#run generation of the modes in parallel (put the tasks in the bg)
	echo "Generating samples..."
	if [ $runNo -eq 1 -o $runNo -eq 12 ]
	then
	- Test_Dpipi gen --dtype QFS --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_Kpi > $logName_genQFS_run1 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi gen --dtype QFS --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_Kpi > $logName_genQFS_run1 2>&1 &
	task=$!
	processes+=( $task )

	- Test_Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_KK > $logName_genCPE_KK_run1 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_KK > $logName_genCPE_KK_run1 2>&1 &
	task=$!
	processes+=( $task )

	- Test_Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_pipi > $logName_genCPE_pipi_run1 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_pipi > $logName_genCPE_pipi_run1 2>&1 &
	task=$!
	processes+=( $task )
	fi

	if [ $runNo -eq 2 -o $runNo -eq 12 ]
	then
	- Test_Dpipi gen --dtype QFS --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_Kpi > $logName_genQFS_run2 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi gen --dtype QFS --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_Kpi > $logName_genQFS_run2 2>&1 &
	task=$!
	processes+=( $task )

	- Test_Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_KK > $logName_genCPE_KK_run2 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_KK > $logName_genCPE_KK_run2 2>&1 &
	task=$!
	processes+=( $task )

	- Test_Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_pipi > $logName_genCPE_pipi_run2 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_pipi > $logName_genCPE_pipi_run2 2>&1 &
	task=$!
	processes+=( $task )
	fi

	wait #wait for the generation to complete
	fi

	echo "Running fit $iFit"
	if [ $runNo -eq 12 ]
	then
	nTasks=6
	else
	nTasks=3
	fi

	SimFitCoordinator $iFit $nExpt $firstExpt $nTasks > $logName_coord 2>&1 &
	task=$!
	echo "SimFitCoordinator process : $task"
	processes+=( $task )

	echo "Initialised coordinator process for fit $iFit"

	sleep 5

	port=""
	while [ "x$port" == "x" ]
	do
	sleep 2
	port=`tail $logName_coord \| grep "Waiting for connection" \| awk '{print $NF}'`
	done

	echo "Coordinator is listening on port $port"
	echo "Initialising tasks..."

	if [ $runNo -eq 1 -o $runNo -eq 12 ]
	then
	- Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run1 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run1 2>&1 &
	task=$!
	echo "QFS Run 1 process : $task"
	processes+=( $task )
	sleep 5
	- Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run1 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run1 2>&1 &
	task=$!
	echo "CP KK Run 1 process : $task"
	processes+=( $task )
	sleep 5
	- Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run1 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run1 2>&1 &
	task=$!
	echo "CP pipi Run 1 process : $task"
	processes+=( $task )
	fi
	if [ $runNo -eq 2 -o $runNo -eq 12 ]
	then
	- Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run2 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run2 2>&1 &
	task=$!
	echo "QFS Run 2 process : $task"
	processes+=( $task )
	sleep 5
	- Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run2 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run2 2>&1 &
	task=$!
	echo "CP KK Run 2 process : $task"
	processes+=( $task )
	sleep 5
	- Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run2 2>&1 &
	+ GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run2 2>&1 &
	task=$!
	echo "CP pipi Run 2 process : $task"
	processes+=( $task )
	fi

	wait #wait for the fits and coordinator to terminate

	echo "Fit completed"

	echo "Warnings/Errors from coordinator:"
	grep -e ERROR -e WARNING -e Error -e Warning $logName_coord

	echo "Warnings/Errors from CPeven KK tasks:"
	if [ $runNo -eq 1 -o $runNo -eq 12 ]
	then
	grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_KK_run1
	fi
	if [ $runNo -eq 2 -o $runNo -eq 12 ]
	then
	grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_KK_run2
	fi

	echo "Warnings/Errors from CPeven pipi tasks:"
	if [ $runNo -eq 1 -o $runNo -eq 12 ]
	then
	grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_pipi_run1
	fi
	if [ $runNo -eq 2 -o $runNo -eq 12 ]
	then
	grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_pipi_run2
	fi

	echo "Warnings/Errors from QFS tasks:"
	if [ $runNo -eq 1 -o $runNo -eq 12 ]
	then
	grep -e ERROR -e WARNING -e Error -e Warning $logName_taskQFS_run1
	fi
	if [ $runNo -eq 2 -o $runNo -eq 12 ]
	then
	grep -e ERROR -e WARNING -e Error -e Warning $logName_taskQFS_run2
	fi

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 5:48 AM (17 h, 52 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4947287
Default Alt Text: (87 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions