diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc
index f5dcd1f..3a45936 100755
--- a/examples/Test_Dpipi.cc
+++ b/examples/Test_Dpipi.cc
@@ -1,405 +1,407 @@
 
 #include <iostream>
 using std::cout;
 using std::cerr;
 using std::endl;
 
 #include <vector>
 #include <map>
 
 #include "TFile.h"
 #include "TH2.h"
 #include "TRandom.h"
 #include "TString.h"
 #include "TSystem.h"
 
 #include "LauDaughters.hh"
 #include "LauDecayTimePdf.hh"
 #include "LauEffModel.hh"
 #include "LauIsobarDynamics.hh"
 #include "LauMagPhaseCoeffSet.hh"
 #include "LauRandom.hh"
 #include "LauRealImagCoeffSet.hh"
 #include "LauTimeDepFitModel.hh"
 #include "LauVetoes.hh"
 #include "LauFlavTag.hh"
 #include "Lau1DHistPdf.hh"
 #include "Lau1DCubicSpline.hh"
 
 void usage(const TString& progName)
 {
 	cerr<<"Usage:"<<endl;
 	cerr<<progName<<" gen D_type [firstExptGen] [nExptGen] [CP eigenvalue]"<<endl;
 	cerr<<"or"<<endl;
 	cerr<<progName<<" fit D_type [port] [iFit] [firstExpt] [nExpt] [nExptGen]"<<endl;
 }
 
 int main(const int argc, const  char ** argv)
 {
 	const TString dtype = argv[2];
 	Int_t port = 0;
 	Int_t iFit = 0;
 	Int_t firstExpt = 0;
 	Int_t firstExptGen = 0;
 	Int_t nExpt = 1;
 	Int_t nExptGen = 1;
 	LauTimeDepFitModel::CPEigenvalue eigenvalue =  LauTimeDepFitModel::CPEven;
 
 	Bool_t fixPhiMix(kFALSE);
 	Bool_t useSinCos(kTRUE);
 
 	// check the command line arguments
 	if (argc<1) {
 		usage(argv[0]);
 		return EXIT_FAILURE;
 	}
 	TString command = argv[1];
 	if (command != "gen" &&  command != "fit") {
 		usage(argv[0]);
 		return EXIT_FAILURE;
 	}
 	if (command == "fit") {
 		if (argc>3) {
 			port = atoi(argv[3]);
 			if (argc>4) {
 				iFit = atoi(argv[4]);
 				if (argc>5) {
 					firstExpt = atoi(argv[5]);
 					if (argc>6) {
 						nExpt = atoi(argv[6]);
 						if (argc>7) {
 							nExptGen = atoi(argv[7]);
 						}
 					}
 				}
 			}
 		}
 
 		for (firstExptGen = 0; firstExptGen<(firstExpt+nExpt); firstExptGen+=nExptGen) {
 		}
 		firstExptGen -= nExptGen;
 
 		if ( (nExpt > nExptGen) || (nExptGen%nExpt != 0) ) {
 			cerr<<"Error, nExpt must be a factor of nExptGen."<<endl;
 			return EXIT_FAILURE;
 		}
 
 	} else {
 		if (argc>3) {
 			firstExptGen = atoi(argv[3]);
 			if (argc>4) {
 				nExptGen = atoi(argv[4]);
 				if (argc>5) {
 					Int_t eigval = atoi(argv[5]);
 					if ( eigval == 1 ) {
 						eigenvalue = LauTimeDepFitModel::CPOdd;
 					} else {
 						eigenvalue = LauTimeDepFitModel::CPEven;
 					}
 				}
 			}
 		}
 	}
 
 	cout<<"dtype "<<dtype<<" port "<<port<<" iFit "<<iFit<<" firstExpt "<<firstExpt<<" nExpt "<<nExpt<<endl;
 
 	Double_t nSigEvents(0);
 	if (dtype=="CPEven"){
 		nSigEvents = 15000;
 	} else {
 		nSigEvents = 10000;
 	}
 
 	LauDaughters* daughtersB0bar(0);
 	LauDaughters* daughtersB0(0);
 	LauEffModel* effModelB0bar(0);
 	LauEffModel* effModelB0(0);
 	LauVetoes* vetoes = new LauVetoes();
 	//vetoes->addMassVeto( 2, 2.00776, 2.01276 );
 
 	daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0");
 	daughtersB0    = new LauDaughters("B0",     "pi+", "pi-", "D0_bar");
 
 	// efficiency
 	effModelB0bar = new LauEffModel(daughtersB0bar, vetoes);
 	effModelB0 = new LauEffModel(daughtersB0, vetoes);
 
 	LauIsobarDynamics* sigModelB0bar(0);
 	LauIsobarDynamics* sigModelB0(0);
 	LauTimeDepFitModel* fitModel(0);
 
 	std::vector<LauParameter *> params;
 
 	LauFlavTag* flavTag = new LauFlavTag(params, kTRUE, "tagFlv", "tagCat");
 	flavTag->setTrueTagVarName("trueTag");
 	
 	// Use alternative tagging calibration parameters?
 	//flavTag->useAveOmegaDeltaOmega();
 
 	TFile* etafile(0);
 	TH1* etahist(0);
 	Lau1DHistPdf* etahistpdf(0);
-	if (command == "gen"){ 
+	//if (command == "gen"){ 
 		etafile = TFile::Open("/Users/mark/analysis/DpipiLaura/laura/examples/histogram.root");
 		etahist = dynamic_cast<TH1*>(etafile->Get("htemp"));
 		etahistpdf = new Lau1DHistPdf("eta",etahist,0,0.54,kFALSE,kFALSE);
-	}
+	//}
 
 	// signal dynamics
 	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);
 
 	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
 	fitModel = new LauTimeDepFitModel(sigModelB0bar,sigModelB0,flavTag);
 	fitModel->setASqMaxValue(1.45);
 
 	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.53*TMath::Cos( 3.00), 0.53*TMath::Sin( 3.00), kFALSE, kFALSE) );
 	coeffset.push_back( new LauRealImagCoeffSet("rho0(770)", 1.22*TMath::Cos( 2.25), 1.22*TMath::Sin( 2.25), kFALSE, kFALSE) );
 	coeffset.push_back( new LauRealImagCoeffSet("f_0(980)",  0.19*TMath::Cos(-2.48), 0.19*TMath::Sin(-2.48), kFALSE, kFALSE) );
 	coeffset.push_back( new LauRealImagCoeffSet("f_2(1270)", 0.75*TMath::Cos( 2.97), 0.75*TMath::Sin( 2.97), kFALSE, kFALSE) );
 
 	for (std::vector<LauAbsCoeffSet*>::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) {
 		fitModel->setAmpCoeffSet(*iter);
 	}
 
 	fitModel->setCPEigenvalue( eigenvalue );
 	fitModel->setPhiMix( 2.0*LauConstants::beta, fixPhiMix, useSinCos );
 
 	fitModel->setAsymmetries(0.0,kFALSE,1.0,kTRUE);
 
 	int tagCat(-1);
 	// Tag cat fractions, dilutions and Delta dilutions
 	if (dtype=="CPEven"){
 		flavTag->addValidTagCat(63);
 		flavTag->setSignalTagCatPars(63, 0.6, 1.0, 0.0, kTRUE, kTRUE,"PerEvtMistag",0.5,0.5,1.0,1.0,1.0);
 		tagCat=63;
 	} else {
 		// Still need tagCat 63 in QFS channels for calibration etc
 		flavTag->addValidTagCat(63);
 		flavTag->setSignalTagCatPars(63, 0.6, 1.0, 0.0, kTRUE, kTRUE,"PerEvtMistag",0.5,0.5,1.0,1.0,1.0);
 		tagCat=63;
 	}
-	if (command == "gen"){ 
+	//if (command == "gen"){ 
 		fitModel->setSignalFlavTagPdfs(63,etahistpdf);
-	}
+	//}
 
 
 	// Delta t PDFs
 
 	const Double_t minDt(0.0);
 	const Double_t maxDt(20.0);
 	const Double_t minDtErr(0.0);
 	const Double_t maxDtErr(2.5);
 	const Int_t nGauss(3);
 	std::vector<Bool_t> scale(nGauss);
 	scale[0] = kTRUE;
 	scale[1] = kTRUE;
 	scale[2] = kFALSE;
 
 	std::vector<LauAbsRValue*> dtPars(10);
 	std::vector<LauAbsRValue*> dtParsExp(9);
 	TString mean0Name("dt_");  mean0Name += "_mean_0";
 	TString mean1Name("dt_");  mean1Name += "_mean_1";
 	TString mean2Name("dt_");  mean2Name += "_mean_2";
 	TString sigma0Name("dt_"); sigma0Name += "_sigma_0";
 	TString sigma1Name("dt_"); sigma1Name += "_sigma_1";
 	TString sigma2Name("dt_"); sigma2Name += "_sigma_2";
 	TString frac1Name("dt_");  frac1Name += "_frac_1";
 	TString frac2Name("dt_");  frac2Name += "_frac_2";
 	TString tauName("dt_");    tauName += "_tau";
 	TString freqName("dt_");   freqName += "_deltaM";
 
 	LauParameter * mean0 = new LauParameter(mean0Name,  -0.181);
 	LauParameter * mean1 = new LauParameter(mean1Name,  -1.27);
 	LauParameter * mean2 = new LauParameter(mean2Name,  0.0);
 	LauParameter * sigma0 = new LauParameter(sigma0Name, 1.067);
 	LauParameter * sigma1 = new LauParameter(sigma1Name, 3.0);
 	LauParameter * sigma2 = new LauParameter(sigma2Name, 8.0);
 	LauParameter * frac1 = new LauParameter(frac1Name,  0.0930);
 	LauParameter * frac2 = new LauParameter(frac2Name,  0.0036);
 	LauParameter * tau = new LauParameter(tauName,    1.520);
 	LauParameter * freq = new LauParameter(freqName,   0.5064);
 
 	TString mean0tagcat63Name("dt_");  mean0tagcat63Name += 63;  mean0tagcat63Name += "_mean_0";
 	TString sigma0tagcat63Name("dt_"); sigma0tagcat63Name += 63; sigma0tagcat63Name += "_sigma_0";
 
 	LauParameter * mean0tagcat63 = new LauParameter(mean0tagcat63Name,  -0.031);
 	LauParameter * sigma0tagcat63 = new LauParameter(sigma0tagcat63Name, 0.972);
 
 	dtPars[0] = mean0;	dtParsExp[0] = mean0;
 	dtPars[1] = mean1;	dtParsExp[1] = mean1;
 	dtPars[2] = mean2;	dtParsExp[2] = mean2;
 	dtPars[3] = sigma0;	dtParsExp[3] = sigma0;
 	dtPars[4] = sigma1;	dtParsExp[4] = sigma1;
 	dtPars[5] = sigma2;	dtParsExp[5] = sigma2;
 	dtPars[6] = frac1;	dtParsExp[6] = frac1;
 	dtPars[7] = frac2;	dtParsExp[7] = frac2;
 	dtPars[8] = tau;	dtParsExp[8] = tau;
 	dtPars[9] = freq;
 
 	//Decay time acceptance spline - same for all tag cats (though doesn't have to be)
 	std::vector<Double_t> dtvals;
 	dtvals.push_back(0.0); dtvals.push_back(2.5); dtvals.push_back(5.0); dtvals.push_back(10.0); dtvals.push_back(20.0);
 	std::vector<Double_t> effvals;
-	effvals.push_back(0.05); effvals.push_back(0.3); effvals.push_back(0.9); effvals.push_back(0.95); effvals.push_back(0.99);
-	Lau1DCubicSpline* dtEffSpline = new Lau1DCubicSpline(dtvals,effvals);
+	effvals.push_back(0.0); effvals.push_back(0.3); effvals.push_back(0.9); effvals.push_back(0.95); effvals.push_back(1.0);
+	
+	//Lau1DCubicSpline* dtEffSpline = new Lau1DCubicSpline(dtvals,effvals);
+	Lau1DCubicSpline* dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::StandardSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural);
 
 	//LauParameters to float Spline y values
 	//Test by floating knot 2 for both tagcats
-	LauParameter* knot_0_0 = new LauParameter("knot_0_0",0.05,0.0,1.0,kTRUE);
+	LauParameter* knot_0_0 = new LauParameter("knot_0_0",0.0,0.0,1.0,kTRUE);
 	LauParameter* knot_0_1 = new LauParameter("knot_0_1",0.30,0.0,1.0,kTRUE);
 	LauParameter* knot_0_2 = new LauParameter("knot_0_2",0.90,0.0,1.0,kFALSE);
 	LauParameter* knot_0_3 = new LauParameter("knot_0_3",0.95,0.0,1.0,kTRUE);
-	LauParameter* knot_0_4 = new LauParameter("knot_0_4",0.99,0.0,1.0,kTRUE);
+	LauParameter* knot_0_4 = new LauParameter("knot_0_4",1.0,0.0,1.0,kTRUE);
 
-	LauParameter* knot_63_0 = new LauParameter("knot_63_0",0.05,0.0,1.0,kTRUE);
-	LauParameter* knot_63_1 = new LauParameter("knot_63_1",0.30,0.0,1.0,kTRUE);
+	LauParameter* knot_63_0 = new LauParameter("knot_63_0",0.0,0.0,1.0,kTRUE);
+	LauParameter* knot_63_1 = new LauParameter("knot_63_1",0.30,0.0,1.0,kFALSE);
 	LauParameter* knot_63_2 = new LauParameter("knot_63_2",0.90,0.0,1.0,kFALSE);
-	LauParameter* knot_63_3 = new LauParameter("knot_63_3",0.95,0.0,1.0,kTRUE);
-	LauParameter* knot_63_4 = new LauParameter("knot_63_4",0.99,0.0,1.0,kTRUE);
+	LauParameter* knot_63_3 = new LauParameter("knot_63_3",0.95,0.0,1.0,kFALSE);
+	LauParameter* knot_63_4 = new LauParameter("knot_63_4",1.0,0.0,1.0,kTRUE);
 
 	std::vector<LauParameter*> effPars0;
 	effPars0.reserve(5);
 	effPars0.push_back(knot_0_0); effPars0.push_back(knot_0_1); effPars0.push_back(knot_0_2); effPars0.push_back(knot_0_3); effPars0.push_back(knot_0_4);
 
 	std::vector<LauParameter*> effPars63;
 	effPars63.reserve(5);
 	effPars63.push_back(knot_63_0); effPars63.push_back(knot_63_1); effPars63.push_back(knot_63_2); effPars63.push_back(knot_63_3); effPars63.push_back(knot_63_4);
 	
 	if (dtype=="CPEven"){
 		LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::ExpTrig, nGauss, scale, LauDecayTimePdf::DecayTime );
 		//LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::SimFitSigBd, nGauss, scale, LauDecayTimePdf::DecayTime );
 		dtPdf->doSmearing(kFALSE);
 		dtPdf->setEffiSpline(dtEffSpline, effPars0);
 		fitModel->setSignalDtPdf( 0, dtPdf );
 	} else {
 		LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::ExpTrig, nGauss, scale, LauDecayTimePdf::DecayTime );
 		//LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::SimFitNormBd, nGauss, scale, LauDecayTimePdf::DecayTime );
 		dtPdf->doSmearing(kFALSE);
 		dtPdf->setEffiSpline(dtEffSpline, effPars0);
 		fitModel->setSignalDtPdf( 0, dtPdf );
 	}
 
 	//if (tagCat==63){
 		dtPars[0] = mean0tagcat63;
 		dtPars[1] = mean1->createClone();
 		dtPars[2] = mean2->createClone();
 		dtPars[3] = sigma0tagcat63;
 		dtPars[4] = sigma1->createClone();
 		dtPars[5] = sigma2->createClone();
 		dtPars[6] = frac1->createClone();
 		dtPars[7] = frac2->createClone();
 		dtPars[8] = tau->createClone();
 		dtPars[9] = freq->createClone();
 
 		//LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::ExpTrig, nGauss, scale, LauDecayTimePdf::DecayTime );
 		if (dtype=="CPEven"){
 			LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::SimFitSigBd, nGauss, scale, LauDecayTimePdf::DecayTime );
 			dtPdf->doSmearing(kFALSE);
 			dtPdf->setEffiSpline(dtEffSpline, effPars63);
 			fitModel->setSignalDtPdf( tagCat, dtPdf );
 		} else {
 			LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::SimFitNormBd, nGauss, scale, LauDecayTimePdf::DecayTime );
 			dtPdf->doSmearing(kFALSE);
 			dtPdf->setEffiSpline(dtEffSpline, effPars63);
 			fitModel->setSignalDtPdf( tagCat, dtPdf );
 		}
 
 	//}
 
 	// set the number of signal events
 	cout<<"nSigEvents = "<<nSigEvents<<endl;
 	LauParameter* nSigPar = new LauParameter("signalEvents", nSigEvents, -2.0*nSigEvents, 2.0*nSigEvents, kTRUE);
 
 	fitModel->setNSigEvents(nSigPar);
 
 	// set the number of experiments
 	if (command == "fit") {
 		fitModel->setNExpts(nExpt, firstExpt);
 	} else {
 		fitModel->setNExpts(nExptGen, firstExptGen);
 	}
 
 	fitModel->useAsymmFitErrors(kFALSE);
 	//fitModel->useRandomInitFitPars(kTRUE);
 	fitModel->useRandomInitFitPars(kFALSE);
 	fitModel->doPoissonSmearing(kFALSE);
 	fitModel->doEMLFit(kFALSE);
 	fitModel->writeLatexTable(kFALSE);
 
 
 	TString dataFile("");
 	TString treeName("fitTree");
 	TString rootFileName("");
 	TString tableFileName("");
 	TString fitToyFileName("");
 	TString splotFileName("");
 	if (command == "fit") {
 		dataFile = "TEST-Dpipi_"+dtype;
 		dataFile += "_expts"; dataFile += firstExptGen; dataFile += "-"; dataFile += firstExptGen+nExptGen-1;
 		dataFile += "_CP"; 
 		if ( eigenvalue == LauTimeDepFitModel::CPEven ) {
 			dataFile += "even";
 		} else {
 			dataFile += "odd";
 		}
 		dataFile += ".root";
 
 		rootFileName = "fits/fit"; rootFileName += iFit;
 		rootFileName += "_expts"; rootFileName += firstExpt; rootFileName += "-"; rootFileName += firstExpt+nExpt-1;
 		rootFileName += ".root";
 
 		tableFileName = "fitResults_"; tableFileName += iFit;
 		tableFileName += "_expts"; tableFileName += firstExpt; tableFileName += "-"; tableFileName += firstExpt+nExpt-1;
 
-		fitToyFileName = "fitToyMC_"; fitToyFileName += iFit;
+		fitToyFileName = "fitToyMC_"+dtype; fitToyFileName += iFit;
 		fitToyFileName += "_expts"; fitToyFileName += firstExpt; fitToyFileName += "-"; fitToyFileName += firstExpt+nExpt-1;
 		fitToyFileName += ".root";
 
 		splotFileName = "splot_"; splotFileName += iFit;
 		splotFileName += "_expts"; splotFileName += firstExpt; splotFileName += "-"; splotFileName += firstExpt+nExpt-1;
 		splotFileName += ".root";
 	} else {
 		dataFile = "TEST-Dpipi_"+dtype;
 		dataFile += "_expts"; dataFile += firstExptGen; dataFile += "-"; dataFile += firstExptGen+nExptGen-1;  dataFile += "_CP";
 		if ( eigenvalue == LauTimeDepFitModel::CPEven ) {
 			dataFile += "even";
 		} else {
 			dataFile += "odd";
 		}
 		dataFile += ".root";
 
 		rootFileName = "dummy.root";
 		tableFileName = "genResults";
 	}
 
 	// Generate toy from the fitted parameters
-	//fitModel->compareFitData(10, fitToyFileName);
+	fitModel->compareFitData(1, fitToyFileName);
 
 	// Write out per-event likelihoods and sWeights
 	//fitModel->writeSPlotData(splotFileName, "splot", kFALSE);
 
 	// Execute the generation/fit
 	if ( command == "fit" ){
 		fitModel->runSlave( dataFile, treeName, rootFileName, tableFileName, "localhost", port );
 	} else {
 		fitModel->run( command, dataFile, treeName, rootFileName, tableFileName );
 	}
 
 	return EXIT_SUCCESS;
 }
diff --git a/inc/Lau1DCubicSpline.hh b/inc/Lau1DCubicSpline.hh
index b3f3302..b1ee48a 100644
--- a/inc/Lau1DCubicSpline.hh
+++ b/inc/Lau1DCubicSpline.hh
@@ -1,212 +1,213 @@
 
 /*
 Copyright 2015 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
 /*! \file Lau1DCubicSpline.hh
     \brief File containing declaration of Lau1DCubicSpline class.
  */
 
 /*! \class Lau1DCubicSpline
     \brief Class for defining a 1D cubic spline based on a set of knots.
 
     Class for defining a 1D cubic spline based on a set of knots.
     Interpolation between the knots is performed either by one of 
     two types of cubic spline (standard or Akima) or by linear interpolation.
     The splines are defined by a piecewise cubic function which, between knots i and i+1, has the form
 
     f_i(x) = (1 - t)*y_i + t*y_i+1 + t*(1 - t)(a*(1 - t) + b*t)
 
     where t = (x - x_i)/(x_i+1 - x_i),
     a =  k_i  *(x_i+1 - x_i) - (y_i+1 - y_i),
     b = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i)
 
     and k_i is (by construction) the first derivative at knot i.
     f(x) and f'(x) are continuous at the internal knots by construction.
     
     For the standard splines, f''(x) is required to be continuous 
     at all internal knots placing n-2 constraints on the n parameters, k_i.
     The final two constraints are set by the boundary conditions.
     At each boundary, the function may be:
 
     (i)   Clamped :     f'(x) = C at the last knot
     (ii)  Natural :     f''(x) = 0 at the last knot
     (iii) Not a knot :  f'''(x) continuous at the second last knot
 
     The algorithms used in these splines can be found on:
     http://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline
     http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
 
     For the Akima splines, the values of k_i are determined from the slopes of the four nearest segments (a_i-1, a_i, a_i+1 and a_i+2) as
 
     k_i = ( | a_i+2 - a_i+1 | a_i  +  | a_i - a_i-1 | a_i+1 ) / ( | a_i+2 - a_i+1 |  +  | a_i - a_i-1 | )
 
     and as
 
     k_i = ( a_i + a_i+1 ) / 2
 
     in the special case a_i-1 == a_i != a_i+1 == a_i+2.
 
     Boundary conditions are specified by the relations
     a_2 - a_1 = a_1 - a_0 = a_0 - a_-1 and 
     a_n-1 - a_n-2 = a_n - a_n-1 = a_n+1 - a_n
 
     The algorithms used in these splines can be found in:
     J.ACM vol. 17 no. 4 pp 589-602
  */
 
 #ifndef LAU_1DCUBICSPLINE
 #define LAU_1DCUBICSPLINE
 
 #include <vector>
-
+#include "LauParameter.hh"
 #include "Rtypes.h"
 
 class Lau1DCubicSpline {
 
 	public:
 		//! Define the allowed interpolation types
 		enum LauSplineType {
 			StandardSpline, /*!< standard cubic splines with f''(x) continuous at all internal knots */
 			AkimaSpline, /*!< Akima cubic splines with f'(x) at each knot defined locally by the positions of only five knots */
 			LinearInterpolation /*! Linear interpolation between each pair of knots */
 		};
 
 		//! Define the allowed boundary condition types
 		/*!
 		    These are only supported by standard splines
 		*/
 		enum LauSplineBoundaryType {
 			Clamped, /*!< clamped boundary - f'(x) = C */
 			Natural, /*!< natural boundary - f''(x) = 0 */
 			NotAKnot /*!< 'not a knot' boundary - f'''(x) continuous at second last knot */
 		};
 
 		//! Constructor
 		/*!
 		    /param [in] xs the x-values of the knots
 		    /param [in] ys the y-values of the knots
 		    /param [in] leftBound the left-hand boundary condition
 		    /param [in] rightBound the right-hand boundary condition
 		    /param [in] dydx0 the gradient at the left-hand boundary of a clamped spline
 		    /param [in] dydxn the gradient at the right-hand boundary of a clamped spline
 		*/
 		Lau1DCubicSpline(const std::vector<Double_t>& xs, const std::vector<Double_t>& ys,
 				 LauSplineType type = Lau1DCubicSpline::StandardSpline,
 				 LauSplineBoundaryType leftBound = Lau1DCubicSpline::NotAKnot, 
 				 LauSplineBoundaryType rightBound = Lau1DCubicSpline::NotAKnot,
 				 Double_t dydx0 = 0.0, Double_t dydxn = 0.0);
 
 		//! Destructor
 		virtual ~Lau1DCubicSpline();
 
 		//! Evaluate the function at given point
 		/*!
 		    \param [in] x the x-coordinate
 		    \return the value of the spline at x
 		*/
 		Double_t evaluate(Double_t x) const;
 
 		//! Update the y-values of the knots
 		/*!
 		    \param [in] ys the y-values of the knots
 		*/
 		void updateYValues(const std::vector<Double_t>& ys);
+		void updateYValues(const std::vector<LauParameter*>& ys);
 
 		//! Update the type of interpolation to perform
 		/*!
 		    \param [in] type the type of interpolation
 		*/
 		void updateType(LauSplineType type);
 
 		//! Update the boundary conditions for the spline
 		/*!
 		    /param [in] leftBound the left-hand boundary condition
 		    /param [in] rightBound the right-hand boundary condition
 		    /param [in] dydx0 the gradient at the left-hand boundary of a clamped spline
 		    /param [in] dydxn the gradient at the right-hand boundary of a clamped spline
 		*/
 		void updateBoundaryConditions(LauSplineBoundaryType leftBound, 
 					      LauSplineBoundaryType rightBound, 
 					      Double_t dydx0 = 0.0,
 					      Double_t dydxn = 0.0);
 
 		//! Return the number of knots
 		UInt_t getnKnots() const {return nKnots_;}
 
 		//! Get y values
 		std::vector<Double_t> getYValues(){return y_;}
 
 	private:
 		//! Copy constructor - not implemented
 		Lau1DCubicSpline( const Lau1DCubicSpline& rhs );
 
 		//! Copy assignment operator - not implemented
 		Lau1DCubicSpline& operator=(const Lau1DCubicSpline& rhs);
 
 		//! Initialise the class
 		void init();
 
 		//! Calculate the first derivative at each knot
 		void calcDerivatives();
 
 		//! Calculate the first derivatives according to the standard method
 		void calcDerivativesStandard();
 		//! Calculate the first derivatives according to the Akima method
 		void calcDerivativesAkima();
 
 		//! The number of knots in the spline
 		const UInt_t nKnots_;
 
 		//! The x-value at each knot
 		std::vector<Double_t> x_;
 		//! The y-value at each knot
 		std::vector<Double_t> y_;
 		//! The first derivative at each knot
 		std::vector<Double_t> dydx_;
 
 		//! The 'a' coefficients used to determine the derivatives
 		std::vector<Double_t> a_;
 		//! The 'b' coefficients used to determine the derivatives
 		std::vector<Double_t> b_;
 		//! The 'c' coefficients used to determine the derivatives
 		std::vector<Double_t> c_;
 		//! The 'd' coefficients used to determine the derivatives
 		std::vector<Double_t> d_;
 
 		//! The type of interpolation to be performed
 		LauSplineType type_;
 
 		//! The left-hand boundary condition on the spline
 		LauSplineBoundaryType leftBound_;
 		//! The right-hand boundary condition on the spline
 		LauSplineBoundaryType rightBound_;
 
 		//! The gradient at the left boundary for a clamped spline
 		Double_t dydx0_;
 		//! The gradient at the right boundary for a clamped spline
 		Double_t dydxn_;
 
 		ClassDef(Lau1DCubicSpline, 0); // Class for defining a 1D cubic spline
 };
 
 #endif
diff --git a/src/Lau1DCubicSpline.cc b/src/Lau1DCubicSpline.cc
index c1eb757..7dc0a0e 100644
--- a/src/Lau1DCubicSpline.cc
+++ b/src/Lau1DCubicSpline.cc
@@ -1,343 +1,351 @@
 
 /*
 Copyright 2015 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
 /*! \file Lau1DCubicSpline.cc
   \brief File containing implementation of Lau1DCubicSpline class.
  */
 
 #include <cmath>
 #include <cstdlib>
 #include <iostream>
 
 #include <TH2.h>
 #include <TSystem.h>
 
 #include "Lau1DCubicSpline.hh"
 
 ClassImp(Lau1DCubicSpline)
 
 
 Lau1DCubicSpline::Lau1DCubicSpline(const std::vector<Double_t>& xs, const std::vector<Double_t>& ys,
 				   LauSplineType type, 
 				   LauSplineBoundaryType leftBound, 
 				   LauSplineBoundaryType rightBound,
 				   Double_t dydx0, Double_t dydxn) :
 	nKnots_(xs.size()),
 	x_(xs),
 	y_(ys),
 	type_(type),
 	leftBound_(leftBound),
 	rightBound_(rightBound),
 	dydx0_(dydx0),
 	dydxn_(dydxn)
 {
 	init();
 }
 
 Lau1DCubicSpline::~Lau1DCubicSpline()
 {
 }
 
 Double_t Lau1DCubicSpline::evaluate(Double_t x) const
 {
 	// do not attempt to extrapolate the spline
 	if( x<x_[0] || x>x_[nKnots_-1] ) {
 		std::cout << "WARNING in Lau1DCubicSpline::evaluate : function is only defined between " << x_[0] << " and " << x_[nKnots_-1] << std::endl;
 		std::cout << "                                        value at " << x << " returned as 0" << std::endl;
 		return 0.;
 	}
 
 	// first determine which 'cell' of the spline x is in
 	// cell i runs from knot i to knot i+1
 	Int_t cell(0);
 
 	while( x > x_[cell+1] ) {
 		++cell;
 	}
 
 	// obtain x- and y-values of the neighbouring knots
 	Double_t xLow  = x_[cell];
 	Double_t xHigh = x_[cell+1];
 	Double_t yLow  = y_[cell];
 	Double_t yHigh = y_[cell+1];
 
 	if(type_ == Lau1DCubicSpline::LinearInterpolation) {
 		return yHigh*(x-xLow)/(xHigh-xLow) + yLow*(xHigh-x)/(xHigh-xLow);
 	}
 
 	// obtain t, the normalised x-coordinate within the cell, 
 	// and the coefficients a and b, which are defined in cell i as:
 	//
 	//      a_i =  k_i  *(x_i+1 - x_i) - (y_i+1 - y_i),
 	//      b_i = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i)
 	//
 	// where k_i is (by construction) the first derivative at knot i
 	Double_t t = (x - xLow) / (xHigh - xLow);
 	Double_t a =     dydx_[cell]   * (xHigh - xLow) - (yHigh - yLow);
 	Double_t b = -1.*dydx_[cell+1] * (xHigh - xLow) + (yHigh - yLow);
 
 	Double_t retVal = (1 - t) * yLow + t * yHigh + t * (1 - t) * ( a * (1 - t) + b * t );
 
 	return retVal;
 }
 
 void Lau1DCubicSpline::updateYValues(const std::vector<Double_t>& ys)
 {
 	y_ = ys;
 	this->calcDerivatives();
 }
 
+void Lau1DCubicSpline::updateYValues(const std::vector<LauParameter*>& ys)
+{
+	for (Int_t i=0; i<nKnots_; ++i){
+		y_[i] = ys[i]->unblindValue();
+	}
+	this->calcDerivatives();
+}
+
 void Lau1DCubicSpline::updateType(LauSplineType type)
 {
 	if(type_ != type) {
 		type_ = type;
 		this->calcDerivatives();
 	}
 }
 
 void Lau1DCubicSpline::updateBoundaryConditions(LauSplineBoundaryType leftBound,
 		LauSplineBoundaryType rightBound,
 		Double_t dydx0, Double_t dydxn)
 {
 	Bool_t updateDerivatives(kFALSE);
 
 	if(leftBound_ != leftBound || rightBound_ != rightBound) {
 		leftBound_ = leftBound;
 		rightBound_ = rightBound;
 		updateDerivatives = kTRUE;
 	}
 
 	if(dydx0_ != dydx0) {
 		dydx0_ = dydx0;
 		if(leftBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = kTRUE;
 	}
 
 	if(dydxn_ != dydxn) {
 		dydxn_ = dydxn;
 		if(rightBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = kTRUE;
 	}
 
 	if(updateDerivatives) {
 		this->calcDerivatives();
 	}
 }
 
 void Lau1DCubicSpline::init()
 {
 	if( y_.size() != x_.size()) {
 		std::cout << "ERROR in Lau1DCubicSpline::init : The number of y-values given does not match the number of x-values" << std::endl;
 		std::cout << "                                  Found " <<  y_.size() << ", expected " << x_.size() << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	dydx_.insert(dydx_.begin(),nKnots_,0.);
 
 	a_.insert(a_.begin(),nKnots_,0.);
 	b_.insert(b_.begin(),nKnots_,0.);
 	c_.insert(c_.begin(),nKnots_,0.);
 	d_.insert(d_.begin(),nKnots_,0.);
 
 	this->calcDerivatives();
 }
 
 void Lau1DCubicSpline::calcDerivatives()
 {
 	switch ( type_ ) {
 		case Lau1DCubicSpline::StandardSpline :
 			this->calcDerivativesStandard();
 			break;
 		case Lau1DCubicSpline::AkimaSpline :
 			this->calcDerivativesAkima();
 			break;
 		case Lau1DCubicSpline::LinearInterpolation :
 			//derivatives not needed for linear interpolation
 			break;
 	}
 }
 
 void Lau1DCubicSpline::calcDerivativesStandard()
 {
 	// derivatives are determined such that the second derivative is continuous at internal knots
 
 	// derivatives, k_i, are the solutions to a set of linear equations of the form:
 	// a_i * k_i-1  +  b_i * k+i  +  c_i * k_i+1  =  d_i with a_0 = 0, c_n-1 = 0
 
 	// this is solved using the tridiagonal matrix algorithm as on en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
 
 	// first and last equations give boundary conditions
 	// - for natural boundary, require f''(x) = 0 at end knot
 	// - for 'not a knot' boundary, require f'''(x) continuous at second knot
 	// - for clamped boundary, require predefined value of f'(x) at end knot
 
 	// non-zero values of a_0 and c_n-1 would give cyclic boundary conditions
 	a_[0] = 0.;
 	c_[nKnots_-1] = 0.;
 
 	// set left boundary condition
 	if(leftBound_ == Lau1DCubicSpline::Natural) {
 		b_[0] = 2./(x_[1]-x_[0]);
 		c_[0] = 1./(x_[1]-x_[0]);
 		d_[0] = 3.*(y_[1]-y_[0])/((x_[1]-x_[0])*(x_[1]-x_[0]));
 
 	} else if(leftBound_ == Lau1DCubicSpline::NotAKnot) {
 		// define the width, h, and the 'slope', delta, of the first cell
 		Double_t h1(x_[1]-x_[0]), h2(x_[2]-x_[0]);
 		Double_t delta1((y_[1]-y_[0])/h1), delta2((y_[2]-y_[1])/h2);
 
 		// these coefficients can be determined by requiring f'''_0(x_1) = f'''_1(x_1)
 		// the requirement f''_0(x_1) = f''_1(x_1) has been used to remove the dependence on k_2
 		b_[0] = h2;
 		c_[0] = h1+h2;
 		d_[0] = delta1*(2.*h2*h2 + 3.*h1*h2)/(h1+h2) + delta2*5.*h1*h1/(h1+h2);
 
 	} else { //Clamped
 		b_[0] = 1.;
 		c_[0] = 0.;
 		d_[0] = dydx0_;
 	}
 
 	// set right boundary condition
 	if(rightBound_ == Lau1DCubicSpline::Natural) {
 		a_[nKnots_-1] = 1./(x_[nKnots_-1]-x_[nKnots_-2]);
 		b_[nKnots_-1] = 2./(x_[nKnots_-1]-x_[nKnots_-2]);
 		d_[nKnots_-1] = 3.*(y_[nKnots_-1]-y_[nKnots_-2])/((x_[nKnots_-1]-x_[nKnots_-2])*(x_[nKnots_-1]-x_[nKnots_-2]));
 
 	} else if(rightBound_ == Lau1DCubicSpline::NotAKnot) {
 		// define the width, h, and the 'slope', delta, of the last cell
 		Double_t hnm1(x_[nKnots_-1]-x_[nKnots_-2]), hnm2(x_[nKnots_-2]-x_[nKnots_-3]);
 		Double_t deltanm1((y_[nKnots_-1]-y_[nKnots_-2])/hnm1), deltanm2((y_[nKnots_-2]-y_[nKnots_-3])/hnm2);
 
 		// these coefficients can be determined by requiring f'''_n-3(x_n-2) = f'''_n-2(x_n-2)
 		// the requirement f''_n-3(x_n-2) = f''_n-2(x_n-2) has been used to remove 
 		// the dependence on k_n-3
 		a_[nKnots_-1] = hnm2 + hnm1;
 		b_[nKnots_-1] = hnm1;
 		d_[nKnots_-1] = deltanm2*hnm1*hnm1/(hnm2+hnm1) + deltanm1*(2.*hnm2*hnm2 + 3.*hnm2*hnm1)/(hnm2+hnm1);
 
 	} else { //Clamped
 		a_[nKnots_-1] = 0.;
 		b_[nKnots_-1] = 1.;
 		d_[nKnots_-1] = dydxn_;
 	}
 
 	// the remaining equations ensure that f_i-1''(x_i) = f''_i(x_i) for all internal knots
 	for(UInt_t i=1; i<nKnots_-1; ++i) {
 		a_[i] = 1./(x_[i]-x_[i-1]);
 		b_[i] = 2./(x_[i]-x_[i-1]) + 2./(x_[i+1]-x_[i]);
 		c_[i] =                      1./(x_[i+1]-x_[i]);
 		d_[i] = 3.*(y_[i]-y_[i-1])/((x_[i]-x_[i-1])*(x_[i]-x_[i-1])) + 3.*(y_[i+1]-y_[i])/((x_[i+1]-x_[i])*(x_[i+1]-x_[i]));
 	}
 
 	// forward sweep - replace c_i and d_i with
 	//
 	// c'_i =  c_i /  b_i                 for i = 0
 	//         c_i / (b_i - a_i * c_i-1)  otherwise
 	//
 	// and 
 	//
 	// d'_i =   d_i                /  b_i                 for i = 0
 	//         (d_i - a_i * d_i-1) / (b_i - a_i * c_i-1)  otherwise
 
 	c_[0] /= b_[0];
 	d_[0] /= b_[0];
 
 	for(UInt_t i=1; i<nKnots_-1; ++i) {
 		c_[i] =  c_[i]                  / (b_[i] - a_[i]*c_[i-1]);
 		d_[i] = (d_[i] - a_[i]*d_[i-1]) / (b_[i] - a_[i]*c_[i-1]);
 	}
 
 	d_[nKnots_-1] = (d_[nKnots_-1] - a_[nKnots_-1]*d_[nKnots_-2]) / (b_[nKnots_-1] - a_[nKnots_-1]*c_[nKnots_-2]);
 
 	// back substitution - calculate k_i = dy/dx at each knot
 	//
 	// k_i = d'_i                   for i = n-1
 	//       d'_i - c'_i * k_i+1    otherwise
 
 	dydx_[nKnots_-1] = d_[nKnots_-1];
 
 	for(Int_t i=nKnots_-2; i>=0; --i) {
 		dydx_[i] = d_[i] - c_[i]*dydx_[i+1];
 	}
 }
 
 void Lau1DCubicSpline::calcDerivativesAkima() {
 	//derivatives are calculated according to the Akima method
 	// J.ACM vol. 17 no. 4 pp 589-602
 
 	Double_t am1(0.), an(0.), anp1(0.);
 
 	// a[i] is the slope of the segment from point i-1 to point i
 	//
 	// n.b. segment 0 is before point 0 and segment n is after point n-1
 	//      internal segments are numbered 1 - n-1
 	for(UInt_t i=1; i<nKnots_; ++i) {
 		a_[i] = (y_[i]-y_[i-1])/(x_[i]-x_[i-1]);
 	}
 
 	// calculate slopes between additional points on each end according to the Akima method
 	// method assumes that the additional points follow a quadratic defined by the last three points
 	// this leads to the relations a[2] - a[1] = a[1] - a[0] = a[0] - a[-1] and a[n-1] - a[n-2] = a[n] - a[n-1] = a[n+1] - a[n]
 	a_[0]         = 2*a_[1]         - a_[2];
 	am1           = 2*a_[0]         - a_[1];
 
 	an            = 2*a_[nKnots_-1] - a_[nKnots_-2];
 	anp1          = 2*an            - a_[nKnots_-1];
 
 	// b[i] is the weight of a[i]   towards dydx[i]
 	// c[i] is the weight of a[i+1] towards dydx[i]
 	// See Appendix A of J.ACM vol. 17 no. 4 pp 589-602 for a justification of these weights
 	b_[0] = TMath::Abs(a_[1] - a_[0]);
 	b_[1] = TMath::Abs(a_[2] - a_[1]);
 
 	c_[0] = TMath::Abs(a_[0] - am1  );
 	c_[1] = TMath::Abs(a_[1] - a_[0]);
 
 	for(UInt_t i=2; i<nKnots_-2; ++i) {
 		b_[i] = TMath::Abs(a_[i+2] - a_[i+1]);
 		c_[i] = TMath::Abs(a_[i]   - a_[i-1]);
 	}
 
 	b_[nKnots_-2] = TMath::Abs(an            - a_[nKnots_-1]);
 	b_[nKnots_-1] = TMath::Abs(anp1          - an           );
 
 	c_[nKnots_-2] = TMath::Abs(a_[nKnots_-2] - a_[nKnots_-3]);
 	c_[nKnots_-1] = TMath::Abs(a_[nKnots_-1] - a_[nKnots_-2]);
 
 	// dy/dx calculated as the weighted average of a[i] and a[i+1]:
 	// dy/dx_i = ( | a_i+2 - a_i+1 | a_i  +  | a_i - a_i-1 | a_i+1 ) / ( | a_i+2 - a_i+1 |  +  | a_i - a_i-1 | )
 	// in the special case a_i-1 == a_i != a_i+1 == a_i+2 this function is undefined so dy/dx is then defined as (a_i + a_i+1) / 2
 	for(UInt_t i=0; i<nKnots_-2; ++i) {
 		if(b_[i]==0 && c_[i]==0) {
 			dydx_[i] = ( a_[i] + a_[i+1] ) / 2.;
 		} else {
 			dydx_[i] = ( b_[i] * a_[i] + c_[i] * a_[i+1] ) / ( b_[i] + c_[i] );
 		}
 	}
 
 	if(b_[nKnots_-1]==0 && c_[nKnots_-1]==0) {
 		dydx_[nKnots_-1] = ( a_[nKnots_-1] + an ) / 2.;
 	} else {
 		dydx_[nKnots_-1] = ( b_[nKnots_-1] * a_[nKnots_-1] + c_[nKnots_-1] * an ) / ( b_[nKnots_-1] + c_[nKnots_-1] );
 	}
 }
diff --git a/src/LauDecayTimePdf.cc b/src/LauDecayTimePdf.cc
index 7581104..c80227a 100644
--- a/src/LauDecayTimePdf.cc
+++ b/src/LauDecayTimePdf.cc
@@ -1,1563 +1,1560 @@
 
 /*
 Copyright 2006 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
 /*! \file LauDecayTimePdf.cc
     \brief File containing implementation of LauDecayTimePdf class.
 */
 
 #include <iostream>
 #include <vector>
 using std::cout;
 using std::cerr;
 using std::endl;
 
 #include <complex>
 using std::complex;
 
 #include "TMath.h"
 #include "TRandom.h"
 #include "TSystem.h"
 
 #include "Lau1DCubicSpline.hh"
 #include "Lau1DHistPdf.hh"
 #include "LauConstants.hh"
 #include "LauComplex.hh"
 #include "LauDecayTimePdf.hh"
 #include "LauFitDataTree.hh"
 #include "LauParameter.hh"
 #include "LauRandom.hh"
 
 ClassImp(LauDecayTimePdf)
 
 
 LauDecayTimePdf::LauDecayTimePdf(const TString& theVarName, const TString& theVarErrName, const std::vector<LauAbsRValue*>& params,
         Double_t minAbscissaVal, Double_t maxAbscissaVal,
 	Double_t minAbscissaErr, Double_t maxAbscissaErr,
 	FuncType type, UInt_t nGauss, const std::vector<Bool_t>& scale, TimeMeasurementMethod method) :
 	varName_(theVarName),
 	varErrName_(theVarErrName),
 	param_(params),
 	smear_(kTRUE),
 	minAbscissa_(minAbscissaVal),
 	maxAbscissa_(maxAbscissaVal),
 	minAbscissaError_(minAbscissaErr),
 	maxAbscissaError_(maxAbscissaErr),
 	abscissaError_(0.0),
 	abscissaErrorGenerated_(kFALSE),
 	errorDistMPV_(0.230),  // for signal 0.234, for qqbar 0.286
 	errorDistSigma_(0.075),  // for signal 0.073, for qqbar 0.102
 	nGauss_(nGauss),
 	mean_(nGauss_,0),
 	sigma_(nGauss_,0),
 	frac_(nGauss_-1,0),
 	tau_(0),
 	deltaM_(0),
 	deltaGamma_(0),
 	fracPrompt_(0),
 	type_(type),
 	method_(method),
 	scaleMeans_(scale),
 	scaleWidths_(scale),
 	expTerm_(0.0),
 	cosTerm_(0.0),
 	sinTerm_(0.0),
 	coshTerm_(0.0),
 	sinhTerm_(0.0),
         normTermExp_(0.0),
         normTermCosh_(0.0),
         normTermSinh_(0.0),
         errTerm_(0.0),
         pdfTerm_(0.0),
         effiTerm_(0.0),
 	state_(Good),
         errHist_(0),
         pdfHist_(0),
         effiFun_(0),
 	effiPars_(0),
 	AProd_(0),
 	ARaw_(0),
 	flavTag_(0),
 	asymTerm_(0),
 	omegaTerm_(0),
 	omegaTermh_(0)
 {
 	this->initialise();
 
 	// Calculate the integrals of the decay time independent of the t
 	this->storeIntegralTimeDep(minAbscissa_, maxAbscissa_);
 }
 
 LauDecayTimePdf::LauDecayTimePdf(const TString& theVarName, const TString& theVarErrName, const std::vector<LauAbsRValue*>& params,
 		Double_t minAbscissaVal, Double_t maxAbscissaVal,
 		Double_t minAbscissaErr, Double_t maxAbscissaErr,
 		FuncType type, UInt_t nGauss, const std::vector<Bool_t>& scaleMeans, const std::vector<Bool_t>& scaleWidths, TimeMeasurementMethod method) :
 	varName_(theVarName),
 	varErrName_(theVarErrName),
 	param_(params),
 	smear_(kTRUE),
 	minAbscissa_(minAbscissaVal),
 	maxAbscissa_(maxAbscissaVal),
 	minAbscissaError_(minAbscissaErr),
 	maxAbscissaError_(maxAbscissaErr),
 	abscissaError_(0.0),
 	abscissaErrorGenerated_(kFALSE),
 	errorDistMPV_(0.230),  // for signal 0.234, for qqbar 0.286
 	errorDistSigma_(0.075),  // for signal 0.073, for qqbar 0.102
 	nGauss_(nGauss),
 	mean_(nGauss_,0),
 	sigma_(nGauss_,0),
 	frac_(nGauss_-1,0),
 	tau_(0),
 	deltaM_(0),
 	deltaGamma_(0),
 	fracPrompt_(0),
 	type_(type),
 	method_(method),
 	scaleMeans_(scaleMeans),
 	scaleWidths_(scaleWidths),
 	expTerm_(0.0),
 	cosTerm_(0.0),
 	sinTerm_(0.0),
 	coshTerm_(0.0),
 	sinhTerm_(0.0),
         normTermExp_(0.0),
         normTermCosh_(0.0),
         normTermSinh_(0.0),
         errTerm_(0.0),
         pdfTerm_(0.0),
         effiTerm_(0.0),
 	state_(Good),
 	errHist_(0),
         pdfHist_(0),
         effiFun_(0),
 	effiPars_(0),
 	AProd_(0),
 	ARaw_(0),
 	flavTag_(0),
 	asymTerm_(0),
 	omegaTerm_(0),
 	omegaTermh_(0)
 {
 	this->initialise();
 
 	// Calculate the integrals of the decay time independent of the t
 	this->storeIntegralTimeDep(minAbscissa_, maxAbscissa_);
 }
 
 LauDecayTimePdf::~LauDecayTimePdf()
 {
 	// Destructor
 	delete errHist_;  errHist_ = 0;
         delete pdfHist_;  pdfHist_ = 0;
         delete effiFun_;  effiFun_ = 0;
 }
 
 void LauDecayTimePdf::initialise()
 {
 	// The parameters are:
 	// - the mean and the sigma (bias and spread in resolution) of the gaussian(s)
 	// - the mean lifetime, denoted tau, of the exponential decay
 	// - the frequency of oscillation, denoted Delta m, of the cosine and sine terms
 	// - the decay width difference, denoted Delta Gamma, of the hyperbolic cosine and sine terms
 	//
 	// The next two arguments specify the range in which the PDF is defined,
 	// and the PDF will be normalised w.r.t. these limits.
 	//
 	// The final three arguments define the type of Delta t PDF (Delta, Exp, ExpTrig or ExpHypTrig ), the number of gaussians
 	// and whether or not the gaussian parameters should be scaled by the per-event errors on Delta t
 
 	// First check whether the scale vector is nGauss in size
 	if (nGauss_ != scaleMeans_.size() || nGauss_ != scaleWidths_.size()) {
 		cerr<<"ERROR in LauDecayTimePdf::initialise : scale vector size not the same as nGauss."<<endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
         if (type_ == Hist){
 	  if (this->nParameters() != 0){
 	    cerr<<"ERROR in LauDecayTimePdf::initialise : Hist PDF should have 0 parameters"<<endl;
 	    gSystem->Exit(EXIT_FAILURE);
 	  }
         }else{
 	  TString meanName("mean_");
 	  TString sigmaName("sigma_");
 	  TString fracName("frac_");
 	  Bool_t foundParams(kTRUE);
 
 	  for (UInt_t i(0); i<nGauss_; ++i) {
 	    TString tempName(meanName); tempName += i;
 	    TString tempName2(sigmaName); tempName2 += i;
 	    TString tempName3(fracName); tempName3 += i;
 	    mean_[i] = this->findParameter(tempName);
 	    foundParams &= (mean_[i] != 0);
 	    sigma_[i] = this->findParameter(tempName2);
 	    foundParams &= (sigma_[i] != 0);
 	    if (i!=0) {
 	      frac_[i-1] = this->findParameter(tempName3);
 	      foundParams &= (frac_[i-1] != 0);
 	    }
 	  }
 
 	  if (type_ == Delta) {
 	    if ((this->nParameters() != (3*nGauss_-1)) || (!foundParams)) {
 	      cerr<<"ERROR in LauDecayTimePdf::initialise : Delta type PDF requires:"<<endl;
 	      cerr<<"	                                     - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\""<<endl;
 	      cerr<<"	                                     - nGauss-1 fractions: \"frac_i\", where i!=0"<<endl;
 	      gSystem->Exit(EXIT_FAILURE);
 	    }
 	  } else if (type_ == Exp) {
 	    tau_ = this->findParameter("tau");
 	    foundParams &= (tau_ != 0);
 	    if ((this->nParameters() != (3*nGauss_-1+1)) || (!foundParams)) {
 	      cerr<<"ERROR in LauDecayTimePdf::initialise : Exp type PDF requires:"<<endl;
 	      cerr<<"	                                     - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\""<<endl;
 	      cerr<<"	                                     - nGauss-1 fractions: \"frac_i\", where i!=0"<<endl;
 	      cerr<<"                                      - the lifetime of the exponential decay: \"tau\"."<<endl;
 	      gSystem->Exit(EXIT_FAILURE);
 	    }
 	  } else if (type_ == ExpTrig) {
 	    tau_    = this->findParameter("tau");
 	    deltaM_ = this->findParameter("deltaM");
 	    foundParams &= (tau_ != 0);
 	    foundParams &= (deltaM_ != 0);
 	    if ((this->nParameters() != (3*nGauss_-1+2)) || (!foundParams)) {
 	      cerr<<"ERROR in LauDecayTimePdf::initialise : ExpTrig type PDF requires:"<<endl;
 	      cerr<<"	                                     - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\""<<endl;
 	      cerr<<"	                                     - nGauss-1 fractions: \"frac_i\", where i!=0"<<endl;
 	      cerr<<"                                      - the lifetime of the exponential decay: \"tau\""<<endl;
 	      cerr<<"                                      - the oscillation frequency: \"deltaM\"."<<endl;
 	      gSystem->Exit(EXIT_FAILURE);
 	    }
 	  } else if (type_ == ExpHypTrig) {
 	    tau_        = this->findParameter("tau");
 	    deltaM_     = this->findParameter("deltaM");
 	    deltaGamma_ = this->findParameter("deltaGamma");
 	    foundParams &= (tau_ != 0);
 	    foundParams &= (deltaM_ != 0);
 	    foundParams &= (deltaGamma_ != 0);
 	    if ((this->nParameters() != (3*nGauss_-1+3)) || (!foundParams)) {
 	      cerr<<"ERROR in LauDecayTimePdf::initialise : ExpHypTrig type PDF requires:"<<endl;
 	      cerr<<"	                                     - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\""<<endl;
 	      cerr<<"	                                     - nGauss-1 fractions: \"frac_i\", where i!=0"<<endl;
 	      cerr<<"                                      - the lifetime of the exponential decay: \"tau\""<<endl;
 	      cerr<<"                                      - the oscillation frequency: \"deltaM\"."<<endl;
 	      cerr<<"                                      - the width difference: \"deltaGamma\"."<<endl;
 	      gSystem->Exit(EXIT_FAILURE);
 	    }
 	  } else if (type_ == DeltaExp) {
 	    tau_ = this->findParameter("tau");
 	    fracPrompt_ = this->findParameter("frac_prompt");
 	    foundParams &= (tau_ != 0);
 	    foundParams &= (fracPrompt_ != 0);
 	    if ((this->nParameters() != (3*nGauss_-1+2)) || (!foundParams)) {
 	      cerr<<"ERROR in LauDecayTimePdf::initialise : DeltaExp type PDF requires:"<<endl;
 	      cerr<<"	                                     - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\""<<endl;
 	      cerr<<"	                                     - nGauss-1 fractions: \"frac_i\", where i!=0"<<endl;
 	      cerr<<"                                      - the lifetime of the exponential decay: \"tau\""<<endl;
 	      cerr<<"                                      - the fraction of the prompt part: \"frac_prompt\"."<<endl;
 	      gSystem->Exit(EXIT_FAILURE);
 	    }
 	  } else if (type_ == SimFitNormBd) {
 	    tau_    = this->findParameter("tau");
 	    deltaM_ = this->findParameter("deltaM");
 	    foundParams &= (tau_ != 0);
 	    foundParams &= (deltaM_ != 0);
 	    if ((this->nParameters() != (3*nGauss_-1+2)) || (!foundParams)) {
 	      cerr<<"ERROR in LauDecayTimePdf::initialise : SimFitNormBd type PDF requires:"<<endl;
 	      cerr<<"	                                     - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\""<<endl;
 	      cerr<<"	                                     - nGauss-1 fractions: \"frac_i\", where i!=0"<<endl;
 	      cerr<<"                                      - the lifetime of the exponential decay: \"tau\""<<endl;
 	      cerr<<"                                      - the oscillation frequency: \"deltaM\"."<<endl;
 	      gSystem->Exit(EXIT_FAILURE);
 	    }
 	  } else if (type_ == SimFitSigBd) {
 	    tau_    = this->findParameter("tau");
 	    deltaM_ = this->findParameter("deltaM");
 	    foundParams &= (tau_ != 0);
 	    foundParams &= (deltaM_ != 0);
 	    if ((this->nParameters() != (3*nGauss_-1+2)) || (!foundParams)) {
 	      cerr<<"ERROR in LauDecayTimePdf::initialise : SimFitSigBd type PDF requires:"<<endl;
 	      cerr<<"	                                     - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\""<<endl;
 	      cerr<<"	                                     - nGauss-1 fractions: \"frac_i\", where i!=0"<<endl;
 	      cerr<<"                                      - the lifetime of the exponential decay: \"tau\""<<endl;
 	      cerr<<"                                      - the oscillation frequency: \"deltaM\"."<<endl;
 	      gSystem->Exit(EXIT_FAILURE);
 	    }
 	  } else if (type_ == SimFitNormBs) {
 	    tau_        = this->findParameter("tau");
 	    deltaM_     = this->findParameter("deltaM");
 	    deltaGamma_ = this->findParameter("deltaGamma");
 	    foundParams &= (tau_ != 0);
 	    foundParams &= (deltaM_ != 0);
 	    foundParams &= (deltaGamma_ != 0);
 	    if ((this->nParameters() != (3*nGauss_-1+3)) || (!foundParams)) {
 	      cerr<<"ERROR in LauDecayTimePdf::initialise : SimFitNormBs type PDF requires:"<<endl;
 	      cerr<<"	                                     - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\""<<endl;
 	      cerr<<"	                                     - nGauss-1 fractions: \"frac_i\", where i!=0"<<endl;
 	      cerr<<"                                      - the lifetime of the exponential decay: \"tau\""<<endl;
 	      cerr<<"                                      - the oscillation frequency: \"deltaM\"."<<endl;
 	      cerr<<"                                      - the width difference: \"deltaGamma\"."<<endl;
 	      gSystem->Exit(EXIT_FAILURE);
 	    }
 	  } else if (type_ == SimFitSigBs) {
 	    tau_        = this->findParameter("tau");
 	    deltaM_     = this->findParameter("deltaM");
 	    deltaGamma_ = this->findParameter("deltaGamma");
 	    foundParams &= (tau_ != 0);
 	    foundParams &= (deltaM_ != 0);
 	    foundParams &= (deltaGamma_ != 0);
 	    if ((this->nParameters() != (3*nGauss_-1+3)) || (!foundParams)) {
 	      cerr<<"ERROR in LauDecayTimePdf::initialise : SimFitSigBs type PDF requires:"<<endl;
 	      cerr<<"	                                     - 2 parameters per Gaussian (i): \"mean_i\" and \"sigma_i\""<<endl;
 	      cerr<<"	                                     - nGauss-1 fractions: \"frac_i\", where i!=0"<<endl;
 	      cerr<<"                                      - the lifetime of the exponential decay: \"tau\""<<endl;
 	      cerr<<"                                      - the oscillation frequency: \"deltaM\"."<<endl;
 	      cerr<<"                                      - the width difference: \"deltaGamma\"."<<endl;
 	      gSystem->Exit(EXIT_FAILURE);
 	    }
 	  }
 	}
 	// Cache the normalisation factor
 	//this->calcNorm();
 }
 
 void LauDecayTimePdf::cacheInfo(const LauFitDataTree& inputData)
 {
 	Bool_t hasBranch = inputData.haveBranch(this->varName());
 	if (!hasBranch) {
 		cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Input data does not contain variable \""<<this->varName()<<"\"."<<endl;
 		return;
 	}
 	hasBranch = inputData.haveBranch(this->varErrName());
 	if (!hasBranch) {
 		cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Input data does not contain variable \""<<this->varErrName()<<"\"."<<endl;
 		return;
 	}
 
 	// Pass the data to the decay-time error PDF for caching
 	if ( errHist_ ) {
 		errHist_->cacheInfo(inputData);
 	}
 
 	if (type_ == Hist){
 	  // Pass the data to the decay-time PDF for caching
 	  if ( pdfHist_ ) {
 	    pdfHist_->cacheInfo(inputData);
 	  }
 	}else{
 
 		// determine whether we are caching our PDF value
 		//TODO
 		//Bool_t doCaching( this->nFixedParameters() == this->nParameters() );
 		//this->cachePDF( doCaching );
 
 		// clear the vectors and reserve enough space
 		const UInt_t nEvents = inputData.nEvents();
 		abscissas_.clear(); abscissas_.reserve(nEvents);
 		abscissaErrors_.clear(); abscissaErrors_.reserve(nEvents);
 		expTerms_.clear(); expTerms_.reserve(nEvents);
 		cosTerms_.clear(); cosTerms_.reserve(nEvents);
 		sinTerms_.clear(); sinTerms_.reserve(nEvents);
 		coshTerms_.clear(); coshTerms_.reserve(nEvents);
 		sinhTerms_.clear(); sinhTerms_.reserve(nEvents);
 		normTermsExp_.clear(); normTermsExp_.reserve(nEvents);
 		normTermsCosh_.clear(); normTermsCosh_.reserve(nEvents);
 		normTermsSinh_.clear(); normTermsSinh_.reserve(nEvents);
         	effiTerms_.clear(); effiTerms_.reserve(nEvents);
 
 		for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
 
 			const LauFitData& dataValues = inputData.getData(iEvt);
 
 			LauFitData::const_iterator iter = dataValues.find(this->varName());
 			const Double_t abscissa = iter->second;
 
 			if (abscissa > this->maxAbscissa() ||  abscissa < this->minAbscissa()) {
 				cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Given value of the decay time: "<<abscissa<<
 					" outside allowed range: ["<<this->minAbscissa()<<","<<this->maxAbscissa()<<"]."<<endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 
 			abscissas_.push_back( abscissa );
 
 			iter = dataValues.find(this->varErrName());
 			Double_t abscissaErr = iter->second;
 
 			if (abscissaErr > this->maxAbscissaError() ||  abscissaErr < this->minAbscissaError()) {
 				cerr<<"ERROR in LauDecayTimePdf::cacheInfo : Given value of the decay-time error: "<<abscissaErr<<
 					" outside allowed range: ["<<this->minAbscissaError()<<","<<this->maxAbscissaError()<<"]."<<endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 
 			abscissaErrors_.push_back(abscissaErr);
 
 			this->calcLikelihoodInfo(abscissa, abscissaErr, iEvt);
 			expTerms_.push_back(expTerm_);
 			cosTerms_.push_back(cosTerm_);
 			sinTerms_.push_back(sinTerm_);
 			coshTerms_.push_back(coshTerm_);
 			sinhTerms_.push_back(sinhTerm_);
 			normTermsExp_.push_back(normTermExp_);
 			normTermsCosh_.push_back(normTermCosh_);
 			normTermsSinh_.push_back(normTermSinh_);
 			
 			if(effiFun_) effiTerms_.push_back(effiFun_->evaluate(abscissa));
 			else effiTerms_.push_back(1.0);
 		}
 	}
 }
 
 void LauDecayTimePdf::calcLikelihoodInfo(UInt_t iEvt)
 {
         if (type_ == Hist){
 	  if ( pdfHist_ ) {
 	    pdfHist_->calcLikelihoodInfo(iEvt);
 	    pdfTerm_ = pdfHist_->getLikelihood();
 	  } else {
 	    pdfTerm_ = 1.0;
 	  }
         }else{
 	expTerm_       = expTerms_[iEvt];
 	cosTerm_       = cosTerms_[iEvt];
 	sinTerm_       = sinTerms_[iEvt];
 	coshTerm_      = coshTerms_[iEvt];
 	sinhTerm_      = sinhTerms_[iEvt];
 	normTermExp_   = normTermsExp_[iEvt];
 	normTermCosh_  = normTermsCosh_[iEvt];
 	normTermSinh_  = normTermsSinh_[iEvt];
         }
 	if ( errHist_ ) {
 		errHist_->calcLikelihoodInfo(iEvt);
 		errTerm_ = errHist_->getLikelihood();
 	} else {
 		errTerm_ = 1.0;
 	}
 	
 	const Double_t abscissa = abscissas_[iEvt];
 
 	//Parameters will change in some cases update things!
 	if (type_ == SimFitNormBd || type_ == SimFitSigBd || type_ == SimFitNormBs || type_ == SimFitSigBs){
 		const Double_t abscissaErr = abscissaErrors_[iEvt];
 		this->calcLikelihoodInfo(abscissa,abscissaErr,iEvt);
 	}
 
 	if ( effiFun_ ) {
 		this->updateEffiSpline(effiPars_);
 		effiTerm_ = effiFun_->evaluate(abscissa);
 		if (effiTerm_>1.0){effiTerm_=1.0;}
 		if (effiTerm_<0.0){effiTerm_=0.0;}
 	} else {
 		effiTerm_ = 1.0;
 	}
 
 
 	// TODO need a check in here that none of the floating parameter values have changed
 	// If they have, then we need to recalculate all or some of the terms
 	/*
 	if ( parsChanged ) {
 		const Double_t abscissa = abscissas_[iEvt][0];
 		const Double_t abscissaErr = abscissaErrors_[iEvt];
 		this->calcLikelihoodInfo(abscissa, abscissaErr);
 	}
 	*/
 }
 
 void LauDecayTimePdf::calcLikelihoodInfo(Double_t abscissa, UInt_t iEvt)
 {
 	// Check whether any of the gaussians should be scaled - if any of them should we need the per-event error
 	Bool_t scale(kFALSE);
 	for (std::vector<Bool_t>::const_iterator iter = scaleMeans_.begin(); iter != scaleMeans_.end(); ++iter) {
 		scale |= (*iter);
 	}
 	for (std::vector<Bool_t>::const_iterator iter = scaleWidths_.begin(); iter != scaleWidths_.end(); ++iter) {
 		scale |= (*iter);
 	}
 	if (scale) {
 		cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfo : Per-event error on Delta t not provided, cannot calculate anything."<<endl;
 		return;
 	} else {
 		this->calcLikelihoodInfo(abscissa, 0.0,iEvt);
 	}
 }
 
 void LauDecayTimePdf::calcLikelihoodInfo(Double_t abscissa, Double_t abscissaErr, UInt_t iEvt)
 {
 	if (abscissa > this->maxAbscissa() ||  abscissa < this->minAbscissa()) {
 		cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfo : Given value of the decay time: "<<abscissa<<
 			" outside allowed range: ["<<this->minAbscissa()<<","<<this->maxAbscissa()<<"]."<<endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	if (abscissaErr > this->maxAbscissaError() ||  abscissaErr < this->minAbscissaError()) {
 		cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfo : Given value of Delta t error: "<<abscissaErr<<
 			" outside allowed range: ["<<this->minAbscissaError()<<","<<this->maxAbscissaError()<<"]."<<endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	// Initialise the various terms to zero
 	if (type_ == Hist){
 	  if ( pdfHist_ ) {
 	    pdfHist_->calcLikelihoodInfo(abscissa);
 	    pdfTerm_ = pdfHist_->getLikelihood();
 	  } else {
 	    pdfTerm_ = 1.0;
 	  }
         }else{
 
 	// Reset the state to Good
 	this->state(Good);
 
 	// If we're not using the resolution function calculate the simple terms and return
 	if (!this->doSmearing()) {
 		this->calcNonSmearedTerms(abscissa, iEvt);
 		return;
 	}
 
 	//TODO how much to be added below for SimFitNormBd/SimFitNormBs/SimFitSigBd/SimFitSigBs
 
 	// Get all the up to date parameter values
 	std::vector<Double_t> frac(nGauss_);
 	std::vector<Double_t> mean(nGauss_);
 	std::vector<Double_t> sigma(nGauss_);
 	Double_t tau(0.0);
 	Double_t deltaM(0.0);
 	Double_t fracPrompt(0.0);
 	Double_t Delta_gamma(0.0);
 	frac[0] = 1.0;
 	for (UInt_t i(0); i<nGauss_; ++i) {
 		mean[i] = mean_[i]->value();
 		sigma[i] = sigma_[i]->value();
 		if (i != 0) {
 			frac[i] = frac_[i-1]->value();
 			frac[0] -= frac[i];
 		}
 	}
 	if (type_ != Delta) {
 		tau = tau_->value();
 		if (type_ == ExpTrig) {
 			deltaM = deltaM_->value();
 		}
 		if (type_ == DeltaExp) {
 			fracPrompt = fracPrompt_->value();
 		}
 		if (type_ == ExpHypTrig){
 			deltaM = deltaM_->value();
 			Delta_gamma  = deltaGamma_->value();
 		}
 	}
 
 	// Scale the gaussian parameters by the per-event error on Delta t (if appropriate)
 	for (UInt_t i(0); i<nGauss_; ++i) {
 		if (scaleMeans_[i]) {
 			mean[i] *= abscissaErr;
 		}
 		if (scaleWidths_[i]) {
 			sigma[i] *= abscissaErr;
 		}
 	}
 
 	// Calculate term needed by every type
 	std::vector<Double_t> x(nGauss_);
 	const Double_t xMax = this->maxAbscissa();
 	const Double_t xMin = this->minAbscissa();
 	for (UInt_t i(0); i<nGauss_; ++i) {
 		x[i] = abscissa - mean[i];
 	}
 
 	// TODO, what to do with this
 	Double_t value(0.0);
 
 	if (type_ == Delta || type_ == DeltaExp) {
 		// Calculate the gaussian function(s)
 		for (UInt_t i(0); i<nGauss_; ++i) {
 			if (TMath::Abs(sigma[i]) > 1e-10) {
 				Double_t exponent(0.0);
 				Double_t norm(0.0);
 				Double_t scale = LauConstants::root2*sigma[i];
 				Double_t scale2 = LauConstants::rootPiBy2*sigma[i];
 				exponent = -0.5*x[i]*x[i]/(sigma[i]*sigma[i]);
 				norm = scale2*(TMath::Erf((xMax - mean[i])/scale)
 						- TMath::Erf((xMin - mean[i])/scale));
 				value += frac[i]*TMath::Exp(exponent)/norm;
 			}
 		}
 	}
 
 	if (type_ != Delta) {
 
 		std::vector<Double_t> expTerms(nGauss_);
 		std::vector<Double_t> cosTerms(nGauss_);
 		std::vector<Double_t> sinTerms(nGauss_);
 		std::vector<Double_t> coshTerms(nGauss_);
 		std::vector<Double_t> sinhTerms(nGauss_);
 
 		std::vector<Double_t> expTermsNorm(nGauss_);
 		// TODO - TEL changed this name to make it compile - please check!
 		std::vector<Double_t> SinhTermsNorm(nGauss_);
 
 		// Calculate values of the PDF convoluated with each Gaussian for a given value of the abscsissa
 		for (UInt_t i(0); i<nGauss_; ++i) {
 
                         // Typical case (1): B0/B0bar
                         if (type_ == ExpTrig) {
 
                                 // LHCb convention, i.e. convolution evaluate between 0 and inf
                                 if (method_ == DecayTime) {
 
                                         // Exponential term
                                         Double_t termExponent = (pow(sigma[i], 2) - 2.0 * tau * x[i])/(2.0 * pow(tau, 2));
                                         Double_t termErfc     = (pow(sigma[i], 2) - tau * x[i])/(LauConstants::root2 * tau * sigma[i]);
 			                expTerms[i] = (1.0/2.0) * TMath::Exp(termExponent) * TMath::Erfc(termErfc);
 
 					Double_t exponentTermRe, exponentTermIm;
 					this->calcTrigExponent(deltaM, tau, x[i], sigma[i], exponentTermRe, exponentTermIm);
 
                                         // Elements related to the trigonometric function, i.e. convolution of Exp*Sin or Cos with Gauss
                                         Double_t sinTrigTermRe, sinTrigTermIm, cosTrigTermRe, cosTrigTermIm;
 					this->calcTrigConv(deltaM, tau, x[i], sigma[i], sinTrigTermRe, sinTrigTermIm, kFALSE);
 					this->calcTrigConv(deltaM, tau, x[i], sigma[i], cosTrigTermRe, cosTrigTermIm, kTRUE);
 
                                         // Combining elements of the full pdf
                                         LauComplex zExp(exponentTermRe, exponentTermIm);
                                         LauComplex zTrigSin(sinTrigTermRe, sinTrigTermIm);
                                         LauComplex zTrigCos(cosTrigTermRe, cosTrigTermIm);
 
                                         LauComplex sinConv = zExp * zTrigSin;
                                         LauComplex cosConv = zExp * zTrigCos;
                                         sinConv.scale(1.0/4.0);
                                         cosConv.scale(1.0/4.0);
 
                                         // Cosine*Exp and Sine*Exp terms
 					cosTerms[i] = cosConv.re();
 					sinTerms[i] = sinConv.im();
 
 			                // Normalisation
 			                Double_t umax = xMax - mean[i];
 			                Double_t umin = xMin - mean[i];
 
                                         expTermsNorm[i] = (1.0/2.0) * tau * (-1.0 + TMath::Erf(umax/(LauConstants::root2 * sigma[i])) + TMath::Erfc(umin/(LauConstants::root2 * sigma[i])) +
                                                       TMath::Exp((pow(sigma[i], 2) - 2.0 * tau * (xMax + xMin - mean[i]))/(2.0 * pow(tau, 2))) *
                                                       (TMath::Exp(xMax/tau) * TMath::Erfc((pow(sigma[i], 2) - xMin)/(LauConstants::root2 * tau))) +
                                                       (TMath::Exp(xMin/tau) * TMath::Erfc((pow(sigma[i], 2) - xMax)/(LauConstants::root2 * tau))));
                                 } else {
 
 
                                 }
                         }
                         // Typical case (2): B0s/B0sbar
                         if (type_ == ExpHypTrig) {
                           // LHCb convention
                                 if (method_ == DecayTime) {
                                   // Convolution of Exp*cosh (Exp*sinh) with a gaussian
                                   //Double_t OverallExpFactor = 0.25*TMath::Exp(-(x[i]-mean[i])*(x[i]-mean[i])/(2*sigma[i]*sigma[i]));
                                   //Double_t ExpFirstTerm = TMath::Exp((2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))*(2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))/(8*sigma[i]*sigma[i]*tau*tau));
                                   //Double_t ExpSecondTerm = TMath::Exp((2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))*(2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))/(8*sigma[i]*sigma[i]*tau*tau));
                                   //Double_t ErfFirstTerm = TMath::Erf((2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
                                   //Double_t ErfSecondTerm = TMath::Erf((2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
 
                                   //Double_t sinhConv = OverallExpFactor*(ExpFirstTerm*(1+ErfFirstTerm) + ExpSecondTerm*(-1+ErfSecondTerm));
                                   //Double_t coshConv = OverallExpFactor*(ExpFirstTerm*(1+ErfFirstTerm) - ExpSecondTerm*(-1+ErfSecondTerm));
 
                                   //cosTerms[i] = sinhConv;
       			          //            sinTerms[i] = coshConv;
 
                                   //TODO: check this formula and try to simplify it!
                                   double OverallExpTerm_max = (1/(2*(-4 + Delta_gamma*Delta_gamma*tau*tau)))*tau*TMath::Exp(-0.5*Delta_gamma*(xMax + mean[i]) - xMax/tau);
                                   double ErfTerm_max = -2*Delta_gamma*tau*TMath::Exp(0.5*Delta_gamma*(xMax+mean[i])+xMax/tau)*TMath::Erf((xMax-mean[i])/(TMath::Sqrt(2)*sigma[i]));
                                   double ExpFirstTerm_max = TMath::Exp(xMax*Delta_gamma+(sigma[i]*sigma[i]*(-2 + Delta_gamma*tau)*(-2 + Delta_gamma*tau))/(8*tau*tau));
                                   double ErfcFirstTerm_max = TMath::Erfc((2*(-xMax + mean[i])*tau + sigma[i]*sigma[i]*(2 - Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
                                   double ExpSecondTerm_max = TMath::Exp(Delta_gamma*mean[i] + (sigma[i]*sigma[i]*(2 + Delta_gamma*tau)*(2 + Delta_gamma*tau))/(8*tau*tau));
                                   double ErfcSecondTerm_max = TMath::Erfc((2*(-xMax + mean[i])*tau + sigma[i]*sigma[i]*(2 + Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
                                   double MaxVal= OverallExpTerm_max*(ErfTerm_max + TMath::Exp(mean[i]/tau)*(ExpFirstTerm_max*(2+Delta_gamma*tau)* ErfcFirstTerm_max + ExpSecondTerm_max*(-2+Delta_gamma*tau)* ErfcSecondTerm_max));
 
                                   double OverallExpTerm_min = (1/(2*(-4 + Delta_gamma*Delta_gamma*tau*tau)))*tau*TMath::Exp(-0.5*Delta_gamma*(xMin + mean[i]) - xMin/tau);
                                   double ErfTerm_min = -2*Delta_gamma*tau*TMath::Exp(0.5*Delta_gamma*(xMin+mean[i])+xMin/tau)*TMath::Erf((xMin-mean[i])/(TMath::Sqrt(2)*sigma[i]));
                                   double ExpFirstTerm_min = TMath::Exp(xMin*Delta_gamma+(sigma[i]*sigma[i]*(-2 + Delta_gamma*tau)*(-2 + Delta_gamma*tau))/(8*tau*tau));
                                   double ErfcFirstTerm_min = TMath::Erfc((2*(-xMin + mean[i])*tau + sigma[i]*sigma[i]*(2 - Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
 				  // TODO - TEL added this (currently identical to ExpSecondTerm_max) to get this to compile - please check!!
                                   double ExpSecondTerm_min = TMath::Exp(Delta_gamma*mean[i] + (sigma[i]*sigma[i]*(2 + Delta_gamma*tau)*(2 + Delta_gamma*tau))/(8*tau*tau));
                                   double ErfcSecondTerm_min = TMath::Erfc((2*(-xMin + mean[i])*tau + sigma[i]*sigma[i]*(2 + Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
                                   double minVal= OverallExpTerm_min*(ErfTerm_min + TMath::Exp(mean[i]/tau)*(ExpFirstTerm_min*(2+Delta_gamma*tau)* ErfcFirstTerm_min + ExpSecondTerm_min*(-2+Delta_gamma*tau)* ErfcSecondTerm_min));
                                   SinhTermsNorm[i] = MaxVal - minVal;
 
                                 } else {
 
                                 }
                         }
 
 		}
 
 		for (UInt_t i(0); i<nGauss_; ++i) {
 			expTerm_     += frac[i]*expTerms[i];
 			cosTerm_     += frac[i]*cosTerms[i];
 			sinTerm_     += frac[i]*sinTerms[i];
       			coshTerm_    += frac[i]*coshTerms[i];
 			sinhTerm_    += frac[i]*sinhTerms[i];
 
 
 
                         normTermExp_ += frac[i]*expTermsNorm[i];
                         //normTermSinh_ += frac[i]*SinhTermsNorm[i];
 		}
 
 		if (type_ == DeltaExp) {
 			value *= fracPrompt;
 			value += (1.0-fracPrompt)*expTerm_;
 		} else {
 			value = expTerm_;
 		}
 	}
 	}
 	if ( errHist_ ) {
 		errHist_->calcLikelihoodInfo(abscissaErr);
 		errTerm_ = errHist_->getLikelihood();
 	} else {
 		errTerm_ = 1.0;
 	}
 
 	if ( effiFun_ ) {
 		this->updateEffiSpline(effiPars_);
 		effiTerm_ = effiFun_->evaluate(abscissa);
 		if (effiTerm_>1.0){effiTerm_=1.0;}
 		if (effiTerm_<0.0){effiTerm_=0.0;}
 	} else {
 		effiTerm_ = 1.0;
 	}
 }
 
 void LauDecayTimePdf::calcLikelihoodInfoGen(Double_t abscissa, Double_t abscissaErr, Int_t curEvtTagFlv, Double_t curEvtMistag, Int_t curEvtTagCat, Int_t curEvtTrueTagFlv)
 {
 	if (abscissa > this->maxAbscissa() ||  abscissa < this->minAbscissa()) {
 		cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfoGen : Given value of the decay time: "<<abscissa<<
 			" outside allowed range: ["<<this->minAbscissa()<<","<<this->maxAbscissa()<<"]."<<endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	if (abscissaErr > this->maxAbscissaError() ||  abscissaErr < this->minAbscissaError()) {
 		cerr<<"ERROR in LauDecayTimePdf::calcLikelihoodInfoGen : Given value of Delta t error: "<<abscissaErr<<
 			" outside allowed range: ["<<this->minAbscissaError()<<","<<this->maxAbscissaError()<<"]."<<endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	// Initialise the various terms to zero
 	if (type_ == Hist){
 	  if ( pdfHist_ ) {
 	    pdfHist_->calcLikelihoodInfo(abscissa);
 	    pdfTerm_ = pdfHist_->getLikelihood();
 	  } else {
 	    pdfTerm_ = 1.0;
 	  }
         }else{
 
 	// Reset the state to Good
 	this->state(Good);
 
 	// If we're not using the resolution function calculate the simple terms and return
 	if (!this->doSmearing()) {
 		this->calcNonSmearedTermsGen(abscissa, curEvtTagFlv, curEvtMistag, curEvtTagCat, curEvtTrueTagFlv);
 		return;
 	}
 
 	//TODO how much to be added below for SimFitNormBd/SimFitNormBs/SimFitSigBd/SimFitSigBs
 
 	// Get all the up to date parameter values
 	std::vector<Double_t> frac(nGauss_);
 	std::vector<Double_t> mean(nGauss_);
 	std::vector<Double_t> sigma(nGauss_);
 	Double_t tau(0.0);
 	Double_t deltaM(0.0);
 	Double_t fracPrompt(0.0);
 	Double_t Delta_gamma(0.0);
 	frac[0] = 1.0;
 	for (UInt_t i(0); i<nGauss_; ++i) {
 		mean[i] = mean_[i]->value();
 		sigma[i] = sigma_[i]->value();
 		if (i != 0) {
 			frac[i] = frac_[i-1]->value();
 			frac[0] -= frac[i];
 		}
 	}
 	if (type_ != Delta) {
 		tau = tau_->value();
 		if (type_ == ExpTrig) {
 			deltaM = deltaM_->value();
 		}
 		if (type_ == DeltaExp) {
 			fracPrompt = fracPrompt_->value();
 		}
 		if (type_ == ExpHypTrig){
 			deltaM = deltaM_->value();
 			Delta_gamma  = deltaGamma_->value();
 		}
 	}
 
 	// Scale the gaussian parameters by the per-event error on Delta t (if appropriate)
 	for (UInt_t i(0); i<nGauss_; ++i) {
 		if (scaleMeans_[i]) {
 			mean[i] *= abscissaErr;
 		}
 		if (scaleWidths_[i]) {
 			sigma[i] *= abscissaErr;
 		}
 	}
 
 	// Calculate term needed by every type
 	std::vector<Double_t> x(nGauss_);
 	const Double_t xMax = this->maxAbscissa();
 	const Double_t xMin = this->minAbscissa();
 	for (UInt_t i(0); i<nGauss_; ++i) {
 		x[i] = abscissa - mean[i];
 	}
 
 	// TODO, what to do with this
 	Double_t value(0.0);
 
 	if (type_ == Delta || type_ == DeltaExp) {
 		// Calculate the gaussian function(s)
 		for (UInt_t i(0); i<nGauss_; ++i) {
 			if (TMath::Abs(sigma[i]) > 1e-10) {
 				Double_t exponent(0.0);
 				Double_t norm(0.0);
 				Double_t scale = LauConstants::root2*sigma[i];
 				Double_t scale2 = LauConstants::rootPiBy2*sigma[i];
 				exponent = -0.5*x[i]*x[i]/(sigma[i]*sigma[i]);
 				norm = scale2*(TMath::Erf((xMax - mean[i])/scale)
 						- TMath::Erf((xMin - mean[i])/scale));
 				value += frac[i]*TMath::Exp(exponent)/norm;
 			}
 		}
 	}
 
 	if (type_ != Delta) {
 
 		std::vector<Double_t> expTerms(nGauss_);
 		std::vector<Double_t> cosTerms(nGauss_);
 		std::vector<Double_t> sinTerms(nGauss_);
 		std::vector<Double_t> coshTerms(nGauss_);
 		std::vector<Double_t> sinhTerms(nGauss_);
 
 		std::vector<Double_t> expTermsNorm(nGauss_);
 		// TODO - TEL changed this name to make it compile - please check!
 		std::vector<Double_t> SinhTermsNorm(nGauss_);
 
 		// Calculate values of the PDF convoluated with each Gaussian for a given value of the abscsissa
 		for (UInt_t i(0); i<nGauss_; ++i) {
 
                         // Typical case (1): B0/B0bar
                         if (type_ == ExpTrig) {
 
                                 // LHCb convention, i.e. convolution evaluate between 0 and inf
                                 if (method_ == DecayTime) {
 
                                         // Exponential term
                                         Double_t termExponent = (pow(sigma[i], 2) - 2.0 * tau * x[i])/(2.0 * pow(tau, 2));
                                         Double_t termErfc     = (pow(sigma[i], 2) - tau * x[i])/(LauConstants::root2 * tau * sigma[i]);
 			                expTerms[i] = (1.0/2.0) * TMath::Exp(termExponent) * TMath::Erfc(termErfc);
 
 					Double_t exponentTermRe, exponentTermIm;
 					this->calcTrigExponent(deltaM, tau, x[i], sigma[i], exponentTermRe, exponentTermIm);
 
                                         // Elements related to the trigonometric function, i.e. convolution of Exp*Sin or Cos with Gauss
                                         Double_t sinTrigTermRe, sinTrigTermIm, cosTrigTermRe, cosTrigTermIm;
 					this->calcTrigConv(deltaM, tau, x[i], sigma[i], sinTrigTermRe, sinTrigTermIm, kFALSE);
 					this->calcTrigConv(deltaM, tau, x[i], sigma[i], cosTrigTermRe, cosTrigTermIm, kTRUE);
 
                                         // Combining elements of the full pdf
                                         LauComplex zExp(exponentTermRe, exponentTermIm);
                                         LauComplex zTrigSin(sinTrigTermRe, sinTrigTermIm);
                                         LauComplex zTrigCos(cosTrigTermRe, cosTrigTermIm);
 
                                         LauComplex sinConv = zExp * zTrigSin;
                                         LauComplex cosConv = zExp * zTrigCos;
                                         sinConv.scale(1.0/4.0);
                                         cosConv.scale(1.0/4.0);
 
                                         // Cosine*Exp and Sine*Exp terms
 					cosTerms[i] = cosConv.re();
 					sinTerms[i] = sinConv.im();
 
 			                // Normalisation
 			                Double_t umax = xMax - mean[i];
 			                Double_t umin = xMin - mean[i];
 
                                         expTermsNorm[i] = (1.0/2.0) * tau * (-1.0 + TMath::Erf(umax/(LauConstants::root2 * sigma[i])) + TMath::Erfc(umin/(LauConstants::root2 * sigma[i])) +
                                                       TMath::Exp((pow(sigma[i], 2) - 2.0 * tau * (xMax + xMin - mean[i]))/(2.0 * pow(tau, 2))) *
                                                       (TMath::Exp(xMax/tau) * TMath::Erfc((pow(sigma[i], 2) - xMin)/(LauConstants::root2 * tau))) +
                                                       (TMath::Exp(xMin/tau) * TMath::Erfc((pow(sigma[i], 2) - xMax)/(LauConstants::root2 * tau))));
                                 } else {
 
 
                                 }
                         }
                         // Typical case (2): B0s/B0sbar
                         if (type_ == ExpHypTrig) {
                           // LHCb convention
                                 if (method_ == DecayTime) {
                                   // Convolution of Exp*cosh (Exp*sinh) with a gaussian
                                   //Double_t OverallExpFactor = 0.25*TMath::Exp(-(x[i]-mean[i])*(x[i]-mean[i])/(2*sigma[i]*sigma[i]));
                                   //Double_t ExpFirstTerm = TMath::Exp((2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))*(2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))/(8*sigma[i]*sigma[i]*tau*tau));
                                   //Double_t ExpSecondTerm = TMath::Exp((2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))*(2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))/(8*sigma[i]*sigma[i]*tau*tau));
                                   //Double_t ErfFirstTerm = TMath::Erf((2*(x[i]-mean[i])*tau+sigma[i]*sigma[i]*(-2+Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
                                   //Double_t ErfSecondTerm = TMath::Erf((2*(-x[i]+mean[i])*tau+sigma[i]*sigma[i]*(2+Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
 
                                   //Double_t sinhConv = OverallExpFactor*(ExpFirstTerm*(1+ErfFirstTerm) + ExpSecondTerm*(-1+ErfSecondTerm));
                                   //Double_t coshConv = OverallExpFactor*(ExpFirstTerm*(1+ErfFirstTerm) - ExpSecondTerm*(-1+ErfSecondTerm));
 
                                   //cosTerms[i] = sinhConv;
       			          //            sinTerms[i] = coshConv;
 
                                   //TODO: check this formula and try to simplify it!
                                   double OverallExpTerm_max = (1/(2*(-4 + Delta_gamma*Delta_gamma*tau*tau)))*tau*TMath::Exp(-0.5*Delta_gamma*(xMax + mean[i]) - xMax/tau);
                                   double ErfTerm_max = -2*Delta_gamma*tau*TMath::Exp(0.5*Delta_gamma*(xMax+mean[i])+xMax/tau)*TMath::Erf((xMax-mean[i])/(TMath::Sqrt(2)*sigma[i]));
                                   double ExpFirstTerm_max = TMath::Exp(xMax*Delta_gamma+(sigma[i]*sigma[i]*(-2 + Delta_gamma*tau)*(-2 + Delta_gamma*tau))/(8*tau*tau));
                                   double ErfcFirstTerm_max = TMath::Erfc((2*(-xMax + mean[i])*tau + sigma[i]*sigma[i]*(2 - Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
                                   double ExpSecondTerm_max = TMath::Exp(Delta_gamma*mean[i] + (sigma[i]*sigma[i]*(2 + Delta_gamma*tau)*(2 + Delta_gamma*tau))/(8*tau*tau));
                                   double ErfcSecondTerm_max = TMath::Erfc((2*(-xMax + mean[i])*tau + sigma[i]*sigma[i]*(2 + Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
                                   double MaxVal= OverallExpTerm_max*(ErfTerm_max + TMath::Exp(mean[i]/tau)*(ExpFirstTerm_max*(2+Delta_gamma*tau)* ErfcFirstTerm_max + ExpSecondTerm_max*(-2+Delta_gamma*tau)* ErfcSecondTerm_max));
 
                                   double OverallExpTerm_min = (1/(2*(-4 + Delta_gamma*Delta_gamma*tau*tau)))*tau*TMath::Exp(-0.5*Delta_gamma*(xMin + mean[i]) - xMin/tau);
                                   double ErfTerm_min = -2*Delta_gamma*tau*TMath::Exp(0.5*Delta_gamma*(xMin+mean[i])+xMin/tau)*TMath::Erf((xMin-mean[i])/(TMath::Sqrt(2)*sigma[i]));
                                   double ExpFirstTerm_min = TMath::Exp(xMin*Delta_gamma+(sigma[i]*sigma[i]*(-2 + Delta_gamma*tau)*(-2 + Delta_gamma*tau))/(8*tau*tau));
                                   double ErfcFirstTerm_min = TMath::Erfc((2*(-xMin + mean[i])*tau + sigma[i]*sigma[i]*(2 - Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
 				  // TODO - TEL added this (currently identical to ExpSecondTerm_max) to get this to compile - please check!!
                                   double ExpSecondTerm_min = TMath::Exp(Delta_gamma*mean[i] + (sigma[i]*sigma[i]*(2 + Delta_gamma*tau)*(2 + Delta_gamma*tau))/(8*tau*tau));
                                   double ErfcSecondTerm_min = TMath::Erfc((2*(-xMin + mean[i])*tau + sigma[i]*sigma[i]*(2 + Delta_gamma*tau))/(2*TMath::Sqrt(2)*sigma[i]*tau));
                                   double minVal= OverallExpTerm_min*(ErfTerm_min + TMath::Exp(mean[i]/tau)*(ExpFirstTerm_min*(2+Delta_gamma*tau)* ErfcFirstTerm_min + ExpSecondTerm_min*(-2+Delta_gamma*tau)* ErfcSecondTerm_min));
                                   SinhTermsNorm[i] = MaxVal - minVal;
 
                                 } else {
 
                                 }
                         }
 
 		}
 
 		for (UInt_t i(0); i<nGauss_; ++i) {
 			expTerm_     += frac[i]*expTerms[i];
 			cosTerm_     += frac[i]*cosTerms[i];
 			sinTerm_     += frac[i]*sinTerms[i];
       			coshTerm_    += frac[i]*coshTerms[i];
 			sinhTerm_    += frac[i]*sinhTerms[i];
 
 
 
                         normTermExp_ += frac[i]*expTermsNorm[i];
                         //normTermSinh_ += frac[i]*SinhTermsNorm[i];
 		}
 
 		if (type_ == DeltaExp) {
 			value *= fracPrompt;
 			value += (1.0-fracPrompt)*expTerm_;
 		} else {
 			value = expTerm_;
 		}
 	}
 	}
 	if ( errHist_ ) {
 		errHist_->calcLikelihoodInfo(abscissaErr);
 		errTerm_ = errHist_->getLikelihood();
 	} else {
 		errTerm_ = 1.0;
 	}
 
 	if ( effiFun_ ) {
 		this->updateEffiSpline(effiPars_);
 		effiTerm_ = effiFun_->evaluate(abscissa);
 		if (effiTerm_>1.0){effiTerm_=1.0;}
 		if (effiTerm_<0.0){effiTerm_=0.0;}
 	} else {
 		effiTerm_ = 1.0;
 	}
 }
 
 void LauDecayTimePdf::calcTrigExponent(Double_t deltaM, Double_t tau, Double_t x, Double_t sigma, Double_t& reTerm, Double_t& imTerm)
 {
 
         Double_t exponentTerm = TMath::Exp(-(2.0 * tau * x + pow(sigma, 2) * (pow(deltaM, 2) * pow(tau, 2) - 1.0))/(2.0 * pow(tau,2)));
         reTerm =   exponentTerm * TMath::Cos(deltaM * (x - pow(sigma,2)/tau));
         imTerm = - exponentTerm * TMath::Sin(deltaM * (x - pow(sigma,2)/tau));
 
 }
 
 void LauDecayTimePdf::calcTrigConv(Double_t deltaM, Double_t tau, Double_t x, Double_t sigma, Double_t& reOutTerm, Double_t& imOutTerm, Bool_t trig)
 {
 
         Double_t reExpTerm, imExpTerm;
         LauComplex zExp;
         LauComplex zTrig1;
         LauComplex zTrig2;
 
         // Calculation for the sine or cosine term
         if (!trig) {
                reExpTerm = TMath::Sin(2.0 * deltaM * (x + pow(sigma,2)/tau));
                imExpTerm = 2.0 * TMath::Sin(pow(deltaM * (x + pow(sigma,2)/tau), 2));
         } else {
                reExpTerm = TMath::Cos(2.0 * deltaM * (x + pow(sigma,2)/tau));
                imExpTerm = TMath::Sin(2.0 * deltaM * (x + pow(sigma,2)/tau));
         }
 
         // Exponential term in front of Erfc/Erfi terms
         zExp.setRealPart(reExpTerm);
         zExp.setImagPart(imExpTerm);
 
         // Nominal Erfc term (common to both sine and cosine expressions
         zTrig1.setRealPart(-(tau * x - pow(sigma,2))/(LauConstants::root2 * tau * sigma));
         zTrig1.setImagPart(-(deltaM * sigma)/ LauConstants::root2);
 
         // Second term for sine (Erfi) or cosine (Erfc) - notice the re-im swap and sign change
         zTrig2.setRealPart(-zTrig1.im());
         zTrig2.setImagPart(-zTrig1.re());
 
         // Calculation of Erfc and Erfi (if necessary)
         LauComplex term1 = ComplexErfc(zTrig1.re(), zTrig1.im());
         LauComplex term2;
         if (!trig) {
                term2 = Erfi(zTrig2.re(), zTrig2.im());
         } else {
                term2 = ComplexErfc(zTrig2.re(), zTrig2.im());
         }
 
         // Multiplying all elemnets of the convolution
         LauComplex output = zExp * term1 + term2;
         reOutTerm = output.re();
         imOutTerm = output.im();
 
 }
 
 LauComplex LauDecayTimePdf::ComplexErf(Double_t x, Double_t y)
 {
   //   Evaluate Erf(x + iy) using an infinite series approximation
   //   From Abramowitz & Stegun (http://people.math.sfu.ca/~cbm/aands/page_299.htm)
   if (x==0){
       // cout << "WARNING: Set x value to 1e-100 to avoid division by 0." << endl;
       x = 1e-100;
   }
   int n = 20; // this cotrols the number of iterations of the sum
   LauComplex ErfTerm(TMath::Erf(x),0.);
   LauComplex CosSineTerm(1-cos(2*x*y), sin(2*x*y));
   CosSineTerm.rescale(TMath::Exp(-x*x)/(2*TMath::Pi()*x));
   LauComplex firstPart = ErfTerm + CosSineTerm;
   LauComplex SumTerm(0,0);
 
   for (int k = 1; k<=n; k++){
       Double_t f_k = 2*x*(1 - cos(2*x*y)*cosh(k*y)) + k*sin(2*x*y)*sinh(k*y);
       Double_t g_k = 2*x*sin(2*x*y)*cosh(k*y) + k*cos(2*x*y)*sinh(k*y);
       LauComplex fgTerm(f_k, g_k);
       fgTerm.rescale(TMath::Exp(-0.25*k*k)/(k*k + 4*x*x));
       SumTerm += fgTerm;
   }
   SumTerm.rescale((2/TMath::Pi())*TMath::Exp(-x*x));
   LauComplex result = firstPart + SumTerm;
   return result;
 }
 
 LauComplex LauDecayTimePdf::Erfi(Double_t x, Double_t y)
 {
   //   Erfi(z) = -I*Erf(I*z) where z = x + iy
   double x_prime = -y;
   double y_prime =  x;
   LauComplex a  = ComplexErf(x_prime, y_prime);
   LauComplex result(a.im(), -a.re());
   return result;
 }
 
 LauComplex LauDecayTimePdf::ComplexErfc(Double_t x, Double_t y)
 {
   //    Erfc(z) = 1 - Erf(z) (z = x + iy)
   LauComplex one(1., 0.);
   LauComplex result = one - ComplexErf(x,y);
   return result;
 }
 
 void LauDecayTimePdf::calcNonSmearedTerms(Double_t abscissa, UInt_t iEvt)
 {
        if (type_ == Hist ){
 	 cerr << "It is a histogrammed PDF" << endl;
          return;
        }
 	if (type_ == Delta) {
 		return;
 	}
 
 	Double_t tau = tau_->value();
 	Double_t deltaM = deltaM_->value();
 
         // Calculate the terms related to cosine and sine not normalised
 	if (type_ == ExpTrig) {
 
                 if (method_ == DecayTime) {
 	                expTerm_  = TMath::Exp(-abscissa/tau)/(2.0*tau);
                 }
                 if (method_ == DecayTimeDiff) {
 	                expTerm_  = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau);
                 }
 		cosTerm_  = TMath::Cos(deltaM*abscissa)*expTerm_;
 		sinTerm_  = TMath::Sin(deltaM*abscissa)*expTerm_;
 		coshTerm_ = expTerm_;
 		sinhTerm_ = 0.0;
 
 	}
 
         // Calculate the terms related to cosine not normalised
 	if (type_ == SimFitNormBd || type_ == SimFitNormBs) {
 
                 if (method_ == DecayTime) {
 	                expTerm_  = TMath::Exp(-abscissa/tau)/(2.0*tau);
                 }
                 if (method_ == DecayTimeDiff) {
 	                expTerm_  = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau);
                 }
 		asymTerm_ = ARaw_->unblindValue();
 
 		Int_t curEvtTagFlv = flavTag_->getEvtTagFlvVals(iEvt);
 		Int_t curEvtTrueTagFlv = flavTag_->getEvtTrueTagVals(iEvt);
 
 		Double_t omega(0);
 		Double_t omegabar(0);
 		if (curEvtTagFlv==1){
 			omegabar = flavTag_->getOmega(iEvt);
 		} else if (curEvtTagFlv==-1){
 			omega = flavTag_->getOmega(iEvt);
 		}
 		
 		Double_t Ap = AProd_->unblindValue(); 
 		omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega);
 		omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega);
 
 		cosTerm_  = TMath::Cos(deltaM*abscissa)*expTerm_*asymTerm_*omegaTerm_*curEvtTrueTagFlv;
 		sinTerm_  = 0.0;
 		coshTerm_ = expTerm_*asymTerm_*omegaTermh_;
 		sinhTerm_ = 0.0;
 
 		if (type_ == SimFitNormBs){
 			Double_t deltaGamma = deltaGamma_->value();
 			coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0);
 		}
 
 	}
         
 	// Calculate the terms related to cosine and sine not normalised
 	if (type_ == SimFitSigBd || type_ == SimFitSigBs) {
 
                 if (method_ == DecayTime) {
 	                expTerm_  = TMath::Exp(-abscissa/tau)/(2.0*tau);
                 }
                 if (method_ == DecayTimeDiff) {
 	                expTerm_  = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau);
                 }
 		Int_t curEvtTagFlv = flavTag_->getEvtTagFlvVals(iEvt);
 
 		Double_t omega(0);
 		Double_t omegabar(0);
 		if (curEvtTagFlv==1){
 			omegabar = flavTag_->getOmega(iEvt);
 		} else if (curEvtTagFlv==-1){
 			omega = flavTag_->getOmega(iEvt);
 		}
 		
 		Double_t Ap = AProd_->unblindValue(); 
 		omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega);
 		omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega);
 
 		cosTerm_  = TMath::Cos(deltaM*abscissa)*expTerm_*omegaTerm_;
 		sinTerm_  = TMath::Sin(deltaM*abscissa)*expTerm_*omegaTerm_;
 		coshTerm_ = expTerm_*omegaTermh_;
 		sinhTerm_ = 0.0;
 		
 		if (type_ == SimFitNormBs){
 			Double_t deltaGamma = deltaGamma_->value();
 			coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0);
 			sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_*omegaTermh_;
 		}
 
 	}
         
         // Calculate the terms related to cosine, sine, cosh and sinh not normalised (no decayTimeDiff implemented)
 	if (type_ == ExpHypTrig) {
 
 		Double_t deltaGamma = deltaGamma_->value();
 	        expTerm_  = TMath::Exp(-abscissa/tau)/(2.0*tau);
 		cosTerm_  = TMath::Cos(deltaM*abscissa)*expTerm_;
 		sinTerm_  = TMath::Sin(deltaM*abscissa)*expTerm_;
 		coshTerm_ = TMath::CosH(deltaGamma*abscissa/2.0)*expTerm_;
 		sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_;
 
 	}
 
 	if ( effiFun_ ) {
 		this->updateEffiSpline(effiPars_);
 		effiTerm_ = effiFun_->evaluate(abscissa);
 		if (effiTerm_>1.0){effiTerm_=1.0;}
 		if (effiTerm_<0.0){effiTerm_=0.0;}
 	} else {
 		effiTerm_ = 1.0;
 	}
 }
 
 void LauDecayTimePdf::calcNonSmearedTermsGen(Double_t abscissa, Int_t curEvtTagFlv, Double_t curEvtMistag, Int_t curEvtTagCat, Int_t curEvtTrueTagFlv)
 {
 	if (type_ == Hist ){
 		cerr << "It is a histogrammed PDF" << endl;
 		return;
 	}
 	if (type_ == Delta) {
 		return;
 	}
 
 	Double_t tau = tau_->value();
 
         // Calculate the terms related to cosine and sine not normalised
 	if (type_ == ExpTrig) {
 
 		Double_t deltaM     = deltaM_->value();
                 if (method_ == DecayTime) {
 	                expTerm_  = TMath::Exp(-abscissa/tau)/(2.0*tau);
                 }
                 if (method_ == DecayTimeDiff) {
 	                expTerm_  = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau);
                 }
 		cosTerm_  = TMath::Cos(deltaM*abscissa)*expTerm_;
 		sinTerm_  = TMath::Sin(deltaM*abscissa)*expTerm_;
 		coshTerm_ = expTerm_;
 		sinhTerm_ = 0.0;
 
 	}
 
         // Calculate the terms related to cosine not normalised
 	if (type_ == SimFitNormBd || type_ == SimFitNormBs) {
 
 		Double_t deltaM     = deltaM_->value();
                 if (method_ == DecayTime) {
 	                expTerm_  = TMath::Exp(-abscissa/tau)/(2.0*tau);
                 }
                 if (method_ == DecayTimeDiff) {
 	                expTerm_  = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau);
                 }
 		asymTerm_ = ARaw_->unblindValue();
 
 		Double_t omega(0);
 		Double_t omegabar(0);
 		if (curEvtTagFlv==1){
 			omegabar = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat);
 		} else if (curEvtTagFlv==-1){
 			omega = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat);
 		}
 		
 		Double_t Ap = AProd_->unblindValue(); 
 		omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega);
 		omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega);
 
 		cosTerm_  = TMath::Cos(deltaM*abscissa)*expTerm_*asymTerm_*omegaTerm_*curEvtTrueTagFlv;
 		sinTerm_  = 0.0;
 		coshTerm_ = expTerm_*asymTerm_*omegaTermh_;
 		sinhTerm_ = 0.0;
 
 		if (type_ == SimFitNormBs){
 			Double_t deltaGamma = deltaGamma_->value();
 			coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0);
 		}
 
 	}
         
 	// Calculate the terms related to cosine and sine not normalised
 	if (type_ == SimFitSigBd || type_ == SimFitSigBs) {
 
 		Double_t deltaM     = deltaM_->value();
                 if (method_ == DecayTime) {
 	                expTerm_  = TMath::Exp(-abscissa/tau)/(2.0*tau);
                 }
                 if (method_ == DecayTimeDiff) {
 	                expTerm_  = TMath::Exp(-TMath::Abs(abscissa)/tau)/(2.0*tau);
                 }
 		Double_t omega(0);
 		Double_t omegabar(0);
 		if (curEvtTagFlv==1){
 			omegabar = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat);
 		} else if (curEvtTagFlv==-1){
 			omega = flavTag_->getOmegaGen(curEvtMistag,curEvtTagFlv,curEvtTagCat);
 		}
 	
 		Double_t Ap = AProd_->unblindValue(); 
 		omegaTerm_ = ((1-Ap)*omegabar-(1+Ap)*omega);
 		omegaTermh_ = ((1-Ap)*omegabar+(1+Ap)*omega);
 		
 		cosTerm_  = TMath::Cos(deltaM*abscissa)*expTerm_*omegaTerm_;
 		sinTerm_  = TMath::Sin(deltaM*abscissa)*expTerm_*omegaTerm_;
 		coshTerm_ = expTerm_*omegaTermh_;
 		sinhTerm_ = 0.0;
 		
 		if (type_ == SimFitNormBs){
 			Double_t deltaGamma = deltaGamma_->value();
 			coshTerm_ *= TMath::CosH(deltaGamma*abscissa/2.0);
 			sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_*omegaTermh_;
 		}
 	}
         
         // Calculate the terms related to cosine, sine, cosh and sinh not normalised (no decayTimeDiff implemented)
 	if (type_ == ExpHypTrig) {
 
 		Double_t deltaM     = deltaM_->value();
 		Double_t deltaGamma = deltaGamma_->value();
 	        expTerm_  = TMath::Exp(-abscissa/tau)/(2.0*tau);
 		cosTerm_  = TMath::Cos(deltaM*abscissa)*expTerm_;
 		sinTerm_  = TMath::Sin(deltaM*abscissa)*expTerm_;
 		coshTerm_ = TMath::CosH(deltaGamma*abscissa/2.0)*expTerm_;
 		sinhTerm_ = TMath::SinH(deltaGamma*abscissa/2.0)*expTerm_;
 
 	}
 
 	if ( effiFun_ ) {
 		this->updateEffiSpline(effiPars_);
 		effiTerm_ = effiFun_->evaluate(abscissa);
 		if (effiTerm_>1.0){effiTerm_=1.0;}
 		if (effiTerm_<0.0){effiTerm_=0.0;}
 	} else {
 		effiTerm_ = 1.0;
 	}
 }
 
 Double_t LauDecayTimePdf::normExpHypTerm(Double_t Abs)
 {
 	Double_t tau         = tau_->value();
 	Double_t deltaGamma  = deltaGamma_->value();
 
 	Double_t y           = tau*deltaGamma/2;
 	Double_t nonTrigTerm = -(TMath::Exp(-Abs/tau))/(1 - y*y);
 
 	Double_t cosHTerm    = TMath::CosH(deltaGamma*Abs/2);
 	Double_t sinHTerm    = TMath::SinH(deltaGamma*Abs/2);
 
 	Double_t normTerm =  nonTrigTerm*(cosHTerm + y*sinHTerm);
 	return normTerm;
 }
 
 Double_t LauDecayTimePdf::normExpHypTermDep(Double_t Abs)
 {
 	Double_t tau         = tau_->value();
 	Double_t deltaGamma  = deltaGamma_->value();
 
 	Double_t y           = tau*deltaGamma/2;
 	Double_t nonTrigTerm = -(TMath::Exp(-Abs/tau))/(1 - y*y);
 
 	Double_t cosHTerm    = TMath::CosH(deltaGamma*Abs/2);
 	Double_t sinHTerm    = TMath::SinH(deltaGamma*Abs/2);
 
 	Double_t normTerm =  nonTrigTerm*(sinHTerm + y*cosHTerm);
 	return normTerm;
 }
 
 void LauDecayTimePdf::storeIntegralTimeDep(Double_t minAbs, Double_t maxAbs)
 {
 	Double_t tau = tau_->value();
 	// Normalisation factor for B0 decays
 	if (type_ == ExpTrig || type_  == SimFitNormBd || type_ == SimFitSigBd ) {
                 if (method_ == DecayTime) {
 		        normTermExp_ = tau * (TMath::Exp(-minAbs/tau) - TMath::Exp(-maxAbs/tau));
                 }
                 if (method_ == DecayTimeDiff) {
 		        normTermExp_ = tau * (2.0 - TMath::Exp(-maxAbs/tau) - TMath::Exp(-minAbs/tau));
                 }
 	}
 
 	// Normalisation factor for Bs decays
 	if (type_ == ExpHypTrig) {
 		normTermCosh_ = normExpHypTerm(maxAbs) - normExpHypTerm(minAbs);
 		normTermSinh_ = normExpHypTermDep( maxAbs) - normExpHypTermDep( minAbs);
 	}
 }
 
 Double_t LauDecayTimePdf::generateError(Bool_t forceNew)
 {
 	if (errHist_ && (forceNew || !abscissaErrorGenerated_)) {
 		LauFitData errData = errHist_->generate(0);
 		abscissaError_ = errData.find(this->varErrName())->second;
 		abscissaErrorGenerated_ = kTRUE;
 	} else {
 		while (forceNew || !abscissaErrorGenerated_) {
 			abscissaError_ = LauRandom::randomFun()->Landau(errorDistMPV_,errorDistSigma_);
 			if (abscissaError_ < maxAbscissaError_ && abscissaError_ > minAbscissaError_) {
 				abscissaErrorGenerated_ = kTRUE;
 				forceNew = kFALSE;
 			}
 		}
 	}
 	return abscissaError_;
 }
 
 /*
 LauFitData LauDecayTimePdf::generate(const LauKinematics* kinematics)
 {
 	// generateError SHOULD have been called before this
 	// function but will call it here just to make sure
 	// (has ns effect if has already been called)
 	abscissaError_ = this->generateError();
 
 	// If the PDF is scaled by the per-event error then need to update the PDF height for each event
 	Bool_t scale(kFALSE);
 	for (std::vector<Bool_t>::const_iterator iter = scaleMeans_.begin(); iter != scaleMeans_.end(); ++iter) {
 		scale |= (*iter);
 	}
 	for (std::vector<Bool_t>::const_iterator iter = scaleWidths_.begin(); iter != scaleWidths_.end(); ++iter) {
 		scale |= (*iter);
 	}
 	if (scale || (!this->heightUpToDate() && !this->cachePDF())) {
 		this->calcPDFHeight(kinematics);
 		this->heightUpToDate(kTRUE);
 	}
 
 	// Generate the value of the abscissa.
 	Bool_t gotAbscissa(kFALSE);
 
 	Double_t genVal(0.0);
 	Double_t genPDFVal(0.0);
 	LauFitData genAbscissa;
 
 	const Double_t xMin = this->minAbscissa();
 	const Double_t xMax = this->maxAbscissa();
 	const Double_t xRange = xMax - xMin;
 
 	while (!gotAbscissa) {
 
 		genVal = LauRandom::randomFun()->Rndm()*xRange + xMin;
 		this->calcLikelihoodInfo(genVal, abscissaError_);
 		genPDFVal = this->getUnNormLikelihood();
 
 		if (LauRandom::randomFun()->Rndm() <= genPDFVal/this->getMaxHeight()) {gotAbscissa = kTRUE;}
 
 		if (genPDFVal > this->getMaxHeight()) {
 			cerr<<"Warning in LauDecayTimePdf::generate()."
 				<<" genPDFVal = "<<genPDFVal<<" is larger than the specified PDF height "
 				<<this->getMaxHeight()<<" for the abscissa = "<<genVal<<". Need to reset height to be larger than "
 				<<genPDFVal<<" by using the setMaxHeight(Double_t) function"
 				<<" and re-run the Monte Carlo generation!"<<endl;
 		}
 	}
 
 	genAbscissa[this->varName()] = genVal;
 
 	// mark that we need a new error to be generated next time
 	abscissaErrorGenerated_ = kFALSE;
 
 	return genAbscissa;
 }
 */
 
 void LauDecayTimePdf::setErrorHisto(const TH1* hist)
 {
 	if ( errHist_ != 0 ) {
 		cerr<<"WARNING in LauDecayTimePdf::setErrorHisto : Error histogram already set, not doing it again."<<endl;
 		return;
 	}
 	errHist_ = new Lau1DHistPdf(this->varErrName(), hist, this->minAbscissaError(), this->maxAbscissaError());
 }
 
 void LauDecayTimePdf::setHistoPdf(const TH1* hist)
 {
 	if ( pdfHist_ != 0 ) {
 		cerr<<"WARNING in LauDecayTimePdf::setHistoPdf : PDF histogram already set, not doing it again."<<endl;
 		return;
 	}
 	pdfHist_ = new Lau1DHistPdf(this->varName(), hist, this->minAbscissa(), this->maxAbscissa());
 }
 
 void LauDecayTimePdf::setEffiSpline(Lau1DCubicSpline* spline)
 {
         if ( effiFun_ != 0 ) {
                 cerr<<"WARNING in LauDecayTimePdf::setEffiPdf : efficiency function already set, not doing it again."<<endl;
                 return;
         }
 	effiFun_ = spline;
 }
 
 void LauDecayTimePdf::setEffiSpline(Lau1DCubicSpline* spline, std::vector<LauParameter*> effiPars)
 {
 	if ( effiFun_ != 0 ) {
 		cerr<<"WARNING in LauDecayTimePdf::setEffiPdf : efficiency function already set, not doing it again."<<endl;
 		return;
 	}
 	effiFun_ = spline;
 
 	if (effiPars.size() != spline->getnKnots()){
 		cerr<<"ERROR in LauDecayTimePdf::setEffiPdf : number of efficiency parameters is not equal to the number of spline knots."<<endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	effiPars_ = effiPars;
 }
 
 LauAbsRValue* LauDecayTimePdf::findParameter(const TString& parName)
 {
 	for ( std::vector<LauAbsRValue*>::iterator iter = param_.begin(); iter != param_.end(); ++iter ) {
 		if ((*iter)->name().Contains(parName)) {
 			return (*iter);
 		}
 	}
 	std::cerr << "ERROR in LauDecayTimePdf::findParameter : Parameter \"" << parName << "\" not found." << std::endl;
 	return 0;
 }
 
 const LauAbsRValue* LauDecayTimePdf::findParameter(const TString& parName) const
 {
 	for ( std::vector<LauAbsRValue*>::const_iterator iter = param_.begin(); iter != param_.end(); ++iter ) {
 		if ((*iter)->name().Contains(parName)) {
 			return (*iter);
 		}
 	}
 	std::cerr << "ERROR in LauDecayTimePdf::findParameter : Parameter \"" << parName << "\" not found." << std::endl;
 	return 0;
 }
 
 void LauDecayTimePdf::updatePulls()
 {
 	for ( std::vector<LauAbsRValue*>::iterator iter = param_.begin(); iter != param_.end(); ++iter ) {
 		std::vector<LauParameter*> params = (*iter)->getPars();
 		for (std::vector<LauParameter*>::iterator params_iter = params.begin(); params_iter != params.end(); ++params_iter ) {
 			if (!(*iter)->fixed()) {
 				(*params_iter)->updatePull();
 			}
 		}
 	}
 }
 
-void LauDecayTimePdf::updateEffiSpline(std::vector<LauParameter*> params)
+void LauDecayTimePdf::updateEffiSpline(std::vector<LauParameter*> effiPars)
 {
-	if (!params.empty()){
-		std::vector<Double_t> yvals;
-		yvals.reserve(params.size());
-		for (ULong_t i(0); i<params.size(); ++i){
-			yvals.push_back(params[i]->unblindValue());
-		}
-		effiFun_->updateYValues(yvals);
+	if (effiPars.size() != spline->getnKnots()){
+		cerr<<"ERROR in LauDecayTimePdf::updateEffiSpline : number of efficiency parameters is not equal to the number of spline knots."<<endl;
+		gSystem->Exit(EXIT_FAILURE);
 	}
+	effiFun_->updateYValues(effiPars);
 }
 
 void LauDecayTimePdf::setAsymmetries(LauParameter* AProd, LauParameter* ARaw)
 {
 	AProd_ = AProd;
 	ARaw_ = ARaw;
 }
 
 void LauDecayTimePdf::setFlavourTagging(LauFlavTag* flavTag)
 {
 	flavTag_ = flavTag;
 }
 
diff --git a/src/LauTimeDepFitModel.cc b/src/LauTimeDepFitModel.cc
index 075c331..0f09839 100644
--- a/src/LauTimeDepFitModel.cc
+++ b/src/LauTimeDepFitModel.cc
@@ -1,2478 +1,2479 @@
 
 /*
 Copyright 2006 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
 /*! \file LauTimeDepFitModel.cc
     \brief File containing implementation of LauTimeDepFitModel class.
 */
 
 #include <iostream>
 #include <iomanip>
 #include <fstream>
 #include <map>
 #include <vector>
 
 #include "TFile.h"
 #include "TMinuit.h"
 #include "TRandom.h"
 #include "TSystem.h"
 #include "TVirtualFitter.h"
 
 #include "LauAbsBkgndDPModel.hh"
 #include "LauAbsCoeffSet.hh"
 #include "LauAbsPdf.hh"
 #include "LauAsymmCalc.hh"
 #include "LauComplex.hh"
 #include "LauConstants.hh"
 #include "LauDPPartialIntegralInfo.hh"
 #include "LauDaughters.hh"
 #include "LauDecayTimePdf.hh"
 #include "LauFitNtuple.hh"
 #include "LauGenNtuple.hh"
 #include "LauIsobarDynamics.hh"
 #include "LauKinematics.hh"
 #include "LauPrint.hh"
 #include "LauRandom.hh"
 #include "LauScfMap.hh"
 #include "LauTimeDepFitModel.hh"
 #include "LauFlavTag.hh"
 
 
 ClassImp(LauTimeDepFitModel)
 
 LauTimeDepFitModel::LauTimeDepFitModel(LauIsobarDynamics* modelB0bar, LauIsobarDynamics* modelB0, LauFlavTag* flavTag) : LauAbsFitModel(),
 	sigModelB0bar_(modelB0bar),
 	sigModelB0_(modelB0),
 	kinematicsB0bar_(modelB0bar ? modelB0bar->getKinematics() : 0),
 	kinematicsB0_(modelB0 ? modelB0->getKinematics() : 0),
 	flavTag_(flavTag),
 	nSigComp_(0),
 	nSigDPPar_(0),
 	nDecayTimePar_(0),
 	nExtraPdfPar_(0),
 	nNormPar_(0),
 	nCalibPar_(0),
 	nEffiPar_(0),
 	nAsymPar_(0),
 	coeffsB0bar_(0),
 	coeffsB0_(0),
 	coeffPars_(0),
 	fitFracB0bar_(0),
 	fitFracB0_(0),
 	fitFracAsymm_(0),
 	acp_(0),
 	meanEffB0bar_("meanEffB0bar",0.0,0.0,1.0),
 	meanEffB0_("meanEffB0",0.0,0.0,1.0),
 	DPRateB0bar_("DPRateB0bar",0.0,0.0,100.0),
 	DPRateB0_("DPRateB0",0.0,0.0,100.0),
 	signalEvents_(0),
 	signalAsym_(0),
 	cpevVarName_(""),
 	cpEigenValue_(CPEven),
 	evtCPEigenVals_(0),
 	deltaM_("deltaM",0.0),
 	deltaGamma_("deltaGamma",0.0),
 	tau_("tau",LauConstants::tauB0),
 	phiMix_("phiMix", 2.0*LauConstants::beta, -LauConstants::threePi, LauConstants::threePi, kFALSE),
 	sinPhiMix_("sinPhiMix", TMath::Sin(2.0*LauConstants::beta), -1.0, 1.0, kFALSE),
 	cosPhiMix_("cosPhiMix", TMath::Cos(2.0*LauConstants::beta), -1.0, 1.0, kFALSE),
 	useSinCos_(kFALSE),
 	phiMixComplex_(TMath::Cos(-2.0*LauConstants::beta),TMath::Sin(-2.0*LauConstants::beta)),
 	signalDecayTimePdfs_(),
 	curEvtDecayTime_(0.0),
 	curEvtDecayTimeErr_(0.0),
 	sigExtraPdf_(),
 	sigFlavTagPdf_(),
 	bkgdFlavTagPdf_(),
 	AProd_("AProd",0.0,-1.0,1.0,kTRUE),
 	ARaw_("ARaw",1.0,0.0,2.0,kTRUE),
 	iterationsMax_(500000),
 	nGenLoop_(0),
 	ASq_(0.0),
 	aSqMaxVar_(0.0),
 	aSqMaxSet_(1.25),
 	storeGenAmpInfo_(kFALSE),
 	signalTree_(),
 	reuseSignal_(kFALSE),
 	sigDPLike_(0.0),
 	sigExtraLike_(0.0),
 	sigFlavTagLike_(0.0),
 	bkgdFlavTagLike_(0.0),
 	sigTotalLike_(0.0)
 {
 	// Set up ftag here?
 
 	// Make sure that the integration scheme will be symmetrised
 	sigModelB0bar_->forceSymmetriseIntegration(kTRUE);
 	sigModelB0_->forceSymmetriseIntegration(kTRUE);
 }
 
 LauTimeDepFitModel::~LauTimeDepFitModel()
 {
    for (LauTagCatEmbDataMap::iterator iter = signalTree_.begin(); iter != signalTree_.end(); ++iter){
       delete iter->second;
    }
  
    for (LauTagCatEmbDataMapList::iterator iterlist = bkgndTree_.begin(); iterlist != bkgndTree_.end(); ++iterlist){
        for (LauTagCatEmbDataMap::iterator iter = (*iterlist).begin(); iter != (*iterlist).end(); ++iter){
            delete iter->second;
        }
    }
 
 }
 
 void LauTimeDepFitModel::setupBkgndVectors()
 {
         UInt_t nBkgnds = this->nBkgndClasses();
         BkgndDPModelsB0_.resize( nBkgnds );
         BkgndDPModelsB0bar_.resize( nBkgnds );
         BkgndPdfsB0_.resize( nBkgnds );
         BkgndPdfsB0bar_.resize( nBkgnds );
         bkgndEvents_.resize( nBkgnds );
         bkgndAsym_.resize( nBkgnds );
         bkgndTree_.resize( nBkgnds );
         reuseBkgnd_.resize( nBkgnds );
         bkgndDPLike_.resize( nBkgnds );
         bkgndExtraLike_.resize( nBkgnds );
         bkgndTotalLike_.resize( nBkgnds );
 }
 
 void LauTimeDepFitModel::setNSigEvents(LauParameter* nSigEvents)
 {
 	if ( nSigEvents == 0 ) {
 		std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The LauParameter pointer is null." << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	if ( signalEvents_ != 0 ) {
 		std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
 		return;
 	}
 	if ( signalAsym_ != 0 ) {
 		std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl;
 		return;
 	}
 
 	signalEvents_ = nSigEvents;
 
 	signalEvents_->name("signalEvents");
 
 	Double_t value = nSigEvents->value();
 	signalEvents_->range(-2.0*(TMath::Abs(value)+1.0),2.0*(TMath::Abs(value)+1.0));
 
 	signalAsym_ = new LauParameter("signalAsym",0.0,-1.0,1.0,kTRUE);
 }
 
 void LauTimeDepFitModel::setNSigEvents(LauParameter* nSigEvents, LauParameter* sigAsym)
 {
 	if ( nSigEvents == 0 ) {
 		std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The event LauParameter pointer is null." << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	if ( sigAsym == 0 ) {
 		std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The asym LauParameter pointer is null." << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	if ( signalEvents_ != 0 ) {
 		std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
 		return;
 	}
 	if ( signalAsym_ != 0 ) {
 		std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl;
 		return;
 	}
 
 	signalEvents_ = nSigEvents;
 	signalEvents_->name("signalEvents");
 	Double_t value = nSigEvents->value();
 	signalEvents_->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
 
 	signalAsym_ = sigAsym;
 	signalAsym_->name("signalAsym");
 	signalAsym_->range(-1.0,1.0);
 }
 
 void LauTimeDepFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents)
 {
          if ( nBkgndEvents == 0 ) {
                 std::cerr << "ERROR in LauTimeDepFitModel::setNBgkndEvents : The background yield LauParameter pointer is null." << std::endl;
                 gSystem->Exit(EXIT_FAILURE);
         } 
 
         if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
                 std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl;
                 std::cerr << "                                        : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
                 gSystem->Exit(EXIT_FAILURE);
         }
 
         UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() );
 
         if ( bkgndEvents_[bkgndID] != 0 ) {
                 std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl;
                 return;
         }
 
         if ( bkgndAsym_[bkgndID] != 0 ) {
                 std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl;
                 return;
         }
 
         nBkgndEvents->name( nBkgndEvents->name()+"Events" );
 	if ( nBkgndEvents->isLValue() ) {
 		Double_t value = nBkgndEvents->value();
 		LauParameter* yield = dynamic_cast<LauParameter*>( nBkgndEvents );
 		yield->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
 	}
         bkgndEvents_[bkgndID] = nBkgndEvents;
 
         bkgndAsym_[bkgndID] = new LauParameter(nBkgndEvents->name()+"Asym",0.0,-1.0,1.0,kTRUE);
 }
 
 void LauTimeDepFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym)
 {
         if ( nBkgndEvents == 0 ) {
                 std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : The background yield LauParameter pointer is null." << std::endl;
                 gSystem->Exit(EXIT_FAILURE);
         }
 
         if ( bkgndAsym == 0 ) {
                 std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : The background asym LauParameter pointer is null." << std::endl;
                 gSystem->Exit(EXIT_FAILURE);
         }
 
         if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
                 std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl;
                 std::cerr << "                                        : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
                 gSystem->Exit(EXIT_FAILURE);
         }
 
         UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() );
 
         if ( bkgndEvents_[bkgndID] != 0 ) {
                 std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl;
                 return;
         }
 
         if ( bkgndAsym_[bkgndID] != 0 ) {
                 std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : You are trying to overwrite the background asymmetry." << std::endl;
                 return;
         }
 
         bkgndEvents_[bkgndID]->name( nBkgndEvents->name()+"Events" );
 	if ( nBkgndEvents->isLValue() ) {
 		Double_t value = nBkgndEvents->value();
 		LauParameter* yield = dynamic_cast<LauParameter*>( nBkgndEvents );
 		yield->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
 	}
         bkgndEvents_[bkgndID] = nBkgndEvents;
 
         bkgndAsym_[bkgndID]->name( nBkgndEvents->name()+"Asym" );
 	if ( bkgndAsym->isLValue() ) {
 		LauParameter* asym = dynamic_cast<LauParameter*>( bkgndAsym );
 		asym->range(-1.0, 1.0);
 	}
         bkgndAsym_[bkgndID] = bkgndAsym;
 }
 
 
 void LauTimeDepFitModel::setSignalDtPdf(Int_t tagCat, LauDecayTimePdf* pdf)
 {
 	if (!flavTag_->validTagCat(tagCat)) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::setSignalDtPdf : Tagging category \""<<tagCat<<"\" not valid, not adding PDF."<<std::endl;
 		return;
 	}
 	if (pdf==0) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::setSignalDtPdf : The PDF pointer is null, not adding it."<<std::endl;
 		return;
 	}
 	signalDecayTimePdfs_[tagCat] = pdf;
 	signalDecayTimePdfs_[tagCat]->setAsymmetries(&AProd_, &ARaw_);
 }
 
 void LauTimeDepFitModel::setSignalPdfs(Int_t tagCat, LauAbsPdf* pdf)
 {
 	// These "extra variables" are assumed to be purely kinematical, like mES and DeltaE
 	//or making use of Rest of Event information, and therefore independent of whether
 	//the parent is a B0 or a B0bar. If this assupmtion doesn't hold, do modify this part!
 	if (!flavTag_->validTagCat(tagCat)) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::setSignalPdfs : Tagging category \""<<tagCat<<"\" not valid, not adding models."<<std::endl;
 		return;
 	}
 	if (pdf==0) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::setSignalPdfs : The PDF pointer is null."<<std::endl;
 		return;
 	}
 
 	sigExtraPdf_[tagCat].push_back(pdf);
 }
 
 void LauTimeDepFitModel::setPhiMix(const Double_t phiMix, const Bool_t fixPhiMix, const Bool_t useSinCos)
 {
 	phiMix_.value(phiMix);  phiMix_.initValue(phiMix);  phiMix_.genValue(phiMix);  phiMix_.fixed(fixPhiMix);
 
 	const Double_t sinPhiMix = TMath::Sin(phiMix);
 	sinPhiMix_.value(sinPhiMix); sinPhiMix_.initValue(sinPhiMix); sinPhiMix_.genValue(sinPhiMix); sinPhiMix_.fixed(fixPhiMix);
 
 	const Double_t cosPhiMix = TMath::Cos(phiMix);
 	cosPhiMix_.value(cosPhiMix); cosPhiMix_.initValue(cosPhiMix); cosPhiMix_.genValue(cosPhiMix); cosPhiMix_.fixed(fixPhiMix);
 
 	useSinCos_ = useSinCos;
 
 	phiMixComplex_.setRealPart(cosPhiMix);
 	phiMixComplex_.setImagPart(-1.0*sinPhiMix);
 }
 
 void LauTimeDepFitModel::initialise()
 {
 	// From the initial parameter values calculate the coefficients
 	// so they can be passed to the signal model
 	this->updateCoeffs();
 
 	// Initialisation
 	if (this->useDP() == kTRUE) {
 		this->initialiseDPModels();
 	}
 
 	//Flavour tagging
 	flavTag_->initialise();
 
 	if (!this->useDP() && sigExtraPdf_.empty()) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::initialise : Signal model doesn't exist for any variable."<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	if (this->useDP() == kTRUE) {
 		// Check that we have all the Dalitz-plot models
 		if ((sigModelB0bar_ == 0) || (sigModelB0_ == 0)) {
 			std::cerr<<"ERROR in LauTimeDepFitModel::initialise : the pointer to one (particle or anti-particle) of the signal DP models is null."<<std::endl;
 			gSystem->Exit(EXIT_FAILURE);
 		}
 	}
 
 	// Clear the vectors of parameter information so we can start from scratch
 	this->clearFitParVectors();
 
 	// Set the fit parameters for signal and background models
 	this->setSignalDPParameters();
 
 	// Set the fit parameters for the decay time models
 	this->setDecayTimeParameters();
 
 	// Set the fit parameters for the extra PDFs
 	this->setExtraPdfParameters();
 
 	// Set the initial bg and signal events
 	this->setFitNEvents();
 
 	// Handle flavour-tagging calibration parameters
 	this->setCalibParams();
 
 	// Add the efficiency parameters
 	this->setEffiParams();
 
 	//Asymmetry terms AProd and ARaw added in setAsymmetries()?
 	this->setAsymParams();
 
 	// Check that we have the expected number of fit variables
 	const LauParameterPList& fitVars = this->fitPars();
 	if (fitVars.size() != (nSigDPPar_ + nDecayTimePar_ + nExtraPdfPar_ + nNormPar_ + nCalibPar_ + nEffiPar_ + nAsymPar_)) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::initialise : Number of fit parameters not of expected size."<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	this->setExtraNtupleVars();
 }
 
 void LauTimeDepFitModel::recalculateNormalisation()
 {
 	sigModelB0bar_->recalculateNormalisation();
 	sigModelB0_->recalculateNormalisation();
 	sigModelB0bar_->modifyDataTree();
 	sigModelB0_->modifyDataTree();
 	this->calcInterferenceTermIntegrals();
 }
 
 void LauTimeDepFitModel::initialiseDPModels()
 {
 	if (sigModelB0bar_ == 0) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::initialiseDPModels : B0bar signal DP model doesn't exist"<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	if (sigModelB0_ == 0) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::initialiseDPModels : B0 signal DP model doesn't exist"<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	// Need to check that the number of components we have and that the dynamics has matches up
 	const UInt_t nAmpB0bar = sigModelB0bar_->getnTotAmp();
 	const UInt_t nAmpB0 = sigModelB0_->getnTotAmp();
 	if ( nAmpB0bar != nAmpB0 ) {
 		std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Unequal number of signal DP components in the particle and anti-particle models: " << nAmpB0bar << " != " << nAmpB0 << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	if ( nAmpB0bar != nSigComp_ ) {
 		std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Number of signal DP components in the model (" << nAmpB0bar << ") not equal to number of coefficients supplied (" << nSigComp_ << ")." << std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 
 	std::cout<<"INFO in LauTimeDepFitModel::initialiseDPModels : Initialising signal DP model"<<std::endl;
 	sigModelB0bar_->initialise(coeffsB0bar_);
 	sigModelB0_->initialise(coeffsB0_);
 
 	fifjEffSum_.clear();
 	fifjEffSum_.resize(nSigComp_);
 	for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) {
 		fifjEffSum_[iAmp].resize(nSigComp_);
 	}
 
 	// calculate the integrals of the A*Abar terms
 	this->calcInterferenceTermIntegrals();
 
 	this->calcInterTermNorm();
 
 }
 
 void LauTimeDepFitModel::calcInterferenceTermIntegrals()
 {
 	const std::vector<LauDPPartialIntegralInfo*>& integralInfoListB0bar = sigModelB0bar_->getIntegralInfos();
 	const std::vector<LauDPPartialIntegralInfo*>& integralInfoListB0 = sigModelB0_->getIntegralInfos();
 
 	// TODO should check (first time) that they match in terms of number of entries in the vectors and that each entry has the same number of points, ranges, weights etc.
 
 	LauComplex A, Abar, fifjEffSumTerm;
 
 	for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) {
 		for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) {
 			fifjEffSum_[iAmp][jAmp].zero();
 		}
 	}
 
 
 	const UInt_t nIntegralRegions = integralInfoListB0bar.size();
 	for ( UInt_t iRegion(0); iRegion < nIntegralRegions; ++iRegion ) {
 
 		const LauDPPartialIntegralInfo* integralInfoB0bar = integralInfoListB0bar[iRegion];
 		const LauDPPartialIntegralInfo* integralInfoB0 = integralInfoListB0[iRegion];
 
 		const UInt_t nm13Points = integralInfoB0bar->getnm13Points();
 		const UInt_t nm23Points = integralInfoB0bar->getnm23Points();
 
 		for (UInt_t m13 = 0; m13 < nm13Points; ++m13) {
 			for (UInt_t m23 = 0; m23 < nm23Points; ++m23) {
 
 				const Double_t weight = integralInfoB0bar->getWeight(m13,m23);
 				const Double_t eff = integralInfoB0bar->getEfficiency(m13,m23);
 				const Double_t effWeight = eff*weight;
 
 				for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) {
 
 					A = integralInfoB0->getAmplitude(m13, m23, iAmp);
 
 					for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) {
 
 						Abar = integralInfoB0bar->getAmplitude(m13, m23, jAmp);
 
 						fifjEffSumTerm = Abar*A.conj();
 
 						fifjEffSumTerm.rescale(effWeight);
 						fifjEffSum_[iAmp][jAmp] += fifjEffSumTerm;
 					}
 				}
 
 			}
 		}
 	}
 
 }
 
 void LauTimeDepFitModel::calcInterTermNorm()
 {
 	const std::vector<Double_t>& fNormB0bar = sigModelB0bar_->getFNorm();
 	const std::vector<Double_t>& fNormB0 = sigModelB0_->getFNorm();
 
 	LauComplex norm;
 	for (UInt_t iAmp = 0; iAmp < nSigComp_; ++iAmp) {
 		for (UInt_t jAmp = 0; jAmp < nSigComp_; ++jAmp) {
 
 			LauComplex coeffTerm = coeffsB0bar_[jAmp]*coeffsB0_[iAmp].conj();
 			coeffTerm *= fifjEffSum_[iAmp][jAmp];
 			coeffTerm.rescale(fNormB0bar[jAmp] * fNormB0[iAmp]);
 
 			norm += coeffTerm;
 		}
 	}
 	norm *= phiMixComplex_;
 
 	interTermReNorm_ = 2.0*norm.re();
 	interTermImNorm_ = 2.0*norm.im();
 }
 
 void LauTimeDepFitModel::setAmpCoeffSet(LauAbsCoeffSet* coeffSet)
 {
 	// Is there a component called compName in the signal models?
 	TString compName = coeffSet->name();
 	TString conjName = sigModelB0bar_->getConjResName(compName);
 
 	const LauDaughters* daughtersB0bar = sigModelB0bar_->getDaughters();
 	const LauDaughters* daughtersB0 = sigModelB0_->getDaughters();
 	const Bool_t conjugate = daughtersB0bar->isConjugate( daughtersB0 );
 
 	if ( ! sigModelB0bar_->hasResonance(compName) ) {
 		if ( ! sigModelB0bar_->hasResonance(conjName) ) {
 			std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0bar signal DP model doesn't contain component \""<<compName<<"\"."<<std::endl;
 			return;
 		}
 		std::cerr<<"WARNING in LauTimeDepFitModel::setAmpCoeffSet : B0bar signal DP model doesn't contain component \""<<compName<<"\" but does contain the conjugate \""<<conjName<<"\", resetting name to use the conjugate."<<std::endl;
 		TString tmp = compName;
 		compName = conjName;
 		conjName = tmp;
 		coeffSet->name( compName );
 	}
 
 	if ( conjugate ) {
 		if ( ! sigModelB0_->hasResonance(conjName) ) {
 			std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<<conjName<<"\"."<<std::endl;
 			return;
 		}
 	} else {
 		if ( ! sigModelB0_->hasResonance(compName) ) {
 			std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<<compName<<"\"."<<std::endl;
 			return;
 		}
 	}
 
 	// Do we already have it in our list of names?
 	for (std::vector<LauAbsCoeffSet*>::const_iterator iter=coeffPars_.begin(); iter!=coeffPars_.end(); ++iter) {
 		if ((*iter)->name() == compName) {
 			std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : Have already set coefficients for \""<<compName<<"\"."<<std::endl;
 			return;
 		}
 	}
 
 	coeffSet->index(nSigComp_);
 	coeffPars_.push_back(coeffSet);
 
 	TString parName = coeffSet->baseName(); parName += "FitFracAsym";
 	fitFracAsymm_.push_back(LauParameter(parName, 0.0, -1.0, 1.0));
 
 	acp_.push_back(coeffSet->acp());
 
 	++nSigComp_;
 
 	std::cout<<"INFO in LauTimeDepFitModel::setAmpCoeffSet : Added coefficients for component \""<<compName<<"\" to the fit model."<<std::endl;
 }
 
 void LauTimeDepFitModel::calcAsymmetries(Bool_t initValues)
 {
 	// Calculate the CP asymmetries
 	// Also calculate the fit fraction asymmetries
 
 	for (UInt_t i = 0; i < nSigComp_; i++) {
 
 		acp_[i] = coeffPars_[i]->acp();
 
 		LauAsymmCalc asymmCalc(fitFracB0bar_[i][i].value(), fitFracB0_[i][i].value());
 		Double_t asym = asymmCalc.getAsymmetry();
 		fitFracAsymm_[i].value(asym);
 		if (initValues) {
 			fitFracAsymm_[i].genValue(asym);
 			fitFracAsymm_[i].initValue(asym);
 		}
 	}
 
 }
 
 void LauTimeDepFitModel::setSignalDPParameters()
 {
 	// Set the fit parameters for the signal model.
 	nSigDPPar_ = 0;
 	if ( ! this->useDP() ) {
 		return;
 	}
 
 	std::cout << "INFO in LauTimeDepFitModel::setSignalDPParameters : Setting the initial fit parameters for the signal DP model." << std::endl;
 
 	// Place isobar coefficient parameters in vector of fit variables
 	LauParameterPList& fitVars = this->fitPars();
 	for (UInt_t i = 0; i < nSigComp_; ++i) {
 		LauParameterPList pars = coeffPars_[i]->getParameters();
 		for (LauParameterPList::iterator iter = pars.begin(); iter != pars.end(); ++iter) {
 			if ( !(*iter)->clone() ) {
 				fitVars.push_back(*iter);
 				++nSigDPPar_;
 			}
 		}
 	}
 
 	// Obtain the resonance parameters and place them in the vector of fit variables and in a separate vector
 	// Need to make sure that they are unique because some might appear in both DP models
 	LauParameterPSet& resVars = this->resPars();
 	resVars.clear();
 
 	LauParameterPList& sigDPParsB0bar = sigModelB0bar_->getFloatingParameters();
 	LauParameterPList& sigDPParsB0 = sigModelB0_->getFloatingParameters();
 
 	for ( LauParameterPList::iterator iter = sigDPParsB0bar.begin(); iter != sigDPParsB0bar.end(); ++iter ) {
 		if ( resVars.insert(*iter).second ) {
 			fitVars.push_back(*iter);
 			++nSigDPPar_;
 		}
 	}
 	for ( LauParameterPList::iterator iter = sigDPParsB0.begin(); iter != sigDPParsB0.end(); ++iter ) {
 		if ( resVars.insert(*iter).second ) {
 			fitVars.push_back(*iter);
 			++nSigDPPar_;
 		}
 	}
 }
 
 UInt_t LauTimeDepFitModel::addParametersToFitList(LauTagCatDtPdfMap& theMap)
 {
 	UInt_t counter(0);
 	LauParameterPList& fitVars = this->fitPars();
 	// loop through the map
 	for (LauTagCatDtPdfMap::iterator iter = theMap.begin(); iter != theMap.end(); ++iter) {
 		// grab the pdf and then its parameters
 		LauDecayTimePdf* thePdf = (*iter).second;	// The first one is the tagging category
 		LauAbsRValuePList& rvalues = thePdf->getParameters();
 		// loop through the parameters
 		for (LauAbsRValuePList::iterator pars_iter = rvalues.begin(); pars_iter != rvalues.end(); ++pars_iter) {
 			LauParameterPList params = (*pars_iter)->getPars();
 			for (LauParameterPList::iterator params_iter = params.begin(); params_iter != params.end(); ++params_iter) {
 				// for each "original" parameter add it to the list of fit parameters and increment the counter
 				if ( !(*params_iter)->clone() && ( !(*params_iter)->fixed() ||
 							(this->twoStageFit() && (*params_iter)->secondStage()) ) ) {
 					fitVars.push_back(*params_iter);
 					++counter;
 				}
 			}
 		}
 	}
 	return counter;
 }
 
 UInt_t LauTimeDepFitModel::addParametersToFitList(LauTagCatPdfListMap& theMap)
 {
 	UInt_t counter(0);
 	// loop through the map
 	for (LauTagCatPdfListMap::iterator iter = theMap.begin(); iter != theMap.end(); ++iter) {
 		counter += this->addFitParameters(iter->second); // first is the tagging category
 	}
 	return counter;
 }
 
 void LauTimeDepFitModel::setDecayTimeParameters()
 {
 	nDecayTimePar_ = 0;
 
 	// Loop over the Dt PDFs
 	nDecayTimePar_ += this->addParametersToFitList(signalDecayTimePdfs_);
 
 	LauParameterPList& fitVars = this->fitPars();
 	if (useSinCos_) {
 		fitVars.push_back(&sinPhiMix_);
 		fitVars.push_back(&cosPhiMix_);
 		nDecayTimePar_ += 2;
 	} else {
 		fitVars.push_back(&phiMix_);
 		++nDecayTimePar_;
 	}
 }
 
 void LauTimeDepFitModel::setExtraPdfParameters()
 {
 	// Include the parameters of the PDF for each tagging category in the fit
 	// NB all of them are passed to the fit, even though some have been fixed through parameter.fixed(kTRUE)
 	// With the new "cloned parameter" scheme only "original" parameters are passed to the fit.
 	// Their clones are updated automatically when the originals are updated.
 
 	nExtraPdfPar_ = 0;
 
 	nExtraPdfPar_ += this->addParametersToFitList(sigExtraPdf_);
 }
 
 void LauTimeDepFitModel::setFitNEvents()
 {
 	nNormPar_ = 0;
 
 	// Initialise the total number of events to be the sum of all the hypotheses
 	Double_t nTotEvts = signalEvents_->value();
 
 	this->eventsPerExpt(TMath::FloorNint(nTotEvts));
 
 	LauParameterPList& fitVars = this->fitPars();
 
 	// if doing an extended ML fit add the signal fraction into the fit parameters
 	if (this->doEMLFit()) {
 		std::cout<<"INFO in LauTimeDepFitModel::setFitNEvents : Initialising number of events for signal and background components..."<<std::endl;
 		fitVars.push_back(signalEvents_);
 		++nNormPar_;
 	} else {
 		std::cout<<"INFO in LauTimeDepFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..."<<std::endl;
 	}
 
 	// if not using the DP in the model we need an explicit signal asymmetry parameter
 	if (this->useDP() == kFALSE) {
 		fitVars.push_back(signalAsym_);
 		++nNormPar_;
 	}
 
 	// TODO arguably should delegate this
 	LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac();
 
 	// tagging-category fractions for signal events
 	for (LauTagCatParamMap::iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) {
 		if (iter == signalTagCatFrac.begin()) {
 			continue;
 		}
 		LauParameter* par = &((*iter).second);
 		fitVars.push_back(par);
 		++nNormPar_;
 	}
 }
 
 void LauTimeDepFitModel::setAsymParams()
 {
 	LauParameterPList& fitVars = this->fitPars();
 	fitVars.push_back(&AProd_);
 	fitVars.push_back(&ARaw_);
 	nAsymPar_+=2;
 }
 
 void LauTimeDepFitModel::setCalibParams()
 {
 	Bool_t useAltPars = flavTag_->getUseAveOmegaDeltaOmega();
 
 	if (useAltPars){
 		LauTagCatParamMap& p0pars_ave = flavTag_->getCalibP0Ave();
 		LauTagCatParamMap& p0pars_delta = flavTag_->getCalibP0Delta();
 		LauTagCatParamMap& p1pars_ave = flavTag_->getCalibP1Ave();
 		LauTagCatParamMap& p1pars_delta = flavTag_->getCalibP1Delta();
 
 		LauParameterPList& fitVars = this->fitPars();
 
 		for(LauTagCatParamMap::iterator iter = p0pars_ave.begin(); iter != p0pars_ave.end(); ++iter){
 			LauParameter* p0 = &((*iter).second);
 			fitVars.push_back(p0);
 			++nCalibPar_;
 		}
 		for(LauTagCatParamMap::iterator iter = p0pars_delta.begin(); iter != p0pars_delta.end(); ++iter){
 			LauParameter* p0 = &((*iter).second);
 			fitVars.push_back(p0);
 			++nCalibPar_;
 		}
 		for(LauTagCatParamMap::iterator iter = p1pars_ave.begin(); iter != p1pars_ave.end(); ++iter){
 			LauParameter* p1 = &((*iter).second);
 			fitVars.push_back(p1);
 			++nCalibPar_;
 		}
 		for(LauTagCatParamMap::iterator iter = p1pars_delta.begin(); iter != p1pars_delta.end(); ++iter){
 			LauParameter* p1 = &((*iter).second);
 			fitVars.push_back(p1);
 			++nCalibPar_;
 		}
 	} else {
 		LauTagCatParamMap& p0pars_b0 = flavTag_->getCalibP0B0();
 		LauTagCatParamMap& p0pars_b0bar = flavTag_->getCalibP0B0bar();
 		LauTagCatParamMap& p1pars_b0 = flavTag_->getCalibP1B0();
 		LauTagCatParamMap& p1pars_b0bar = flavTag_->getCalibP1B0bar();
 
 		LauParameterPList& fitVars = this->fitPars();
 
 		for(LauTagCatParamMap::iterator iter = p0pars_b0.begin(); iter != p0pars_b0.end(); ++iter){
 			LauParameter* p0 = &((*iter).second);
 			fitVars.push_back(p0);
 			++nCalibPar_;
 		}
 		for(LauTagCatParamMap::iterator iter = p0pars_b0bar.begin(); iter != p0pars_b0bar.end(); ++iter){
 			LauParameter* p0 = &((*iter).second);
 			fitVars.push_back(p0);
 			++nCalibPar_;
 		}
 		for(LauTagCatParamMap::iterator iter = p1pars_b0.begin(); iter != p1pars_b0.end(); ++iter){
 			LauParameter* p1 = &((*iter).second);
 			fitVars.push_back(p1);
 			++nCalibPar_;
 		}
 		for(LauTagCatParamMap::iterator iter = p1pars_b0bar.begin(); iter != p1pars_b0bar.end(); ++iter){
 			LauParameter* p1 = &((*iter).second);
 			fitVars.push_back(p1);
 			++nCalibPar_;
 		}
 	}
 
 }
 
 void LauTimeDepFitModel::setEffiParams()
 {
 	LauParameterPList& fitVars = this->fitPars();
 	LauTagCatDtPdfMap& signalDecayTimePdfs = this->getSignalDecayTimePdfs();
 	
 	// Loop over DecayTimePdfs
 	for (LauTagCatDtPdfMap::iterator iter = signalDecayTimePdfs.begin(); iter != signalDecayTimePdfs.end(); ++iter){
 		LauDecayTimePdf* pdf = (*iter).second;
 		std::vector<LauParameter*>& effiPars = pdf->getEffiPars(); 
 		for(std::vector<LauParameter*>::iterator iter2 = effiPars.begin(); iter2 != effiPars.end(); ++iter2){
 			LauParameter* par = *iter2;
 			fitVars.push_back(par);
 			++nEffiPar_;
 		}
 	}
 }
 
 void LauTimeDepFitModel::setExtraNtupleVars()
 {
 	// Set-up other parameters derived from the fit results, e.g. fit fractions.
 
 	if (this->useDP() != kTRUE) {
 		return;
 	}
 
 	// First clear the vectors so we start from scratch
 	this->clearExtraVarVectors();
 
 	LauParameterList& extraVars = this->extraPars();
 
 	// Add the B0 and B0bar fit fractions for each signal component
 	fitFracB0bar_ = sigModelB0bar_->getFitFractions();
 	if (fitFracB0bar_.size() != nSigComp_) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0bar_.size()<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	for (UInt_t i(0); i<nSigComp_; ++i) {
 		if (fitFracB0bar_[i].size() != nSigComp_) {
 			std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0bar_[i].size()<<std::endl;
 			gSystem->Exit(EXIT_FAILURE);
 		}
 	}
 
 	for (UInt_t i(0); i<nSigComp_; ++i) {
 		for (UInt_t j = i; j < nSigComp_; j++) {
 			TString name = fitFracB0bar_[i][j].name();
 			name.Insert( name.Index("FitFrac"), "B0bar" );
 			fitFracB0bar_[i][j].name(name);
 			extraVars.push_back(fitFracB0bar_[i][j]);
 		}
 	}
 
 	fitFracB0_ = sigModelB0_->getFitFractions();
 	if (fitFracB0_.size() != nSigComp_) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0_.size()<<std::endl;
 		gSystem->Exit(EXIT_FAILURE);
 	}
 	for (UInt_t i(0); i<nSigComp_; ++i) {
 		if (fitFracB0_[i].size() != nSigComp_) {
 			std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0_[i].size()<<std::endl;
 			gSystem->Exit(EXIT_FAILURE);
 		}
 	}
 
 	for (UInt_t i(0); i<nSigComp_; ++i) {
 		for (UInt_t j = i; j < nSigComp_; j++) {
 			TString name = fitFracB0_[i][j].name();
 			name.Insert( name.Index("FitFrac"), "B0" );
 			fitFracB0_[i][j].name(name);
 			extraVars.push_back(fitFracB0_[i][j]);
 		}
 	}
 
 	// Calculate the ACPs and FitFrac asymmetries
 	this->calcAsymmetries(kTRUE);
 
 	// Add the Fit Fraction asymmetry for each signal component
 	for (UInt_t i = 0; i < nSigComp_; i++) {
 		extraVars.push_back(fitFracAsymm_[i]);
 	}
 
 	// Add the calculated CP asymmetry for each signal component
 	for (UInt_t i = 0; i < nSigComp_; i++) {
 		extraVars.push_back(acp_[i]);
 	}
 
 	// Now add in the DP efficiency values
 	Double_t initMeanEffB0bar = sigModelB0bar_->getMeanEff().initValue();
 	meanEffB0bar_.value(initMeanEffB0bar); meanEffB0bar_.initValue(initMeanEffB0bar); meanEffB0bar_.genValue(initMeanEffB0bar);
 	extraVars.push_back(meanEffB0bar_);
 
 	Double_t initMeanEffB0 = sigModelB0_->getMeanEff().initValue();
 	meanEffB0_.value(initMeanEffB0); meanEffB0_.initValue(initMeanEffB0); meanEffB0_.genValue(initMeanEffB0);
 	extraVars.push_back(meanEffB0_);
 
 	// Also add in the DP rates
 	Double_t initDPRateB0bar = sigModelB0bar_->getDPRate().initValue();
 	DPRateB0bar_.value(initDPRateB0bar); DPRateB0bar_.initValue(initDPRateB0bar); DPRateB0bar_.genValue(initDPRateB0bar);
 	extraVars.push_back(DPRateB0bar_);
 
 	Double_t initDPRateB0 = sigModelB0_->getDPRate().initValue();
 	DPRateB0_.value(initDPRateB0); DPRateB0_.initValue(initDPRateB0); DPRateB0_.genValue(initDPRateB0);
 	extraVars.push_back(DPRateB0_);
 }
 
 void LauTimeDepFitModel::setAsymmetries(const Double_t AProd, const Bool_t AProdFix, const Double_t ARaw, const Bool_t ARawFix){
 	AProd_.value(AProd);
 	AProd_.fixed(AProdFix);
 	ARaw_.value(ARaw);
 	ARaw_.fixed(ARawFix);
 
 	//this->setAsymParams();
 }
 
 void LauTimeDepFitModel::finaliseFitResults(const TString& tablePrefixName)
 {
 	// Retrieve parameters from the fit results for calculations and toy generation
 	// and eventually store these in output root ntuples/text files
 
 	// Now take the fit parameters and update them as necessary
 	// i.e. to make mag > 0.0, phase in the right range.
 	// This function will also calculate any other values, such as the
 	// fit fractions, using any errors provided by fitParErrors as appropriate.
 	// Also obtain the pull values: (measured - generated)/(average error)
 
 	if (this->useDP() == kTRUE) {
 		for (UInt_t i = 0; i < nSigComp_; ++i) {
 			// Check whether we have "a > 0.0", and phases in the right range
 			coeffPars_[i]->finaliseValues();
 		}
 	}
 
 
 	// update the pulls on the event fractions and asymmetries
 	if (this->doEMLFit()) {
 		signalEvents_->updatePull();
 	}
 	if (this->useDP() == kFALSE) {
 		signalAsym_->updatePull();
 	}
 
 	// Finalise the pulls on the decay time parameters
 	for (LauTagCatDtPdfMap::iterator iter = signalDecayTimePdfs_.begin(); iter != signalDecayTimePdfs_.end(); ++iter) {
 		LauDecayTimePdf* pdf = (*iter).second;
 		pdf->updatePulls();
 	}
 
 	if (useSinCos_) {
 		cosPhiMix_.updatePull();
 		sinPhiMix_.updatePull();
 	} else {
 		this->checkMixingPhase();
 	}
 
 	// Update the pulls on all the extra PDFs' parameters
 	for (LauTagCatPdfListMap::iterator iter = sigExtraPdf_.begin(); iter != sigExtraPdf_.end(); ++iter) {
 		this->updateFitParameters(iter->second);
 	}
 
 	LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac();
 
 	// Tagging-category fractions for signal and background events
 	Double_t firstCatFrac(1.0);
 	Int_t firstCat(0);
 	for (LauTagCatParamMap::iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) {
 		if (iter == signalTagCatFrac.begin()) {
 			firstCat = iter->first;
 			continue;
 		}
 		LauParameter& par = (*iter).second;
 		firstCatFrac -= par.value();
 
 		// update the parameter pull
 		par.updatePull();
 	}
 	signalTagCatFrac[firstCat].value(firstCatFrac);
 	signalTagCatFrac[firstCat].updatePull();
 
 	// Fill the fit results to the ntuple
 
 	// update the coefficients and then calculate the fit fractions and ACP's
 	if (this->useDP() == kTRUE) {
 		this->updateCoeffs();
 		sigModelB0bar_->updateCoeffs(coeffsB0bar_);  sigModelB0bar_->calcExtraInfo();
 		sigModelB0_->updateCoeffs(coeffsB0_);        sigModelB0_->calcExtraInfo();
 
 		LauParArray fitFracB0bar = sigModelB0bar_->getFitFractions();
 		if (fitFracB0bar.size() != nSigComp_) {
 			std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0bar.size()<<std::endl;
 			gSystem->Exit(EXIT_FAILURE);
 		}
 		for (UInt_t i(0); i<nSigComp_; ++i) {
 			if (fitFracB0bar[i].size() != nSigComp_) {
 				std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0bar[i].size()<<std::endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 		}
 		LauParArray fitFracB0 = sigModelB0_->getFitFractions();
 		if (fitFracB0.size() != nSigComp_) {
 			std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0.size()<<std::endl;
 			gSystem->Exit(EXIT_FAILURE);
 		}
 		for (UInt_t i(0); i<nSigComp_; ++i) {
 			if (fitFracB0[i].size() != nSigComp_) {
 				std::cerr<<"ERROR in LauTimeDepFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: "<<fitFracB0[i].size()<<std::endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 		}
 
 		for (UInt_t i(0); i<nSigComp_; ++i) {
 			for (UInt_t j(i); j<nSigComp_; ++j) {
 				fitFracB0bar_[i][j].value(fitFracB0bar[i][j].value());
 				fitFracB0_[i][j].value(fitFracB0[i][j].value());
 			}
 		}
 
 		meanEffB0bar_.value(sigModelB0bar_->getMeanEff().value());
 		meanEffB0_.value(sigModelB0_->getMeanEff().value());
 		DPRateB0bar_.value(sigModelB0bar_->getDPRate().value());
 		DPRateB0_.value(sigModelB0_->getDPRate().value());
 
 		this->calcAsymmetries();
 
 		// Then store the final fit parameters, and any extra parameters for
 		// the signal model (e.g. fit fractions, FF asymmetries, ACPs, mean efficiency and DP rate)
 
 		this->clearExtraVarVectors();
 		LauParameterList& extraVars = this->extraPars();
 
 		for (UInt_t i(0); i<nSigComp_; ++i) {
 			for (UInt_t j(i); j<nSigComp_; ++j) {
 				extraVars.push_back(fitFracB0bar_[i][j]);
 			}
 		}
 		for (UInt_t i(0); i<nSigComp_; ++i) {
 			for (UInt_t j(i); j<nSigComp_; ++j) {
 				extraVars.push_back(fitFracB0_[i][j]);
 			}
 		}
 		for (UInt_t i(0); i<nSigComp_; ++i) {
 			extraVars.push_back(fitFracAsymm_[i]);
 		}
 		for (UInt_t i(0); i<nSigComp_; ++i) {
 			extraVars.push_back(acp_[i]);
 		}
 		extraVars.push_back(meanEffB0bar_);
 		extraVars.push_back(meanEffB0_);
 		extraVars.push_back(DPRateB0bar_);
 		extraVars.push_back(DPRateB0_);
 
 		this->printFitFractions(std::cout);
 		this->printAsymmetries(std::cout);
 	}
 
 	const LauParameterPList& fitVars = this->fitPars();
 	const LauParameterList& extraVars = this->extraPars();
 	LauFitNtuple* ntuple = this->fitNtuple();
 	ntuple->storeParsAndErrors(fitVars, extraVars);
 
 	// find out the correlation matrix for the parameters
 	ntuple->storeCorrMatrix(this->iExpt(), this->fitStatus(), this->covarianceMatrix());
 
 	// Fill the data into ntuple
 	ntuple->updateFitNtuple();
 
 	// Print out the partial fit fractions, phases and the
 	// averaged efficiency, reweighted by the dynamics (and anything else)
 	if (this->writeLatexTable()) {
 		TString sigOutFileName(tablePrefixName);
 		sigOutFileName += "_"; sigOutFileName += this->iExpt(); sigOutFileName += "Expt.tex";
 		this->writeOutTable(sigOutFileName);
 	}
 }
 
 void LauTimeDepFitModel::printFitFractions(std::ostream& output)
 {
 	// Print out Fit Fractions, total DP rate and mean efficiency
 	// First for the B0bar events
 	for (UInt_t i = 0; i < nSigComp_; i++) {
 		const TString compName(coeffPars_[i]->name());
 		output<<"B0bar FitFraction for component "<<i<<" ("<<compName<<") = "<<fitFracB0bar_[i][i]<<std::endl;
 	}
 	output<<"B0bar overall DP rate (integral of matrix element squared) = "<<DPRateB0bar_<<std::endl;
 	output<<"B0bar average efficiency weighted by whole DP dynamics = "<<meanEffB0bar_<<std::endl;
 
 	// Then for the B0 sample
 	for (UInt_t i = 0; i < nSigComp_; i++) {
 		const TString compName(coeffPars_[i]->name());
 		const TString conjName(sigModelB0bar_->getConjResName(compName));
 		output<<"B0 FitFraction for component "<<i<<" ("<<conjName<<") = "<<fitFracB0_[i][i]<<std::endl;
 	}
 	output<<"B0 overall DP rate (integral of matrix element squared) = "<<DPRateB0_<<std::endl;
 	output<<"B0 average efficiency weighted by whole DP dynamics = "<<meanEffB0_<<std::endl;
 }
 
 void LauTimeDepFitModel::printAsymmetries(std::ostream& output)
 {
 	for (UInt_t i = 0; i < nSigComp_; i++) {
 		const TString compName(coeffPars_[i]->name());
 		output<<"Fit Fraction asymmetry for component "<<i<<" ("<<compName<<") = "<<fitFracAsymm_[i]<<std::endl;
 	}
 	for (UInt_t i = 0; i < nSigComp_; i++) {
 		const TString compName(coeffPars_[i]->name());
 		output<<"ACP for component "<<i<<" ("<<compName<<") = "<<acp_[i].value()<<" +- "<<acp_[i].error()<<std::endl;
 	}
 }
 
 void LauTimeDepFitModel::writeOutTable(const TString& outputFile)
 {
 	// Write out the results of the fit to a tex-readable table
 	std::ofstream fout(outputFile);
 	LauPrint print;
 
 	std::cout<<"INFO in LauTimeDepFitModel::writeOutTable : Writing out results of the fit to the tex file "<<outputFile<<std::endl;
 
 	if (this->useDP() == kTRUE) {
 		// print the fit coefficients in one table
 		coeffPars_.front()->printTableHeading(fout);
 		for (UInt_t i = 0; i < nSigComp_; i++) {
 			coeffPars_[i]->printTableRow(fout);
 		}
 		fout<<"\\hline"<<std::endl;
 		fout<<"\\end{tabular}"<<std::endl<<std::endl;
 
 		// print the fit fractions in another
 		fout<<"\\begin{tabular}{|l|c|c|c|c|}"<<std::endl;
 		fout<<"\\hline"<<std::endl;
 		fout<<"Component & \\Bzb\\ Fit Fraction & \\Bz\\ Fit Fraction & Fit Fraction Asymmetry & $A_{\\CP}$ \\\\"<<std::endl;
 		fout<<"\\hline"<<std::endl;
 		Double_t fitFracSumB0bar(0.0);
 		Double_t fitFracSumB0(0.0);
 		for (UInt_t i = 0; i < nSigComp_; i++) {
 			TString resName = coeffPars_[i]->name();
 			resName = resName.ReplaceAll("_", "\\_");
 			fout<<resName<<"  &  $";
 			Double_t fitFracB0bar = fitFracB0bar_[i][i].value();
 			fitFracSumB0bar += fitFracB0bar;
 			print.printFormat(fout, fitFracB0bar);
 			fout << "$  &  $" << std::endl;
 			Double_t fitFracB0 = fitFracB0_[i][i].value();
 			fitFracSumB0 += fitFracB0;
 			print.printFormat(fout, fitFracB0);
 			fout << "$  &  $" << std::endl;
 			Double_t fitFracAsymm = fitFracAsymm_[i].value();
 			print.printFormat(fout, fitFracAsymm);
 			fout << "$  &  $" << std::endl;
 			Double_t acp = acp_[i].value();
 			Double_t acpErr = acp_[i].error();
 			print.printFormat(fout, acp);
 			fout<<" \\pm ";
 			print.printFormat(fout, acpErr);
 			fout<<"$ \\\\"<<std::endl;
 		}
 		fout<<"\\hline"<<std::endl;
 
 
 		// Also print out sum of fit fractions
 		fout << "Fit Fraction Sum = & $";
 		print.printFormat(fout, fitFracSumB0bar);
 		fout << "$  &  $";
 		print.printFormat(fout, fitFracSumB0);
 		fout << "$  &   &   \\\\" << std::endl;
 
 		fout << "\\hline \n\\hline" << std::endl;
 
 		fout << "DP rate = & $";
 		print.printFormat(fout, DPRateB0bar_.value());
 		fout << "$  &  $";
 		print.printFormat(fout, DPRateB0_.value());
 		fout << "$  &   &   \\\\" << std::endl;
 
 		fout << "$< \\varepsilon > =$ & $";
 		print.printFormat(fout, meanEffB0bar_.value());
 		fout << "$  &  $";
 		print.printFormat(fout, meanEffB0_.value());
 		fout << "$  &   &   \\\\" << std::endl;
 
 		if (useSinCos_) {
 			fout << "$\\sinPhiMix =$ & $";
 			print.printFormat(fout, sinPhiMix_.value());
 			fout << " \\pm ";
 			print.printFormat(fout, sinPhiMix_.error());
 			fout << "$ &  &  &  &  &  &  &  \\\\" << std::endl;
 
 			fout << "$\\cosPhiMix =$ & $";
 			print.printFormat(fout, cosPhiMix_.value());
 			fout << " \\pm ";
 			print.printFormat(fout, cosPhiMix_.error());
 			fout << "$ &  &  &  &  &  &  &  \\\\" << std::endl;
 		} else {
 			fout << "$\\phiMix =$ & $";
 			print.printFormat(fout, phiMix_.value());
 			fout << " \\pm ";
 			print.printFormat(fout, phiMix_.error());
 			fout << "$ &  &  &  &  &  &  &  \\\\" << std::endl;
 		}
 
 		fout << "\\hline \n\\end{tabular}" << std::endl;
 	}
 
 	if (!sigExtraPdf_.empty()) {
 		fout<<"\\begin{tabular}{|l|c|}"<<std::endl;
 		fout<<"\\hline"<<std::endl;
 		fout<<"\\Extra Signal PDFs' Parameters: & \\\\"<<std::endl;
 		for (LauTagCatPdfListMap::const_iterator iter = sigExtraPdf_.begin(); iter != sigExtraPdf_.end(); ++iter) {
 			this->printFitParameters(iter->second, fout);
 		}
 		fout<<"\\hline \n\\end{tabular}"<<std::endl<<std::endl;
 	}
 }
 
 void LauTimeDepFitModel::checkInitFitParams()
 {
 	// Update the number of signal events to be total-sum(background events)
 	this->updateSigEvents();
 
 	// Check whether we want to have randomised initial fit parameters for the signal model
 	if (this->useRandomInitFitPars() == kTRUE) {
 		this->randomiseInitFitPars();
 	}
 }
 
 void LauTimeDepFitModel::randomiseInitFitPars()
 {
 	// Only randomise those parameters that are not fixed!
 	std::cout<<"INFO in LauTimeDepFitModel::randomiseInitFitPars : Randomising the initial values of the coefficients of the DP components (and phiMix)..."<<std::endl;
 
 	for (UInt_t i = 0; i < nSigComp_; i++) {
 		coeffPars_[i]->randomiseInitValues();
 	}
 
 	phiMix_.randomiseValue(-LauConstants::pi, LauConstants::pi);
 	if (useSinCos_) {
 		sinPhiMix_.initValue(TMath::Sin(phiMix_.initValue()));
 		cosPhiMix_.initValue(TMath::Cos(phiMix_.initValue()));
 	}
 }
 
 LauTimeDepFitModel::LauGenInfo LauTimeDepFitModel::eventsToGenerate()
 {
 	// Determine the number of events to generate for each hypothesis
 	// If we're smearing then smear each one individually
 	// NB this individual smearing has to be done individually per tagging category as well
 
 	LauGenInfo nEvtsGen;
 
 	LauTagCatGenInfo eventsB0, eventsB0bar;
 
 	// Signal
 	// If we're including the DP and decay time we can't decide on the tag
 	// yet, since it depends on the whole DP+dt PDF, however, if
 	// we're not then we need to decide.
 	Double_t evtWeight(1.0);
 	Double_t nEvts = signalEvents_->genValue();
 	if ( nEvts < 0.0 ) {
 		evtWeight = -1.0;
 		nEvts = TMath::Abs( nEvts );
 	}
 
         LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac();
 
 	Double_t sigAsym(0.0);
 	if (this->useDP() == kFALSE) {
 		sigAsym = signalAsym_->genValue();
 
 		for (LauTagCatParamMap::const_iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) {
 			const LauParameter& par = iter->second;
 			Double_t eventsbyTagCat = par.value() * nEvts;
 			Double_t eventsB0byTagCat = TMath::Nint(eventsbyTagCat/2.0 * (1.0 - sigAsym));
 			Double_t eventsB0barbyTagCat  = TMath::Nint(eventsbyTagCat/2.0 * (1.0 + sigAsym));
 			if (this->doPoissonSmearing()) {
 				eventsB0byTagCat = LauRandom::randomFun()->Poisson(eventsB0byTagCat);
 				eventsB0barbyTagCat = LauRandom::randomFun()->Poisson(eventsB0barbyTagCat);
 			}
 			eventsB0[iter->first] = std::make_pair( TMath::Nint(eventsB0byTagCat), evtWeight );
 			eventsB0bar[iter->first] = std::make_pair( TMath::Nint(eventsB0barbyTagCat), evtWeight );
 		}
         // CONVENTION WARNING
 		nEvtsGen[std::make_pair("signal",-1)] = eventsB0;
 		nEvtsGen[std::make_pair("signal",+1)] = eventsB0bar;
 	} else {
 		Double_t rateB0bar = sigModelB0bar_->getDPRate().value();
 		Double_t rateB0 = sigModelB0_->getDPRate().value();
 		if ( rateB0bar+rateB0 > 1e-30) {
 			sigAsym = (rateB0bar-rateB0)/(rateB0bar+rateB0);
 		}
 
 		for (LauTagCatParamMap::const_iterator iter = signalTagCatFrac.begin(); iter != signalTagCatFrac.end(); ++iter) {
 			const LauParameter& par = iter->second;
 			Double_t eventsbyTagCat = par.value() * nEvts;
 			if (this->doPoissonSmearing()) {
 				eventsbyTagCat = LauRandom::randomFun()->Poisson(eventsbyTagCat);
 			}
 			eventsB0[iter->first] = std::make_pair( TMath::Nint(eventsbyTagCat), evtWeight );
 		}
 		nEvtsGen[std::make_pair("signal",0)] = eventsB0; // generate signal event, decide tag later.
 	}
 
 	std::cout<<"INFO in LauTimeDepFitModel::eventsToGenerate : Generating toy MC with:"<<std::endl;
 	std::cout<<"                                             : Signal asymmetry  = "<<sigAsym<<" and number of signal events = "<<signalEvents_->genValue()<<std::endl;
 
 	return nEvtsGen;
 }
 
 Bool_t LauTimeDepFitModel::genExpt()
 {
 	// Routine to generate toy Monte Carlo events according to the various models we have defined.
 
 	// Determine the number of events to generate for each hypothesis
 	LauGenInfo nEvts = this->eventsToGenerate();
 
 	Bool_t genOK(kTRUE);
 	Int_t evtNum(0);
 
 	// Loop over the hypotheses and generate the appropriate number of
 	// events for each one
 	for (LauGenInfo::const_iterator iter = nEvts.begin(); iter != nEvts.end(); ++iter) {
 
 		// find the category of events (e.g. signal)
 		const TString& evtCategory(iter->first.first);
 
 		// find the true flavour
 		curEvtTagFlv_ = iter->first.second;
 
 		// loop through each tagging category
 		const LauTagCatGenInfo& eventsByTagCat = iter->second;
 		for (LauTagCatGenInfo::const_iterator tagCatIter = eventsByTagCat.begin(); tagCatIter != eventsByTagCat.end(); ++tagCatIter) {
 
 			// get the tagging category
 			curEvtTagCat_ = tagCatIter->first;
 
 			// get the event weight for this category
 			const Double_t evtWeight( tagCatIter->second.second );
 
 			// get the number of events to generate in this configuration
 			const Int_t nEvtsGen( tagCatIter->second.first );
 
 			for (Int_t iEvt(0); iEvt<nEvtsGen; ++iEvt) {
 
 				this->setGenNtupleDoubleBranchValue( "evtWeight", evtWeight );
 
 				if (evtCategory == "signal") {
 					this->setGenNtupleIntegerBranchValue("genSig",1);
 					// All the generate*Event() methods have to fill in curEvtDecayTime_ and curEvtDecayTimeErr_
 					// In addition, generateSignalEvent has to decide on the tag and fill in curEvtTagFlv_
 					genOK = this->generateSignalEvent();
 				} else {
 					genOK = kFALSE;
 				}
 
 				if (!genOK) {
 					// If there was a problem with the generation then break out and return.
 					// The problem model will have adjusted itself so that all should be OK next time.
 					break;
 				}
 
 				if (this->useDP() == kTRUE) {
 					this->setDPDtBranchValues();	// store DP, decay time and tagging variables in the ntuple
 				}
 
 				// If the tagging category is 0, we don't know what the flavour is...
 				if (curEvtTagCat_==0){
 					curEvtTagFlv_=0;
 				}
 
 				// Store the event's tag and tagging category
 				this->setGenNtupleIntegerBranchValue("cpEigenvalue", cpEigenValue_);
 				this->setGenNtupleIntegerBranchValue("tagCat",curEvtTagCat_);
 				this->setGenNtupleIntegerBranchValue("tagFlv",curEvtTagFlv_);
 				this->setGenNtupleDoubleBranchValue(flavTag_->getMistagVarName(),curEvtMistag_);
 				this->setGenNtupleDoubleBranchValue(flavTag_->getTrueTagVarName(),curEvtTrueTagFlv_);
 
 				//Store calibration parameters and asymmetries during debugging at least
 				//this->setGenNtupleDoubleBranchValue("AProd", AProd_->unblindValue());
 				//this->setGenNtupleDoubleBranchValue("ARaw", ARaw_->unblindValue());
 
 				// Store the event number (within this experiment)
 				// and then increment it
 				this->setGenNtupleIntegerBranchValue("iEvtWithinExpt",evtNum);
 				++evtNum;
 
 				// Write the values into the tree
 				this->fillGenNtupleBranches();
 
 				// Print an occasional progress message
 				if (iEvt%1000 == 0) {std::cout<<"INFO in LauTimeDepFitModel::genExpt : Generated event number "<<iEvt<<" out of "<<nEvtsGen<<" "<<evtCategory<<" events."<<std::endl;}
 
 			}	// end of loop over events for given species, tagCat and tagFlv
 
 			if (!genOK) {
 				break;
 			}
 
 		}	// end of loop over tagCats for the given species and tagFlv
 		if (!genOK) {
 			break;
 		}
 
 	}	//end of loop over species and tagFlv.
 
 	if (this->useDP() && genOK) {
 
 		sigModelB0bar_->checkToyMC(kTRUE);
 		sigModelB0_->checkToyMC(kTRUE);
 
 		std::cout<<"aSqMaxSet = "<<aSqMaxSet_<<" and aSqMaxVar = "<<aSqMaxVar_<<std::endl;
 
 		// Get the fit fractions if they're to be written into the latex table
 		if (this->writeLatexTable()) {
 			LauParArray fitFracB0bar = sigModelB0bar_->getFitFractions();
 			if (fitFracB0bar.size() != nSigComp_) {
 				std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0bar.size()<<std::endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 			for (UInt_t i(0); i<nSigComp_; ++i) {
 				if (fitFracB0bar[i].size() != nSigComp_) {
 					std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0bar[i].size()<<std::endl;
 					gSystem->Exit(EXIT_FAILURE);
 				}
 			}
 			LauParArray fitFracB0 = sigModelB0_->getFitFractions();
 			if (fitFracB0.size() != nSigComp_) {
 				std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0.size()<<std::endl;
 				gSystem->Exit(EXIT_FAILURE);
 			}
 			for (UInt_t i(0); i<nSigComp_; ++i) {
 				if (fitFracB0[i].size() != nSigComp_) {
 					std::cerr<<"ERROR in LauTimeDepFitModel::generate : Fit Fraction array of unexpected dimension: "<<fitFracB0[i].size()<<std::endl;
 					gSystem->Exit(EXIT_FAILURE);
 				}
 			}
 			for (UInt_t i(0); i<nSigComp_; ++i) {
 				for (UInt_t j(i); j<nSigComp_; ++j) {
 					fitFracB0bar_[i][j].value(fitFracB0bar[i][j].value());
 					fitFracB0_[i][j].value(fitFracB0[i][j].value());
 				}
 			}
 			meanEffB0bar_.value(sigModelB0bar_->getMeanEff().value());
 			meanEffB0_.value(sigModelB0_->getMeanEff().value());
 			DPRateB0bar_.value(sigModelB0bar_->getDPRate().value());
 			DPRateB0_.value(sigModelB0_->getDPRate().value());
 		}
 	}
 
 	// If we're reusing embedded events or if the generation is being
 	// reset then clear the lists of used events
 	//if (!signalTree_.empty() && (reuseSignal_ || !genOK)) {
 	if (reuseSignal_ || !genOK) {
 		for(LauTagCatEmbDataMap::const_iterator iter = signalTree_.begin(); iter != signalTree_.end(); ++iter) {
 			(iter->second)->clearUsedList();
 		}
 	}
 
 	return genOK;
 }
 
 Bool_t LauTimeDepFitModel::generateSignalEvent()
 {
 	// Generate signal event, including SCF if necessary.
 	// DP:DecayTime generation follows.
 	// If it's ok, we then generate mES, DeltaE, Fisher/NN...
 	Bool_t genOK(kTRUE);
 	Bool_t generatedEvent(kFALSE);
 
 	Bool_t doSquareDP = kinematicsB0bar_->squareDP();
 	doSquareDP &= kinematicsB0_->squareDP();
 
 	LauKinematics* kinematics(kinematicsB0bar_);
 
 	// find the right decay time PDF for the current tagging category
 	LauTagCatDtPdfMap::const_iterator dt_iter = signalDecayTimePdfs_.find(curEvtTagCat_);
 	LauDecayTimePdf* decayTimePdf = (dt_iter != signalDecayTimePdfs_.end()) ? dt_iter->second : 0;
 
 	// add flavour tagging for time pdf
 	decayTimePdf->setFlavourTagging(flavTag_);
 
 	// find the right embedded data for the current tagging category
 	LauTagCatEmbDataMap::const_iterator emb_iter = signalTree_.find(curEvtTagCat_);
 	LauEmbeddedData* embeddedData = (emb_iter != signalTree_.end()) ? emb_iter->second : 0;
 
 	// find the right extra PDFs for the current tagging category
 	LauTagCatPdfListMap::iterator extra_iter = sigExtraPdf_.find(curEvtTagCat_);
 	LauPdfList* extraPdfs = (extra_iter != sigExtraPdf_.end()) ? &(extra_iter->second) : 0;
 
 	if (this->useDP()) {
 		if (embeddedData) {
 			embeddedData->getEmbeddedEvent(kinematics);
 			curEvtTagFlv_ = TMath::Nint(embeddedData->getValue("tagFlv"));
 			curEvtDecayTimeErr_ = embeddedData->getValue(decayTimePdf->varErrName());
 			curEvtDecayTime_ = embeddedData->getValue(decayTimePdf->varName());
 			if (embeddedData->haveBranch("mcMatch")) {
 				Int_t match = TMath::Nint(embeddedData->getValue("mcMatch"));
 				if (match) {
 					this->setGenNtupleIntegerBranchValue("genTMSig",1);
 					this->setGenNtupleIntegerBranchValue("genSCFSig",0);
 				} else {
 					this->setGenNtupleIntegerBranchValue("genTMSig",0);
 					this->setGenNtupleIntegerBranchValue("genSCFSig",1);
 				}
 			}
 		} else {
 
 			nGenLoop_ = 0;
 			// generate the decay time error (NB the kTRUE forces the generation of a new value)
 			curEvtDecayTimeErr_ = decayTimePdf->generateError(kTRUE);
 
 			while (generatedEvent == kFALSE && nGenLoop_ < iterationsMax_) {
 
 				// Calculate the unnormalised truth-matched signal likelihood
 
 				// First let define the tag flavour CONVENTION WARNING
 				Double_t randNo = LauRandom::randomFun()->Rndm();
 				if (randNo < 0.5) {
 					curEvtTagFlv_ = +1; // B0bar tag
 				} else {
 					curEvtTagFlv_ = -1; // B0 tag
 				}
 
 				curEvtTrueTagFlv_ = 0;
 
 				if (decayTimePdf->getFuncType() == LauDecayTimePdf::SimFitNormBd || decayTimePdf->getFuncType() == LauDecayTimePdf::SimFitNormBs || curEvtTagCat_==0){
 					// First let define the tag flavour CONVENTION WARNING
 					Double_t rand = LauRandom::randomFun()->Rndm();
 					if (rand < 0.5) {
 						curEvtTrueTagFlv_ = +1; // B0bar tag
 					} else {
 						curEvtTrueTagFlv_ = -1; // B0 tag
 					}
 				}
 
 				// Calculate event quantities that depend only on the tagCat and tagFlv
 				Double_t qD(0.);
 				Double_t qDDo2(0.);
 	
 				if(flavTag_->getUsePerEvtMistag() && curEvtTagCat_!=0 ){
 					curEvtMistag_ = flavTag_->getEtaGen(sigFlavTagPdf_[curEvtTagCat_]);
 					Double_t omega = flavTag_->getOmegaGen(curEvtMistag_,curEvtTagFlv_,curEvtTagCat_);
 					qD = curEvtTagFlv_*(1.0-2.0*omega);
 					//qD = curEvtTagFlv_*(1-2*curEvtMistag_);
 				}else{
 					curEvtMistag_ = 0.5;
                                         LauTagCatParamMap dilution_ = flavTag_->getDilution();
                                         LauTagCatParamMap deltaDilution_ = flavTag_->getDeltaDilution();
 
 					qD = curEvtTagFlv_*dilution_[curEvtTagCat_].unblindValue();
 					qDDo2 = curEvtTagFlv_*0.5*deltaDilution_[curEvtTagCat_].unblindValue();
 				}
 
 				// Generate the DP position
 				Double_t m13Sq(0.0), m23Sq(0.0);
 				kinematicsB0bar_->genFlatPhaseSpace(m13Sq, m23Sq);
 
 				// Next, calculate the total A and Abar for the given DP position
 				sigModelB0_->calcLikelihoodInfo(m13Sq, m23Sq);
 				sigModelB0bar_->calcLikelihoodInfo(m13Sq, m23Sq);
 
 				// Retrieve the amplitudes and efficiency from the dynamics
 				const LauComplex& Abar = sigModelB0bar_->getEvtDPAmp();
 				const LauComplex& A = sigModelB0_->getEvtDPAmp();
 				Double_t eff = sigModelB0bar_->getEvtEff();
 
 				// Next calculate the DP terms
 				Double_t aSqSum = A.abs2() + Abar.abs2();
 				Double_t aSqDif = A.abs2() - Abar.abs2();
                                 
 				LauComplex inter = Abar * A.conj() * phiMixComplex_;
 				Double_t interTermIm = 2.0 * inter.im();
 				Double_t interTermRe = 2.0 * inter.re();
 
 				// Generate decay time
 				const Double_t tMin = decayTimePdf->minAbscissa();
 				const Double_t tMax = decayTimePdf->maxAbscissa();
 				curEvtDecayTime_ = LauRandom::randomFun()->Rndm()*(tMax-tMin) + tMin;
 
 				// Calculate all the decay time info
 				//decayTimePdf->calcLikelihoodInfo(curEvtDecayTime_,curEvtDecayTimeErr_);
 				decayTimePdf->calcLikelihoodInfoGen(curEvtDecayTime_,curEvtDecayTimeErr_,curEvtTagFlv_,curEvtMistag_,curEvtTagCat_,curEvtTrueTagFlv_);
 				// ...and check that the calculation went ok, otherwise loop again
 				if (decayTimePdf->state() != LauDecayTimePdf::Good) {
 					std::cout<<"decayTimePdf state is bad"<<std::endl;
 					++nGenLoop_;
 					continue;
 				}
 
 				// Get the decay time acceptance
 				Double_t dtEff = decayTimePdf->getEffiTerm();
 
 				// First get all the decay time terms
 				Double_t dtCos = decayTimePdf->getCosTerm();
 				Double_t dtSin = decayTimePdf->getSinTerm();
 				Double_t dtCosh = decayTimePdf->getCoshTerm();
 				Double_t dtSinh = decayTimePdf->getSinhTerm();
 
 				// Combine all terms
 				Double_t cosTerm  = dtCos * qD * aSqDif;
 				Double_t sinTerm  = dtSin * qD * interTermIm;
 				Double_t coshTerm = dtCosh * (1.0 + qDDo2) * aSqSum;
 				Double_t sinhTerm = dtSinh * (1.0 + qDDo2) * interTermRe;
 
 				if ( cpEigenValue_ == CPOdd ) {
 					sinTerm  *= -1.0;
 					sinhTerm *= -1.0;
 				}
 				// ... to get the total and multiply by the efficiency
 				Double_t ASq = coshTerm + cosTerm - sinTerm + sinhTerm;
 				//if (curEvtTrueTagFlv_ !=0 && curEvtTagCat_>0 && curEvtTagFlv_==-1){
 				//	std::cout<<"##### Event details #####"<<std::endl;
                                 //	std::cout<< "dtCos dtSin dtCosh dtSinh " << dtCos << " " << dtSin << " " << dtCosh << " " << dtSinh << std::endl;
                                 //	std::cout<< "qD aSqDif dtSin interTermIm "  << qD << " " << aSqDif << " " << dtSin << " " << interTermIm << std::endl;
 				//	std::cout<<"Cosh Cos Sin Sinh "<< coshTerm<< " " << cosTerm << " " << sinTerm << " " << sinhTerm << std::endl;
                                 //	std::cout<< "Total Amplitude : " << ASq << std::endl;
 				//}
 				//ASq /= decayTimePdf->getNormTerm();
 				ASq *= eff;
 				ASq *= dtEff;
                                 //std::cout << "Total Amplitude Eff: " << ASq << std::endl;
 
 				//Finally we throw the dice to see whether this event should be generated
 				//We make a distinction between the likelihood of TM and SCF to tag the SCF events as such
 				randNo = LauRandom::randomFun()->Rndm();
 				if (randNo <= ASq/aSqMaxSet_ ) {
 					generatedEvent = kTRUE;
 					nGenLoop_ = 0;
 					if (ASq > aSqMaxVar_) {aSqMaxVar_ = ASq;}
 				} else {
 					nGenLoop_++;
 				}
 
 			} // end of while !generatedEvent loop
 
 		} // end of if (embeddedData) else control
 	} else {
 		if ( embeddedData ) {
 			embeddedData->getEmbeddedEvent(0);
 			curEvtTagFlv_ = TMath::Nint(embeddedData->getValue("tagFlv"));
 			curEvtDecayTimeErr_ = embeddedData->getValue(decayTimePdf->varErrName());
 			curEvtDecayTime_ = embeddedData->getValue(decayTimePdf->varName());
 		}
 	}
 
 	// Check whether we have generated the toy MC OK.
 	if (nGenLoop_ >= iterationsMax_) {
 		aSqMaxSet_ = 1.01 * aSqMaxVar_;
 		genOK = kFALSE;
 		std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Hit max iterations: setting aSqMaxSet_ to "<<aSqMaxSet_<<std::endl;
 	} else if (aSqMaxVar_ > aSqMaxSet_) {
 		aSqMaxSet_ = 1.01 * aSqMaxVar_;
 		genOK = kFALSE;
 		std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Found a larger ASq value: setting aSqMaxSet_ to "<<aSqMaxSet_<<std::endl;
 	}
 
 	if (genOK) {
 		//Some variables, like Fisher or NN, might use m13Sq and m23Sq from the kinematics
 		//kinematicsB0bar_ is up to date, update kinematicsB0_
 		kinematicsB0_->updateKinematics(kinematicsB0bar_->getm13Sq(), kinematicsB0bar_->getm23Sq() );
 		this->generateExtraPdfValues(extraPdfs, embeddedData);
 	}
 	// Check for problems with the embedding
 	if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) {
 		std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events."<<std::endl;
 		embeddedData->clearUsedList();
 	}
 
 	return genOK;
 }
 
 void LauTimeDepFitModel::setupGenNtupleBranches()
 {
 	// Setup the required ntuple branches
 	this->addGenNtupleDoubleBranch("evtWeight");
 	this->addGenNtupleIntegerBranch("genSig");
 
 	this->addGenNtupleIntegerBranch("cpEigenvalue");
 	this->addGenNtupleIntegerBranch("tagFlv");
 	this->addGenNtupleIntegerBranch("tagCat");
 
 	if (this->useDP() == kTRUE) {
 		// Let's add the decay time variables.
 		if (signalDecayTimePdfs_.begin() != signalDecayTimePdfs_.end()) {
 			LauDecayTimePdf* pdf = signalDecayTimePdfs_.begin()->second;
 			this->addGenNtupleDoubleBranch(pdf->varName());
 			this->addGenNtupleDoubleBranch(pdf->varErrName());
 		}
 		this->addGenNtupleDoubleBranch("m12");
 		this->addGenNtupleDoubleBranch("m23");
 		this->addGenNtupleDoubleBranch("m13");
 		this->addGenNtupleDoubleBranch("m12Sq");
 		this->addGenNtupleDoubleBranch("m23Sq");
 		this->addGenNtupleDoubleBranch("m13Sq");
 		this->addGenNtupleDoubleBranch("cosHel12");
 		this->addGenNtupleDoubleBranch("cosHel23");
 		this->addGenNtupleDoubleBranch("cosHel13");
 		if (kinematicsB0bar_->squareDP() && kinematicsB0_->squareDP()) {
 			this->addGenNtupleDoubleBranch("mPrime");
 			this->addGenNtupleDoubleBranch("thPrime");
 		}
 
 		// Can add the real and imaginary parts of the B0 and B0bar total
 		// amplitudes seen in the generation (restrict this with a flag
 		// that defaults to false)
 		if ( storeGenAmpInfo_ ) {
 			this->addGenNtupleDoubleBranch("reB0Amp");
 			this->addGenNtupleDoubleBranch("imB0Amp");
 			this->addGenNtupleDoubleBranch("reB0barAmp");
 			this->addGenNtupleDoubleBranch("imB0barAmp");
 		}
 	}
 
 	// Let's look at the extra variables for signal in one of the tagging categories
 	if ( ! sigExtraPdf_.empty() ) {
 		LauPdfList oneTagCatPdfList = sigExtraPdf_.begin()->second;
 		for (LauPdfList::const_iterator pdf_iter = oneTagCatPdfList.begin(); pdf_iter != oneTagCatPdfList.end(); ++pdf_iter) {
 			for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
 				if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
 					this->addGenNtupleDoubleBranch( (*var_iter) );
 				}
 			}
 		}
 	}
 }
 
 void LauTimeDepFitModel::setDPDtBranchValues()
 {
 	// Store the decay time variables.
 	if (signalDecayTimePdfs_.begin() != signalDecayTimePdfs_.end()) {
 		LauDecayTimePdf* pdf = signalDecayTimePdfs_.begin()->second;
 		this->setGenNtupleDoubleBranchValue(pdf->varName(),curEvtDecayTime_);
 		this->setGenNtupleDoubleBranchValue(pdf->varErrName(),curEvtDecayTimeErr_);
 	}
 
     // CONVENTION WARNING
 	LauKinematics* kinematics(0);
 	if (curEvtTagFlv_<0) {
 		kinematics = kinematicsB0_;
 	} else {
 		kinematics = kinematicsB0bar_;
 	}
 
 	// Store all the DP information
 	this->setGenNtupleDoubleBranchValue("m12", kinematics->getm12());
 	this->setGenNtupleDoubleBranchValue("m23", kinematics->getm23());
 	this->setGenNtupleDoubleBranchValue("m13", kinematics->getm13());
 	this->setGenNtupleDoubleBranchValue("m12Sq", kinematics->getm12Sq());
 	this->setGenNtupleDoubleBranchValue("m23Sq", kinematics->getm23Sq());
 	this->setGenNtupleDoubleBranchValue("m13Sq", kinematics->getm13Sq());
 	this->setGenNtupleDoubleBranchValue("cosHel12", kinematics->getc12());
 	this->setGenNtupleDoubleBranchValue("cosHel23", kinematics->getc23());
 	this->setGenNtupleDoubleBranchValue("cosHel13", kinematics->getc13());
 	if (kinematics->squareDP()) {
 		this->setGenNtupleDoubleBranchValue("mPrime", kinematics->getmPrime());
 		this->setGenNtupleDoubleBranchValue("thPrime", kinematics->getThetaPrime());
 	}
 
 	// Can add the real and imaginary parts of the B0 and B0bar total
 	// amplitudes seen in the generation (restrict this with a flag
 	// that defaults to false)
 	if ( storeGenAmpInfo_ ) {
 		if ( this->getGenNtupleIntegerBranchValue("genSig")==1 ) {
 			LauComplex Abar = sigModelB0bar_->getEvtDPAmp();
 			LauComplex A = sigModelB0_->getEvtDPAmp();
 
 			this->setGenNtupleDoubleBranchValue("reB0Amp", A.re());
 			this->setGenNtupleDoubleBranchValue("imB0Amp", A.im());
 			this->setGenNtupleDoubleBranchValue("reB0barAmp", Abar.re());
 			this->setGenNtupleDoubleBranchValue("imB0barAmp", Abar.im());
 		} else {
 			this->setGenNtupleDoubleBranchValue("reB0Amp", 0.0);
 			this->setGenNtupleDoubleBranchValue("imB0Amp", 0.0);
 			this->setGenNtupleDoubleBranchValue("reB0barAmp", 0.0);
 			this->setGenNtupleDoubleBranchValue("imB0barAmp", 0.0);
 		}
 	}
 }
 
 void LauTimeDepFitModel::generateExtraPdfValues(LauPdfList* extraPdfs, LauEmbeddedData* embeddedData)
 {
     // CONVENTION WARNING
 	LauKinematics* kinematics(0);
 	if (curEvtTagFlv_<0) {
 		kinematics = kinematicsB0_;
 	} else {
 		kinematics = kinematicsB0bar_;
 	}
 
 	// Generate from the extra PDFs
 	if (extraPdfs) {
 		for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
 			LauFitData genValues;
 			if (embeddedData) {
 				genValues = embeddedData->getValues( (*pdf_iter)->varNames() );
 			} else {
 				genValues = (*pdf_iter)->generate(kinematics);
 			}
 			for ( LauFitData::const_iterator var_iter = genValues.begin(); var_iter != genValues.end(); ++var_iter ) {
 				TString varName = var_iter->first;
 				if ( varName != "m13Sq" && varName != "m23Sq" ) {
 					Double_t value = var_iter->second;
 					this->setGenNtupleDoubleBranchValue(varName,value);
 				}
 			}
 		}
 	}
 }
 
 void LauTimeDepFitModel::propagateParUpdates()
 {
 	// Update the complex mixing phase
 	if (useSinCos_) {
 		phiMixComplex_.setRealPart(cosPhiMix_.unblindValue());
 		phiMixComplex_.setImagPart(-1.0*sinPhiMix_.unblindValue());
 	} else {
 		phiMixComplex_.setRealPart(TMath::Cos(-1.0*phiMix_.unblindValue()));
 		phiMixComplex_.setImagPart(TMath::Sin(-1.0*phiMix_.unblindValue()));
 	}
 
 	// Update the total normalisation for the signal likelihood
 	if (this->useDP() == kTRUE) {
 		this->updateCoeffs();
 		sigModelB0bar_->updateCoeffs(coeffsB0bar_);
 		sigModelB0_->updateCoeffs(coeffsB0_);
 		this->calcInterTermNorm();
 	}
 
 	// Update the signal events from the background numbers if not doing an extended fit
 	// And update the tagging category fractions
 	this->updateSigEvents();
 }
 
 void LauTimeDepFitModel::updateSigEvents()
 {
 	// The background parameters will have been set from Minuit.
 	// We need to update the signal events using these.
 	if (!this->doEMLFit()) {
 		Double_t nTotEvts = this->eventsPerExpt();
 		Double_t signalEvents = nTotEvts;
 		// TODO loop over the background yields and subtract
 		signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts);
 		if ( ! signalEvents_->fixed() ) {
 			signalEvents_->value(signalEvents);
 		}
 	}
 
 	// tagging-category fractions for signal events
 	flavTag_->setFirstTagCatFractions();
 }
 
 void LauTimeDepFitModel::cacheInputFitVars()
 {
 	// Fill the internal data trees of the signal and background models.
 	// Note that we store the events of both charges in both the
 	// negative and the positive models. It's only later, at the stage
 	// when the likelihood is being calculated, that we separate them.
 
 	LauFitDataTree* inputFitData = this->fitData();
 
 	evtCPEigenVals_.clear();
 
 	const Bool_t hasCPEV = ( (cpevVarName_ != "") && inputFitData->haveBranch( cpevVarName_ ) );
 
 	UInt_t nEvents = inputFitData->nEvents();
 	evtCPEigenVals_.reserve( nEvents );
 	LauFitData::const_iterator fitdata_iter;
 	for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
 		const LauFitData& dataValues = inputFitData->getData(iEvt);
 
 		// if the CP-eigenvalue is in the data use those, otherwise use the default
 		if ( hasCPEV ) {
 			fitdata_iter = dataValues.find( cpevVarName_ );
 			const Int_t cpEV = static_cast<Int_t>( fitdata_iter->second );
 			if ( cpEV == 1 ) {
 				cpEigenValue_ = CPEven;
 			} else if ( cpEV == -1 ) {
 				cpEigenValue_ = CPOdd;
 			} else {
 				std::cerr<<"WARNING in LauTimeDepFitModel::cacheInputFitVars : Unknown value: "<<cpEV<<" for CP eigenvalue, setting to CP-even"<<std::endl;
 				cpEigenValue_ = CPEven;
 			}
 		}
 		evtCPEigenVals_.push_back( cpEigenValue_ );
 
 	}
 
 	// We'll cache the DP amplitudes at the end because we'll
 	// append some points that the other PDFs won't deal with.
 
 	// Flavour tagging information
     	flavTag_->cacheInputFitVars(inputFitData);
 
 	if (this->useDP() == kTRUE) {
 		// DecayTime and SigmaDecayTime
 		for (LauTagCatDtPdfMap::iterator dt_iter = signalDecayTimePdfs_.begin(); dt_iter != signalDecayTimePdfs_.end(); ++dt_iter) {
 			(*dt_iter).second->setFlavourTagging(flavTag_);
 			(*dt_iter).second->cacheInfo(*inputFitData);
 		}
 	}
 
 	// ...and then the extra PDFs
 	for (LauTagCatPdfListMap::iterator pdf_iter = sigExtraPdf_.begin(); pdf_iter != sigExtraPdf_.end(); ++pdf_iter) {
 		this->cacheInfo(pdf_iter->second, *inputFitData);
 	}
 
 	if (this->useDP() == kTRUE) {
 		sigModelB0bar_->fillDataTree(*inputFitData);
 		sigModelB0_->fillDataTree(*inputFitData);
 	}
 }
 
 Double_t LauTimeDepFitModel::getTotEvtLikelihood(const UInt_t iEvt)
 {
 
 	// Find out whether the tag-side B was a B0 or a B0bar.
 	curEvtTagFlv_ = flavTag_->getEvtTagFlvVals(iEvt);
 
 	// Also get the tagging category.
 	curEvtTagCat_ = flavTag_->getEvtTagCatVals(iEvt);
 
 	// Get the CP eigenvalue of the current event
 	cpEigenValue_ = evtCPEigenVals_[iEvt];
 
 	// Get the DP and DecayTime likelihood for signal (TODO and eventually backgrounds)
 	this->getEvtDPDtLikelihood(iEvt);
 
-	// Get the flavour tagging likelihood from eta PDFs (per tagging category - TODO backgrounds to come later) 
-	this->getEvtFlavTagLikelihood(iEvt);
+	// Get the flavour tagging likelihood from eta PDFs (per tagging category - TODO backgrounds to come later)
+	sigFlavTagLike_ = 1.0;
+	//this->getEvtFlavTagLikelihood(iEvt);
 
 	// Get the combined extra PDFs likelihood for signal (TODO and eventually backgrounds)
 	this->getEvtExtraLikelihoods(iEvt);
 
 	// Construct the total likelihood for signal, qqbar and BBbar backgrounds
 	Double_t sigLike = sigDPLike_ * sigFlavTagLike_ * sigExtraLike_;
 	//std::cout << "DP like = " << sigDPLike_ << std::endl;
 	//std::cout << "flav tag like = " << sigFlavTagLike_ << std::endl;
 	//std::cout << "extra like = " << sigExtraLike_ << std::endl;
 
 	Double_t signalEvents = signalEvents_->unblindValue();
 	if (this->useDP() == kFALSE) {
 		signalEvents *= 0.5 * (1.0 + curEvtTagFlv_ * signalAsym_->unblindValue());
 	}
 
 	// TODO better to store this info for each event
 	LauTagCatParamMap& signalTagCatFrac = flavTag_->getSignalTagCatFrac();
 	const Double_t sigTagCatFrac = signalTagCatFrac[curEvtTagCat_].unblindValue();
 	//std::cout << "Signal tag cat frac = " << sigTagCatFrac << std::endl;
 
 	// Construct the total event likelihood
 	Double_t likelihood(sigLike*sigTagCatFrac);
 	//std::cout << "Likelihoood = " << likelihood << std::endl;
 
 	if ( ! signalEvents_->fixed() ) {
 		likelihood *= signalEvents;
 	}
 
 	return likelihood;
 }
 
 Double_t LauTimeDepFitModel::getEventSum() const
 {
 	Double_t eventSum(0.0);
 	eventSum += signalEvents_->unblindValue();
 	return eventSum;
 }
 
 void LauTimeDepFitModel::getEvtDPDtLikelihood(const UInt_t iEvt)
 {
 	// Function to return the signal and background likelihoods for the
 	// Dalitz plot for the given event evtNo.
 	sigDPLike_ = 1.0;	//There's always a likelihood term for signal, so we better not zero it.
 
 	if ( this->useDP() == kFALSE ) {
 		return;
 	}
 
 	// Mistag probabilities. Defined as: omega = prob of the tagging B0 being reported as B0bar
 	// Whether we want omega or omegaBar depends on q_tag, hence curEvtTagFlv_*... in the previous lines
 	//Double_t misTagFrac = 0.5 * (1.0 - dilution_[curEvtTagCat_] - qDDo2);
 	//Double_t misTagFracBar = 0.5 * (1.0 - dilution_[curEvtTagCat_] + qDDo2);
 
 	// Calculate event quantities
 	Double_t qD(0);
 	Double_t qDDo2(0);
 	if (flavTag_->getUsePerEvtMistag()){
 		//qDDo2 term accounted for automatically with per event information
 		Double_t omega = flavTag_->getOmega(iEvt);
 		qD = curEvtTagFlv_*(1.0-2.0*omega);
 	} else {
 		// TODO need to sort out references here
 		LauTagCatParamMap dilution_ = flavTag_->getDilution();
 		LauTagCatParamMap deltaDilution_ = flavTag_->getDeltaDilution();
 
 		qD = curEvtTagFlv_*dilution_[curEvtTagCat_].unblindValue();
 		qDDo2 = curEvtTagFlv_*0.5*deltaDilution_[curEvtTagCat_].unblindValue();
 	}
 
 	// Get the dynamics to calculate everything required for the likelihood calculation
 	sigModelB0bar_->calcLikelihoodInfo(iEvt);
 	sigModelB0_->calcLikelihoodInfo(iEvt);
 
 	// Retrieve the amplitudes and efficiency from the dynamics
 	const LauComplex& Abar = sigModelB0bar_->getEvtDPAmp();
 	const LauComplex& A = sigModelB0_->getEvtDPAmp();
 	Double_t eff = sigModelB0bar_->getEvtEff();
 
 	// Next calculate the DP terms
 	Double_t aSqSum = A.abs2() + Abar.abs2();
 	Double_t aSqDif = A.abs2() - Abar.abs2();
 	LauComplex inter = Abar * A.conj() * phiMixComplex_;
 	Double_t interTermIm = 2.0 * inter.im();
 	Double_t interTermRe = 2.0 * inter.re();
 
 	// First get all the decay time terms
 	LauDecayTimePdf* decayTimePdf = signalDecayTimePdfs_[curEvtTagCat_];
 	decayTimePdf->calcLikelihoodInfo(iEvt);
 
 	// Get the decay time acceptance
 	Double_t dtEff = decayTimePdf->getEffiTerm();
 
 	// First get all the decay time terms
 	Double_t dtCos = decayTimePdf->getCosTerm();
 	Double_t dtSin = decayTimePdf->getSinTerm();
 	Double_t dtCosh = decayTimePdf->getCoshTerm();
 	Double_t dtSinh = decayTimePdf->getSinhTerm();
                                 
 	//std::cout << dtCos << " " << dtSin << " " << dtCosh << " " << dtSinh << std::endl;
 
 	Double_t cosTerm  = dtCos * qD * aSqDif;
 	Double_t sinTerm  = dtSin * qD * interTermIm;
 	Double_t coshTerm = dtCosh * (1.0 + qDDo2) * aSqSum;
 	Double_t sinhTerm = dtSinh * (1.0 + qDDo2) * interTermRe;
 				
 	if ( cpEigenValue_ == CPOdd ) {
 		sinTerm  *= -1.0;
 		sinhTerm *= -1.0;
 	}
 
 	// ... to get the total and multiply by the efficiency (DP and decay time)
 	Double_t ASq = coshTerm + cosTerm - sinTerm + sinhTerm;
 	ASq *= eff;
 	ASq *= dtEff;
 
 	//std::cout<<"Cosh Cos Sin Sinh "<< coshTerm<< " " << cosTerm << " " << sinTerm << " " << sinhTerm << std::endl;
 	//std::cout << "Total Amplitude : " << ASq << std::endl;
 
 	// Calculate the DP and time normalisation
 	Double_t normTermIndep = sigModelB0bar_->getDPNorm() + sigModelB0_->getDPNorm();
 	Double_t normTermCosh(0.0); 
 	Double_t norm(0.0); 
 
         if (decayTimePdf->getFuncType() == LauDecayTimePdf::ExpTrig || decayTimePdf->getFuncType() == LauDecayTimePdf::SimFitNormBd || decayTimePdf->getFuncType() == LauDecayTimePdf::SimFitSigBd){ 
                 normTermCosh = decayTimePdf->getNormTermExp();
                 norm = normTermIndep*normTermCosh; 
         }
         if (decayTimePdf->getFuncType() == LauDecayTimePdf::ExpHypTrig){             
 	        normTermCosh = decayTimePdf->getNormTermCosh();
 		Double_t normTermDep   = interTermReNorm_;
 		Double_t normTermSinh  = decayTimePdf->getNormTermSinh();
                 norm = normTermIndep*normTermCosh + normTermDep*normTermSinh;
         }
 
 
 	// Calculate the normalised signal likelihood
 	//std::cout << "ASq = " << ASq << std::endl;
 	//std::cout << "norm = " << norm << std::endl;
 	sigDPLike_ = ASq / norm;
 }
 
 void LauTimeDepFitModel::getEvtExtraLikelihoods(const UInt_t iEvt)
 {
 	// Function to return the signal and background likelihoods for the
 	// extra variables for the given event evtNo.
 	sigExtraLike_ = 1.0;	//There's always a likelihood term for signal, so we better not zero it.
 
 	// First, those independent of the tagging of the event:
 	// signal
 	LauTagCatPdfListMap::iterator sig_iter = sigExtraPdf_.find(curEvtTagCat_);
 	LauPdfList* pdfList = (sig_iter != sigExtraPdf_.end())? &(sig_iter->second) : 0;
 	if (pdfList) {
 		sigExtraLike_ = this->prodPdfValue( *pdfList, iEvt );
 	}
 
 }
 
 void LauTimeDepFitModel::getEvtFlavTagLikelihood(const UInt_t iEvt)
 {
 	// Function to return the signal and background likelihoods for the
 	// extra variables for the given event evtNo.
 	sigFlavTagLike_ = 1.0;	//There's always a likelihood term for signal, so we better not zero it.
 
 	// First, those independent of the tagging of the event:
 	// signal
 	LauAbsPdf* pdf = sigFlavTagPdf_[curEvtTagCat_];
 	if (pdf) {
 		pdf->calcLikelihoodInfo(iEvt);
 		sigFlavTagLike_ = pdf->getLikelihood();
 	}
 
 	if (sigFlavTagLike_<=0){
 		std::cout<<"INFO in LauTimeDepFitModel::getEvtFlavTagLikelihood : Event with 0 FlavTag Liklihood"<<std::endl;
 	}
 	// TODO Add in the background components too
 
 }
 
 void LauTimeDepFitModel::updateCoeffs()
 {
 	coeffsB0bar_.clear(); coeffsB0_.clear();
 	coeffsB0bar_.reserve(nSigComp_); coeffsB0_.reserve(nSigComp_);
 	for (UInt_t i = 0; i < nSigComp_; ++i) {
 		coeffsB0bar_.push_back(coeffPars_[i]->antiparticleCoeff());
 		coeffsB0_.push_back(coeffPars_[i]->particleCoeff());
 	}
 }
 
 void LauTimeDepFitModel::checkMixingPhase()
 {
 	Double_t phase = phiMix_.value();
 	Double_t genPhase = phiMix_.genValue();
 
 	// Check now whether the phase lies in the right range (-pi to pi).
 	Bool_t withinRange(kFALSE);
 	while (withinRange == kFALSE) {
 		if (phase > -LauConstants::pi && phase < LauConstants::pi) {
 			withinRange = kTRUE;
 		} else {
 			// Not within the specified range
 			if (phase > LauConstants::pi) {
 				phase -= LauConstants::twoPi;
 			} else if (phase < -LauConstants::pi) {
 				phase += LauConstants::twoPi;
 			}
 		}
 	}
 
 	// A further problem can occur when the generated phase is close to -pi or pi.
 	// The phase can wrap over to the other end of the scale -
 	// this leads to artificially large pulls so we wrap it back.
 	Double_t diff = phase - genPhase;
 	if (diff > LauConstants::pi) {
 		phase -= LauConstants::twoPi;
 	} else if (diff < -LauConstants::pi) {
 		phase += LauConstants::twoPi;
 	}
 
 	// finally store the new value in the parameter
 	// and update the pull
 	phiMix_.value(phase);
 	phiMix_.updatePull();
 }
 
 void LauTimeDepFitModel::embedSignal(Int_t tagCat, const TString& fileName, const TString& treeName,
 					Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment)
 {
 	if (signalTree_[tagCat]) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::embedSignal : Already embedding signal from file for tagging category "<<tagCat<<"."<<std::endl;
 		return;
 	}
 
 	if (!reuseEventsWithinEnsemble && reuseEventsWithinExperiment) {
 		std::cerr<<"WARNING in LauTimeDepFitModel::embedSignal : Conflicting options provided, will not reuse events at all."<<std::endl;
 		reuseEventsWithinExperiment = kFALSE;
 	}
 
 	signalTree_[tagCat] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);
 	Bool_t dataOK = signalTree_[tagCat]->findBranches();
 	if (!dataOK) {
 		delete signalTree_[tagCat]; signalTree_[tagCat] = 0;
 		std::cerr<<"ERROR in LauTimeDepFitModel::embedSignal : Problem creating data tree for embedding."<<std::endl;
 		return;
 	}
 
 	reuseSignal_ = reuseEventsWithinEnsemble;
 }
 
 void LauTimeDepFitModel::setupSPlotNtupleBranches()
 {
 	// add branches for storing the experiment number and the number of
 	// the event within the current experiment
 	this->addSPlotNtupleIntegerBranch("iExpt");
 	this->addSPlotNtupleIntegerBranch("iEvtWithinExpt");
 
 	// Store the efficiency of the event (for inclusive BF calculations).
 	if (this->storeDPEff()) {
 		this->addSPlotNtupleDoubleBranch("efficiency");
 	}
 
 	// Store the total event likelihood for each species.
 	this->addSPlotNtupleDoubleBranch("sigTotalLike");
 
 	// Store the DP likelihoods
 	if (this->useDP()) {
 		this->addSPlotNtupleDoubleBranch("sigDPLike");
 	}
 
 	// Store the likelihoods for each extra PDF
 	const LauPdfList* pdfList( &(sigExtraPdf_.begin()->second) );
 	this->addSPlotNtupleBranches(pdfList, "sig");
 }
 
 void LauTimeDepFitModel::addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix)
 {
 	if (!extraPdfs) {
 		return;
 	}
 	// Loop through each of the PDFs
 	for (LauPdfList::const_iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
 
 		// Count the number of input variables that are not
 		// DP variables (used in the case where there is DP
 		// dependence for e.g. MVA)
 		UInt_t nVars(0);
 		for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
 			if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
 				++nVars;
 			}
 		}
 
 		if ( nVars == 1 ) {
 			// If the PDF only has one variable then
 			// simply add one branch for that variable
 			TString varName = (*pdf_iter)->varName();
 			TString name(prefix);
 			name += varName;
 			name += "Like";
 			this->addSPlotNtupleDoubleBranch(name);
 		} else if ( nVars == 2 ) {
 			// If the PDF has two variables then we
 			// need a branch for them both together and
 			// branches for each
 			TString allVars("");
 			for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
 				allVars += (*var_iter);
 				TString name(prefix);
 				name += (*var_iter);
 				name += "Like";
 				this->addSPlotNtupleDoubleBranch(name);
 			}
 			TString name(prefix);
 			name += allVars;
 			name += "Like";
 			this->addSPlotNtupleDoubleBranch(name);
 		} else {
 			std::cerr<<"WARNING in LauTimeDepFitModel::addSPlotNtupleBranches : Can't yet deal with 3D PDFs."<<std::endl;
 		}
 	}
 }
 
 Double_t LauTimeDepFitModel::setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt)
 {
 	// Store the PDF value for each variable in the list
 	Double_t totalLike(1.0);
 	Double_t extraLike(0.0);
 	if ( !extraPdfs ) {
 		return totalLike;
 	}
 
 	for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) {
 
 		// calculate the likelihood for this event
 		(*pdf_iter)->calcLikelihoodInfo(iEvt);
 		extraLike = (*pdf_iter)->getLikelihood();
 		totalLike *= extraLike;
 
 		// Count the number of input variables that are not
 		// DP variables (used in the case where there is DP
 		// dependence for e.g. MVA)
 		UInt_t nVars(0);
 		for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
 			if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
 				++nVars;
 			}
 		}
 
 		if ( nVars == 1 ) {
 			// If the PDF only has one variable then
 			// simply store the value for that variable
 			TString varName = (*pdf_iter)->varName();
 			TString name(prefix);
 			name += varName;
 			name += "Like";
 			this->setSPlotNtupleDoubleBranchValue(name, extraLike);
 		} else if ( nVars == 2 ) {
 			// If the PDF has two variables then we
 			// store the value for them both together
 			// and for each on their own
 			TString allVars("");
 			for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
 				allVars += (*var_iter);
 				TString name(prefix);
 				name += (*var_iter);
 				name += "Like";
 				Double_t indivLike = (*pdf_iter)->getLikelihood( (*var_iter) );
 				this->setSPlotNtupleDoubleBranchValue(name, indivLike);
 			}
 			TString name(prefix);
 			name += allVars;
 			name += "Like";
 			this->setSPlotNtupleDoubleBranchValue(name, extraLike);
 		} else {
 			std::cerr<<"WARNING in LauAllFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs."<<std::endl;
 		}
 	}
 	return totalLike;
 }
 
 LauSPlot::NameSet LauTimeDepFitModel::variableNames() const
 {
 	LauSPlot::NameSet nameSet;
 	if (this->useDP()) {
 		nameSet.insert("DP");
 	}
 	LauPdfList pdfList( (sigExtraPdf_.begin()->second) );
 	for (LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
 		// Loop over the variables involved in each PDF
 		for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
 			// If they are not DP coordinates then add them
 			if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
 				nameSet.insert( (*var_iter) );
 			}
 		}
 	}
 	return nameSet;
 }
 
 LauSPlot::NumbMap LauTimeDepFitModel::freeSpeciesNames() const
 {
 	LauSPlot::NumbMap numbMap;
 	if (!signalEvents_->fixed() && this->doEMLFit()) {
 		numbMap["sig"] = signalEvents_->genValue();
 	}
 	return numbMap;
 }
 
 LauSPlot::NumbMap LauTimeDepFitModel::fixdSpeciesNames() const
 {
 	LauSPlot::NumbMap numbMap;
 	if (signalEvents_->fixed() && this->doEMLFit()) {
 		numbMap["sig"] = signalEvents_->genValue();
 	}
 	return numbMap;
 }
 
 LauSPlot::TwoDMap LauTimeDepFitModel::twodimPDFs() const
 {
 	LauSPlot::TwoDMap twodimMap;
 
 	const LauPdfList* pdfList =  &(sigExtraPdf_.begin()->second);
 	for (LauPdfList::const_iterator pdf_iter = pdfList->begin(); pdf_iter != pdfList->end(); ++pdf_iter) {
 		// Count the number of input variables that are not DP variables
 		UInt_t nVars(0);
 		for ( std::vector<TString>::const_iterator var_iter = (*pdf_iter)->varNames().begin(); var_iter != (*pdf_iter)->varNames().end(); ++var_iter ) {
 			if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
 				++nVars;
 			}
 		}
 		if ( nVars == 2 ) {
 			twodimMap.insert( std::make_pair( "sig", std::make_pair( (*pdf_iter)->varNames()[0], (*pdf_iter)->varNames()[1] ) ) );
 		}
 	}
 
 	return twodimMap;
 }
 
 void LauTimeDepFitModel::storePerEvtLlhds()
 {
 	std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Storing per-event likelihood values..."<<std::endl;
 
 	LauFitDataTree* inputFitData = this->fitData();
 
 	// if we've not been using the DP model then we need to cache all
 	// the info here so that we can get the efficiency from it
 	if (!this->useDP() && this->storeDPEff()) {
 		sigModelB0bar_->initialise(coeffsB0bar_);
 		sigModelB0_->initialise(coeffsB0_);
 		sigModelB0bar_->fillDataTree(*inputFitData);
 		sigModelB0_->fillDataTree(*inputFitData);
 	}
 
 	UInt_t evtsPerExpt(this->eventsPerExpt());
 
 	LauIsobarDynamics* sigModel(sigModelB0bar_);
 
 	for (UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt) {
 		// Find out whether we have B0bar or B0
 		curEvtTagFlv_ = flavTag_->getEvtTagFlvVals(iEvt);
 		curEvtTagCat_ = flavTag_->getEvtTagCatVals(iEvt);
 		curEvtMistag_ = flavTag_->getOmega(iEvt);
 
 		LauTagCatPdfListMap::iterator sig_iter = sigExtraPdf_.find(curEvtTagCat_);
 		LauPdfList* sigPdfs = (sig_iter != sigExtraPdf_.end())? &(sig_iter->second) : 0;
 
 		// the DP information
 		this->getEvtDPDtLikelihood(iEvt);
 		if (this->storeDPEff()) {
 			if (!this->useDP()) {
 				sigModel->calcLikelihoodInfo(iEvt);
 			}
 			this->setSPlotNtupleDoubleBranchValue("efficiency",sigModel->getEvtEff());
 		}
 		if (this->useDP()) {
 			sigTotalLike_ = sigDPLike_;
 			this->setSPlotNtupleDoubleBranchValue("sigDPLike",sigDPLike_);
 		} else {
 			sigTotalLike_ = 1.0;
 		}
 
 		// the signal PDF values
 		sigTotalLike_ *= this->setSPlotNtupleBranchValues(sigPdfs, "sig", iEvt);
 
 		// the total likelihoods
 		this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_);
 
 		// fill the tree
 		this->fillSPlotNtupleBranches();
 	}
 	std::cout<<"INFO in LauTimeDepFitModel::storePerEvtLlhds : Finished storing per-event likelihood values."<<std::endl;
 }
 
 void LauTimeDepFitModel::weightEvents( const TString& /*dataFileName*/, const TString& /*dataTreeName*/ )
 {
 	std::cerr << "ERROR in LauTimeDepFitModel::weightEvents : Method not available for this fit model." << std::endl;
 	return;
 }
 
 void LauTimeDepFitModel::savePDFPlots(const TString& /*label*/)
 {
 }
 
 void LauTimeDepFitModel::savePDFPlotsWave(const TString& /*label*/, const Int_t& /*spin*/)
 {
 }
 
 void LauTimeDepFitModel::setSignalFlavTagPdfs(const Int_t tagCat, LauAbsPdf* pdf)
 {
 	if (!flavTag_->validTagCat(tagCat)) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::setSignalFlavTagPdfs : Tagging category \""<<tagCat<<"\" not valid, not adding models."<<std::endl;
 		return;
 	}
 	if (pdf==0) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::setSignalFlavTagPdfs : The PDF pointer is null."<<std::endl;
 		return;
 	}
 
 	sigFlavTagPdf_[tagCat] = pdf;
 }
 
 void LauTimeDepFitModel::setBkgdFlavTagPdfs(const TString name, LauAbsPdf* pdf)
 {
 	// TODO - when backgrounds are added - implement a check that "name" really is the name of a background category
 	if (pdf==0) {
 		std::cerr<<"ERROR in LauTimeDepFitModel::setBkgdFlavTagPdfs : The PDF pointer is null."<<std::endl;
 		return;
 	}
 
 	bkgdFlavTagPdf_[name] = pdf;
 }