LauSimpleFitModel.cc

Go to the documentation of this file.

/*
Copyright 2004 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
*/

#include <iostream>
#include <iomanip>
#include <fstream>
#include <typeinfo>

#include "TFile.h"
#include "TH2.h"
#include "TMinuit.h"
#include "TRandom.h"
#include "TSystem.h"
#include "TVirtualFitter.h"

#include "LauAbsBkgndDPModel.hh"
#include "LauAbsCoeffSet.hh"
#include "LauIsobarDynamics.hh"
#include "LauAbsPdf.hh"
#include "LauComplex.hh"
#include "LauConstants.hh"
#include "LauDaughters.hh"
#include "LauEmbeddedData.hh"
#include "LauEffModel.hh"
#include "LauFitNtuple.hh"
#include "LauGenNtuple.hh"
#include "LauKinematics.hh"
#include "LauPrint.hh"
#include "LauRandom.hh"
#include "LauScfMap.hh"
#include "LauSimpleFitModel.hh"
#include "TGraph2D.h"
#include "TGraph.h"
#include "TStyle.h"
#include "TCanvas.h"

ClassImp(LauSimpleFitModel)


LauSimpleFitModel::LauSimpleFitModel(LauIsobarDynamics* sigModel) : LauAbsFitModel(),
        sigDPModel_(sigModel),
        kinematics_(sigModel ? sigModel->getKinematics() : 0),
        usingBkgnd_(kFALSE),
        nSigComp_(0),
        nSigDPPar_(0),
        nExtraPdfPar_(0),
        nNormPar_(0),
        meanEff_("meanEff",0.0,0.0,1.0),
        dpRate_("dpRate",0.0,0.0,100.0),
        //signalEvents_("signalEvents",1.0,0.0,1.0),
        signalEvents_(0),
        useSCF_(kFALSE),
        useSCFHist_(kFALSE),
        scfFrac_("scfFrac",0.0,0.0,1.0),
        scfFracHist_(0),
        scfMap_(0),
        compareFitData_(kFALSE),
        signalTree_(0),
        reuseSignal_(kFALSE),
        useReweighting_(kFALSE),
        sigDPLike_(0.0),
        scfDPLike_(0.0),
        sigExtraLike_(0.0),
        scfExtraLike_(0.0),
        sigTotalLike_(0.0),
        scfTotalLike_(0.0)
{
}

LauSimpleFitModel::~LauSimpleFitModel()
{
        delete signalTree_;
        for (LauBkgndEmbDataList::iterator iter = bkgndTree_.begin(); iter != bkgndTree_.end(); ++iter) {
                delete (*iter);
        }
        delete scfFracHist_;
}

void LauSimpleFitModel::setupBkgndVectors()
{
        UInt_t nBkgnds = this->nBkgndClasses();
        bkgndDPModels_.resize( nBkgnds );
        bkgndPdfs_.resize( nBkgnds );
        bkgndEvents_.resize( nBkgnds );
        bkgndTree_.resize( nBkgnds );
        reuseBkgnd_.resize( nBkgnds );
        bkgndDPLike_.resize( nBkgnds );
        bkgndExtraLike_.resize( nBkgnds );
        bkgndTotalLike_.resize( nBkgnds );
}

void LauSimpleFitModel::setNSigEvents(LauParameter* nSigEvents)
{
        if ( nSigEvents == 0 ) {
                std::cerr << "ERROR in LauSimpleFitModel::setNSigEvents : The signal yield LauParameter pointer is null." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }
        if ( signalEvents_ != 0 ) {
                std::cerr << "ERROR in LauSimpleFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
                return;
        }

        signalEvents_ = nSigEvents;
        TString name = signalEvents_->name();
        if ( ! name.Contains("signalEvents") && !( name.BeginsWith("signal") && name.EndsWith("Events") ) ) {
                signalEvents_->name("signalEvents");
        }
        Double_t value = nSigEvents->value();
        signalEvents_->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0));
}

void LauSimpleFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents)
{
        if ( nBkgndEvents == 0 ) {
                std::cerr << "ERROR in LauSimpleFitModel::setNBkgndEvents : The background yield LauParameter pointer is null." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }

        if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) {
                std::cerr << "ERROR in LauSimpleFitModel::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 LauSimpleFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << 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;
}

void LauSimpleFitModel::splitSignalComponent( const TH2* dpHisto, const Bool_t upperHalf, const Bool_t fluctuateBins, LauScfMap* scfMap )
{
        if ( useSCF_ == kTRUE ) {
                std::cerr << "ERROR in LauSimpleFitModel::splitSignalComponent : Have already setup SCF." << std::endl;
                return;
        }

        if ( dpHisto == 0 ) {
                std::cerr << "ERROR in LauSimpleFitModel::splitSignalComponent : The histogram pointer is null." << std::endl;
                return;
        }

        const LauDaughters* daughters = sigDPModel_->getDaughters();
        scfFracHist_ = new LauEffModel( daughters, 0 );
        scfFracHist_->setEffHisto( dpHisto, kTRUE, fluctuateBins, 0.0, 0.0, upperHalf, daughters->squareDP() );

        scfMap_ = scfMap;

        useSCF_ = kTRUE;
        useSCFHist_ = kTRUE;
}

void LauSimpleFitModel::splitSignalComponent( const Double_t scfFrac, const Bool_t fixed )
{
        if ( useSCF_ == kTRUE ) {
                std::cerr << "ERROR in LauSimpleFitModel::splitSignalComponent : Have already setup SCF." << std::endl;
                return;
        }

        scfFrac_.range( 0.0, 1.0 );
        scfFrac_.value( scfFrac ); scfFrac_.initValue( scfFrac ); scfFrac_.genValue( scfFrac );
        scfFrac_.fixed( fixed );

        useSCF_ = kTRUE;
        useSCFHist_ = kFALSE;
}

void LauSimpleFitModel::setBkgndDPModel(const TString& bkgndClass, LauAbsBkgndDPModel* bkgndDPModel)
{
        if (bkgndDPModel == 0) {
                std::cerr << "ERROR in LauSimpleFitModel::setBkgndDPModel : The model pointer is null." << std::endl;
                return;
        }

        // check that this background name is valid
        if ( ! this->validBkgndClass( bkgndClass ) ) {
                std::cerr << "ERROR in LauSimpleFitModel::setBkgndDPModel : Invalid background class \"" << bkgndClass << "\"." << std::endl;
                std::cerr << "                                            : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
                return;
        }

        UInt_t bkgndID = this->bkgndClassID( bkgndClass );
        bkgndDPModels_[bkgndID] = bkgndDPModel;

        usingBkgnd_ = kTRUE;
}

void LauSimpleFitModel::setSignalPdf(LauAbsPdf* pdf)
{
        if (pdf==0) {
                std::cerr << "ERROR in LauSimpleFitModel::setSignalPdf : The PDF pointer is null." << std::endl;
                return;
        }
        signalPdfs_.push_back(pdf);
}

void LauSimpleFitModel::setSCFPdf(LauAbsPdf* pdf)
{
        if (pdf==0) {
                std::cerr << "ERROR in LauSimpleFitModel::setSCFPdf : The PDF pointer is null." << std::endl;
                return;
        }
        scfPdfs_.push_back(pdf);
}

void LauSimpleFitModel::setBkgndPdf(const TString& bkgndClass, LauAbsPdf* pdf)
{
        if (pdf == 0) {
                std::cerr << "ERROR in LauSimpleFitModel::setBkgndPdf : The PDF pointer is null." << std::endl;
                return;
        }

        // check that this background name is valid
        if ( ! this->validBkgndClass( bkgndClass ) ) {
                std::cerr << "ERROR in LauSimpleFitModel::setBkgndPdf : Invalid background class \"" << bkgndClass << "\"." << std::endl;
                std::cerr << "                                        : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
                return;
        }

        UInt_t bkgndID = this->bkgndClassID( bkgndClass );
        bkgndPdfs_[bkgndID].push_back(pdf);

        usingBkgnd_ = kTRUE;
}

void LauSimpleFitModel::setAmpCoeffSet(LauAbsCoeffSet* coeffSet)
{
        // Resize the coeffPars vector if not already done
        if ( coeffPars_.empty() ) {
                coeffPars_.resize( sigDPModel_->getnTotAmp() );
                for (std::vector<LauAbsCoeffSet*>::iterator iter = coeffPars_.begin(); iter != coeffPars_.end(); ++iter) {
                        (*iter) = 0;
                }
        }

        // Is there a component called compName in the signal model?
        TString compName(coeffSet->name());
        const Int_t index = sigDPModel_->resonanceIndex(compName);
        if ( index < 0 ) {
                std::cerr << "ERROR in LauSimpleFitModel::setAmpCoeffSet : Signal DP model doesn't contain component \"" << compName << "\"." << std::endl;
                return;
        }

        // Do we already have it in our list of names?
        if ( coeffPars_[index] != 0 && coeffPars_[index]->name() == compName) {
                std::cerr << "ERROR in LauSimpleFitModel::setAmpCoeffSet : Have already set coefficients for \"" << compName << "\"." << std::endl;
                return;
        }

        coeffSet->index(index);
        coeffPars_[index] = coeffSet;

        ++nSigComp_;

        std::cout << "INFO in LauSimpleFitModel::setAmpCoeffSet : Added coefficients for component A"<< index << ": \"" << compName << "\" to the fit model." << std::endl;
        coeffSet->printParValues();
}


void LauSimpleFitModel::initialise()
{
        // Initialisation
        if (!this->useDP() && signalPdfs_.empty()) {
                std::cerr << "ERROR in LauSimpleFitModel::initialise : Signal model doesn't exist for any variable." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }

        if (this->useDP()) {
                // Check that we have all the Dalitz-plot models
                if ( sigDPModel_ == 0 ) {
                        std::cerr << "ERROR in LauSimpleFitModel::initialise : The pointer to the signal DP model is null.\n";
                        std::cerr << "                                       : Removing the Dalitz Plot from the model." << std::endl;
                        this->useDP(kFALSE);
                }
                if ( usingBkgnd_ ) {
                        for (LauBkgndDPModelList::const_iterator dpmodel_iter = bkgndDPModels_.begin(); dpmodel_iter != bkgndDPModels_.end(); ++dpmodel_iter ) {
                                if ( (*dpmodel_iter) == 0 ) {
                                        std::cerr << "ERROR in LauSimpleFitModel::initialise : The pointer to one of the background DP models is null.\n";
                                        std::cerr << "                                       : Removing the Dalitz Plot from the model." << std::endl;
                                        this->useDP(kFALSE);
                                        break;
                                }
                        }
                }

                // Need to check that the number of components we have and that the dynamics has matches up
                UInt_t nAmp = sigDPModel_->getnTotAmp();
                if (nAmp != nSigComp_) {
                        std::cerr << "ERROR in LauSimpleFitModel::initialise : Number of signal DP components with magnitude and phase set not right." << std::endl;
                        gSystem->Exit(EXIT_FAILURE);
                }

                // From the initial parameter values calculate the coefficients
                // so they can be passed to the signal model
                this->updateCoeffs();

                // If all is well, go ahead and initialise them
                this->initialiseDPModels();
        }

        // Next check that, if a given component is being used we've got the
        // right number of PDFs for all the variables involved
        // TODO - should probably check variable names and so on as well

        UInt_t nsigpdfvars(0);
        for ( LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter ) {
                std::vector<TString> varNames = (*pdf_iter)->varNames();
                for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
                        if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
                                ++nsigpdfvars;
                        }
                }
        }
        if (useSCF_) {
                UInt_t nscfpdfvars(0);
                for ( LauPdfList::const_iterator pdf_iter = scfPdfs_.begin(); pdf_iter != scfPdfs_.end(); ++pdf_iter ) {
                        std::vector<TString> varNames = (*pdf_iter)->varNames();
                        for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
                                if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
                                        ++nscfpdfvars;
                                }
                        }
                }
                if (nscfpdfvars != nsigpdfvars) {
                        std::cerr << "ERROR in LauSimpleFitModel::initialise : There are " << nsigpdfvars << " TM signal PDF variables but " << nscfpdfvars << " SCF signal PDF variables." << std::endl;
                        gSystem->Exit(EXIT_FAILURE);
                }
        }
        if (usingBkgnd_) {
                for (LauBkgndPdfsList::const_iterator bgclass_iter = bkgndPdfs_.begin(); bgclass_iter != bkgndPdfs_.end(); ++bgclass_iter) {
                        UInt_t nbkgndpdfvars(0);
                        const LauPdfList& pdfList = (*bgclass_iter);
                        for ( LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter ) {
                                std::vector<TString> varNames = (*pdf_iter)->varNames();
                                for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
                                        if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
                                                ++nbkgndpdfvars;
                                        }
                                }
                        }
                        if (nbkgndpdfvars != nsigpdfvars) {
                                std::cerr << "ERROR in LauSimpleFitModel::initialise : There are " << nsigpdfvars << " signal PDF variables but " << nbkgndpdfvars << " bkgnd PDF variables." << 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 various extra PDFs
        this->setExtraPdfParameters();

        // Set the initial bg and signal events
        this->setFitNEvents();

        // Check that we have the expected number of fit variables
        const LauParameterPList& fitVars = this->fitPars();
        if (fitVars.size() != (nSigDPPar_ + nExtraPdfPar_ + nNormPar_)) {
                std::cerr << "ERROR in LauSimpleFitModel::initialise : Number of fit parameters not of expected size." << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }

        this->setExtraNtupleVars();
}

void LauSimpleFitModel::initialiseDPModels()
{
        std::cout << "INFO in LauSimpleFitModel::initialiseDPModels : Initialising DP models" << std::endl;

        sigDPModel_->initialise(coeffs_);

        if (usingBkgnd_ == kTRUE) {
                for (LauBkgndDPModelList::iterator iter = bkgndDPModels_.begin(); iter != bkgndDPModels_.end(); ++iter) {
                        (*iter)->initialise();
                }
        }
}

void LauSimpleFitModel::recalculateNormalisation()
{
        //std::cout << "INFO in LauSimpleFitModel::recalculateNormalizationInDPModels : Recalc Norm in DP model" << std::endl;
        sigDPModel_->recalculateNormalisation();
        sigDPModel_->modifyDataTree();
}

void LauSimpleFitModel::setSignalDPParameters()
{
        // Set the fit parameters for the signal model.

        nSigDPPar_ = 0;
        if ( ! this->useDP() ) {
                return;
        }

        std::cout << "INFO in LauSimpleFitModel::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
        LauParameterPSet& resVars = this->resPars();
        resVars.clear();

        LauParameterPList& sigDPPars = sigDPModel_->getFloatingParameters();

        for ( LauParameterPList::iterator iter = sigDPPars.begin(); iter != sigDPPars.end(); ++iter ) {
                if ( resVars.insert(*iter).second ) {
                        fitVars.push_back(*iter);
                        ++nSigDPPar_;
                }
        }
}

void LauSimpleFitModel::setExtraPdfParameters()
{
        // Include all the parameters of the various PDFs 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->addFitParameters(signalPdfs_);

        if (useSCF_ == kTRUE) {
                nExtraPdfPar_ += this->addFitParameters(scfPdfs_);
        }

        if (usingBkgnd_ == kTRUE) {
                for (LauBkgndPdfsList::iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) {
                        nExtraPdfPar_ += this->addFitParameters(*iter);
                }
        }
}

void LauSimpleFitModel::setFitNEvents()
{
        if ( signalEvents_ == 0 ) {
                std::cerr << "ERROR in LauSimpleFitModel::setFitNEvents : Signal yield not defined." << std::endl;
                return;
        }
        nNormPar_ = 0;

        // initialise the total number of events to be the sum of all the hypotheses
        Double_t nTotEvts = signalEvents_->value();
        for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
                nTotEvts += (*iter)->value();
                if ( (*iter) == 0 ) {
                        std::cerr << "ERROR in LauSimpleFitModel::setFitNEvents : Background yield not defined." << std::endl;
                        return;
                }
        }
        this->eventsPerExpt(TMath::FloorNint(nTotEvts));

        LauParameterPList& fitVars = this->fitPars();

        // if doing an extended ML fit add the number of signal events into the fit parameters
        if (this->doEMLFit()) {
                std::cout << "INFO in LauSimpleFitModel::setFitNEvents : Initialising number of events for signal and background components..." << std::endl;
                // add the signal events to the list of fit parameters
                fitVars.push_back(signalEvents_);
                ++nNormPar_;
        } else {
                std::cout << "INFO in LauSimpleFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..." << std::endl;
        }

        if (useSCF_ && !useSCFHist_) {
                fitVars.push_back(&scfFrac_);
                ++nNormPar_;
        }

        if (usingBkgnd_ == kTRUE) {
                for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
                        std::vector<LauParameter*> parameters = (*iter)->getPars();
                        for ( LauParameter* parameter : parameters ) {
                                if ( ! parameter->clone() ) {
                                        fitVars.push_back(parameter);
                                        ++nNormPar_;
                                }
                        }
                }
        }
}

void LauSimpleFitModel::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 a fit fraction for each signal component
        fitFrac_ = sigDPModel_->getFitFractions();
        if (fitFrac_.size() != nSigComp_) {
                std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFrac_.size() << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }
        for (UInt_t i(0); i<nSigComp_; ++i) {
                if (fitFrac_[i].size() != nSigComp_) {
                        std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFrac_[i].size() << std::endl;
                        gSystem->Exit(EXIT_FAILURE);
                }
        }

        // Add the fit fraction that has not been corrected for the efficiency for each signal component
        fitFracEffUnCorr_ = sigDPModel_->getFitFractionsEfficiencyUncorrected();
        if (fitFracEffUnCorr_.size() != nSigComp_) {
                std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr_.size() << std::endl;
                gSystem->Exit(EXIT_FAILURE);
        }
        for (UInt_t i(0); i<nSigComp_; ++i) {
                if (fitFracEffUnCorr_[i].size() != nSigComp_) {
                        std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr_[i].size() << std::endl;
                        gSystem->Exit(EXIT_FAILURE);
                }
        }

        for (UInt_t i = 0; i < nSigComp_; i++) {
                for (UInt_t j = i; j < nSigComp_; j++) {
                        extraVars.push_back(fitFrac_[i][j]);
                        extraVars.push_back(fitFracEffUnCorr_[i][j]);
                }
        }

        // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
        std::vector<LauParameter> extraParams = sigDPModel_->getExtraParameters();
        std::vector<LauParameter>::iterator extraIter;
        for (extraIter = extraParams.begin(); extraIter != extraParams.end(); ++extraIter) {
                LauParameter extraParameter = (*extraIter);
                extraVars.push_back(extraParameter);
        }

        // Now add in the DP efficiency value
        Double_t initMeanEff = sigDPModel_->getMeanEff().initValue();
        meanEff_.value(initMeanEff); meanEff_.initValue(initMeanEff); meanEff_.genValue(initMeanEff);
        extraVars.push_back(meanEff_);

        // Also add in the DP rate
        Double_t initDPRate = sigDPModel_->getDPRate().initValue();
        dpRate_.value(initDPRate); dpRate_.initValue(initDPRate); dpRate_.genValue(initDPRate);
        extraVars.push_back(dpRate_);
}

void LauSimpleFitModel::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
        // e.g. to make mag > 0.0 and 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 mag > 0.0, and phase in the right range
                        coeffPars_[i]->finaliseValues();
                }
        }

        // update the pulls on the events
        if (this->doEMLFit()) {
                signalEvents_->updatePull();
        }
        if (useSCF_ && !useSCFHist_) {
                scfFrac_.updatePull();
        }
        if (usingBkgnd_ == kTRUE) {
                for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
                        std::vector<LauParameter*> parameters = (*iter)->getPars();
                        for ( LauParameter* parameter : parameters ) {
                                parameter->updatePull();
                        }
                }
        }

        // Update the pulls on all the extra PDFs' parameters
        this->updateFitParameters(signalPdfs_);
        if (useSCF_ == kTRUE) {
                this->updateFitParameters(scfPdfs_);
        }
        if (usingBkgnd_ == kTRUE) {
                for (LauBkgndPdfsList::iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) {
                        this->updateFitParameters(*iter);
                }
        }

        // Fill the fit results to the ntuple for current experiment

        // update the coefficients and then calculate the fit fractions
        if (this->useDP() == kTRUE) {
                this->updateCoeffs();
                sigDPModel_->updateCoeffs(coeffs_);
                sigDPModel_->calcExtraInfo();

                LauParArray fitFrac = sigDPModel_->getFitFractions();
                if (fitFrac.size() != nSigComp_) {
                        std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFrac.size() << std::endl;
                        gSystem->Exit(EXIT_FAILURE);
                }
                for (UInt_t i(0); i<nSigComp_; ++i) {
                        if (fitFrac[i].size() != nSigComp_) {
                                std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFrac[i].size() << std::endl;
                                gSystem->Exit(EXIT_FAILURE);
                        }
                }

                LauParArray fitFracEffUnCorr = sigDPModel_->getFitFractionsEfficiencyUncorrected();
                if (fitFracEffUnCorr.size() != nSigComp_) {
                        std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr.size() << std::endl;
                        gSystem->Exit(EXIT_FAILURE);
                }
                for (UInt_t i(0); i<nSigComp_; ++i) {
                        if (fitFracEffUnCorr[i].size() != nSigComp_) {
                                std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr[i].size() << std::endl;
                                gSystem->Exit(EXIT_FAILURE);
                        }
                }

                for (UInt_t i(0); i<nSigComp_; ++i) {
                        for (UInt_t j(i); j<nSigComp_; ++j) {
                                fitFrac_[i][j].value(fitFrac[i][j].value());
                                fitFracEffUnCorr_[i][j].value(fitFracEffUnCorr[i][j].value());
                        }
                }

                meanEff_.value(sigDPModel_->getMeanEff().value());
                dpRate_.value(sigDPModel_->getDPRate().value());

                this->clearExtraVarVectors();
                LauParameterList& extraVars = this->extraPars();

                // Then store the final fit parameters, and any extra parameters for
                // the signal model (e.g. fit fractions)
                for (UInt_t i(0); i<nSigComp_; ++i) {
                        for (UInt_t j(i); j<nSigComp_; ++j) {
                                extraVars.push_back(fitFrac_[i][j]);
                                extraVars.push_back(fitFracEffUnCorr_[i][j]);
                        }
                }

                // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions)
                std::vector<LauParameter> extraParams = sigDPModel_->getExtraParameters();
                std::vector<LauParameter>::iterator extraIter;
                for (extraIter = extraParams.begin(); extraIter != extraParams.end(); ++extraIter) {
                        LauParameter extraParameter = (*extraIter);
                        extraVars.push_back(extraParameter);
                }               

                // Now add in the DP efficiency value
                extraVars.push_back(meanEff_);

                // Also add in the DP rate
                extraVars.push_back(dpRate_);

                this->printFitFractions(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 LauSimpleFitModel::printFitFractions(std::ostream& output)
{
        // Print out Fit Fractions, total DP rate and mean efficiency
        for (UInt_t i = 0; i < nSigComp_; i++) {
                output << "FitFraction for component " << i << " (" << coeffPars_[i]->name() << ") = " << fitFrac_[i][i] << std::endl;
        }
        output << "Overall DP rate (integral of matrix element squared) = " << dpRate_ << std::endl;
        output << "Average efficiency weighted by whole DP dynamics = " << meanEff_ << std::endl;
}

void LauSimpleFitModel::writeOutTable(const TString& outputFile)
{
        // Write out the results of the fit to a tex-readable table
        // TODO - need to include the yields in this table
        std::ofstream fout(outputFile);
        LauPrint print;

        std::cout << "INFO in LauSimpleFitModel::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|}" << std::endl;
                fout << "\\hline" << std::endl;
                fout << "Component & FitFraction \\\\" << std::endl;
                fout << "\\hline" << std::endl;
                Double_t fitFracSum(0.0);
                for (UInt_t i = 0; i < nSigComp_; i++) {
                        TString resName = coeffPars_[i]->name();
                        resName = resName.ReplaceAll("_", "\\_");
                        fout << resName << "  &  $";
                        Double_t fitFrac = fitFrac_[i][i].value();
                        fitFracSum += fitFrac;
                        print.printFormat(fout, fitFrac);
                        fout << "$ \\\\" << std::endl;
                }
                fout << "\\hline" << std::endl;

                // Also print out sum of fit fractions
                fout << "Fit Fraction Sum  &  $";
                print.printFormat(fout, fitFracSum);
                fout << "$ \\\\" << std::endl;
                fout << "\\hline" << std::endl;

                fout << "DP rate  &  $";
                print.printFormat(fout, dpRate_.value());
                fout << "$ \\\\" << std::endl;
                fout << "\\hline" << std::endl;

                fout << "$< \\varepsilon >$  &  $";
                print.printFormat(fout, meanEff_.value());
                fout << "$ \\\\" << std::endl;
                fout << "\\hline" << std::endl;
                fout << "\\end{tabular}" << std::endl << std::endl;
        }

        if (!signalPdfs_.empty()) {
                fout << "\\begin{tabular}{|l|c|}" << std::endl;
                fout << "\\hline" << std::endl;
                if (useSCF_ == kTRUE) {
                        fout << "\\Extra TM Signal PDFs' Parameters: & \\\\" << std::endl;
                } else {
                        fout << "\\Extra Signal PDFs' Parameters: & \\\\" << std::endl;
                }
                this->printFitParameters(signalPdfs_, fout);
                if (useSCF_ == kTRUE && !scfPdfs_.empty()) {
                        fout << "\\hline" << std::endl;
                        fout << "\\Extra SCF Signal PDFs' Parameters: & \\\\" << std::endl;
                        this->printFitParameters(scfPdfs_, fout);
                }
                if (usingBkgnd_ == kTRUE && !bkgndPdfs_.empty()) {
                        fout << "\\hline" << std::endl;
                        fout << "\\Extra Background PDFs' Parameters: & \\\\" << std::endl;
                        for (LauBkgndPdfsList::const_iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) {
                                this->printFitParameters(*iter, fout);
                        }
                }
                fout << "\\hline \n\\end{tabular}" << std::endl << std::endl;
        }
}

void LauSimpleFitModel::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) {
                std::cout << "INFO in LauSimpleFitModel::checkInitFitParams : Setting random parameters for the signal DP model" << std::endl;
                this->randomiseInitFitPars();
        }
}

void LauSimpleFitModel::randomiseInitFitPars()
{
        // Only randomise those parameters that are not fixed!
        std::cout << "INFO in LauSimpleFitModel::randomiseInitFitPars : Randomising the initial fit magnitudes and phases of the resonances..." << std::endl;

        for (UInt_t i = 0; i < nSigComp_; i++) {
                coeffPars_[i]->randomiseInitValues();
        }
}

std::pair<LauSimpleFitModel::LauGenInfo,Bool_t> LauSimpleFitModel::eventsToGenerate()
{
        // Determine the number of events to generate for each hypothesis
        // If we're smearing then smear each one individually

        LauGenInfo nEvtsGen;

        // Keep track of whether any yield or asymmetry parameters are blinded
        Bool_t blind = kFALSE;

        // Signal
        Double_t evtWeight(1.0);
        Int_t nEvts = TMath::FloorNint(signalEvents_->genValue());
        if ( nEvts < 0 ) {
                evtWeight = -1.0;
                nEvts = TMath::Abs( nEvts );
        }
        if (this->doPoissonSmearing()) {
                nEvts = LauRandom::randomFun()->Poisson(nEvts);
        }
        nEvtsGen["signal"] = std::make_pair( nEvts, evtWeight );
        if ( signalEvents_->blind() ) {
                blind = kTRUE;
        }

        // Backgrounds
        const UInt_t nBkgnds = this->nBkgndClasses();
        for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
                const TString& bkgndClass = this->bkgndClassName(bkgndID);
                const LauAbsRValue* evtsPar = bkgndEvents_[bkgndID];
                evtWeight = 1.0;
                nEvts = TMath::FloorNint( evtsPar->genValue() );
                if ( nEvts < 0 ) {
                        evtWeight = -1.0;
                        nEvts = TMath::Abs( nEvts );
                }
                if (this->doPoissonSmearing()) {
                        nEvts = LauRandom::randomFun()->Poisson(nEvts);
                }
                nEvtsGen[bkgndClass] = std::make_pair( nEvts, evtWeight );
                if ( evtsPar->blind() ) {
                        blind = kTRUE;
                }
        }

        return std::make_pair( nEvtsGen, blind );
}

Bool_t LauSimpleFitModel::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
        std::pair<LauGenInfo,Bool_t> info = this->eventsToGenerate();
        LauGenInfo nEvts = info.first;
        const Bool_t blind = info.second;

        Bool_t genOK(kTRUE);
        Int_t evtNum(0);

        const UInt_t nBkgnds = this->nBkgndClasses();
        std::vector<TString> bkgndClassNames(nBkgnds);
        std::vector<TString> bkgndClassNamesGen(nBkgnds);
        for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                TString name( this->bkgndClassName(iBkgnd) );
                bkgndClassNames[iBkgnd] = name;
                bkgndClassNamesGen[iBkgnd] = "gen"+name;
        }

        const Bool_t storeSCFTruthInfo = ( useSCF_ || ( this->enableEmbedding() && signalTree_ != 0 && signalTree_->haveBranch("mcMatch") ) );

        // Loop over the hypotheses and generate the requested number of events for each one
        for (LauGenInfo::const_iterator iter = nEvts.begin(); iter != nEvts.end(); ++iter) {

                const TString& type(iter->first);
                Int_t nEvtsGen( iter->second.first );
                Double_t evtWeight( iter->second.second );

                for (Int_t iEvt(0); iEvt<nEvtsGen; ++iEvt) {

                        this->setGenNtupleDoubleBranchValue( "evtWeight", evtWeight );
                        // Add efficiency information
                        this->setGenNtupleDoubleBranchValue( "efficiency", 1 );

                        if (type == "signal") {
                                this->setGenNtupleIntegerBranchValue("genSig",1);
                                for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                                        this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], 0 );
                                }
                                genOK = this->generateSignalEvent();
                                this->setGenNtupleDoubleBranchValue( "efficiency", sigDPModel_->getEvtEff() );
                        } else {
                                this->setGenNtupleIntegerBranchValue("genSig",0);
                                if ( storeSCFTruthInfo ) {
                                        this->setGenNtupleIntegerBranchValue("genTMSig",0);
                                        this->setGenNtupleIntegerBranchValue("genSCFSig",0);
                                }
                                UInt_t bkgndID(0);
                                for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                                        Int_t gen(0);
                                        if ( bkgndClassNames[iBkgnd] == type ) {
                                                gen = 1;
                                                bkgndID = iBkgnd;
                                        }
                                        this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], gen );
                                }
                                genOK = this->generateBkgndEvent(bkgndID);
                        }

                        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->setDPBranchValues();
                        }

                        // Store the event number (within this experiment)
                        this->setGenNtupleIntegerBranchValue("iEvtWithinExpt",evtNum);
                        // and then increment it
                        ++evtNum;

                        this->fillGenNtupleBranches();
                        if ( !blind && (iEvt%500 == 0) ) {
                                std::cout << "INFO in LauSimpleFitModel::genExpt : Generated event number " << iEvt << " out of " << nEvtsGen << " " << type << " events." << std::endl;
                        }
                }

                if (!genOK) {
                        break;
                }
        }

        if (this->useDP() && genOK) {

                sigDPModel_->checkToyMC(kTRUE,kTRUE);

                // Get the fit fractions if they're to be written into the latex table
                if (this->writeLatexTable()) {
                        LauParArray fitFrac = sigDPModel_->getFitFractions();
                        if (fitFrac.size() != nSigComp_) {
                                std::cerr << "ERROR in LauSimpleFitModel::genExpt : Fit Fraction array of unexpected dimension: " << fitFrac.size() << std::endl;
                                gSystem->Exit(EXIT_FAILURE);
                        }
                        for (UInt_t i(0); i<nSigComp_; ++i) {
                                if (fitFrac[i].size() != nSigComp_) {
                                        std::cerr << "ERROR in LauSimpleFitModel::genExpt : Fit Fraction array of unexpected dimension: " << fitFrac[i].size() << std::endl;
                                        gSystem->Exit(EXIT_FAILURE);
                                }
                        }

                        for (UInt_t i(0); i<nSigComp_; ++i) {
                                for (UInt_t j = i; j < nSigComp_; j++) {
                                        fitFrac_[i][j].value(fitFrac[i][j].value());
                                }
                        }
                        meanEff_.value(sigDPModel_->getMeanEff().value());
                        dpRate_.value(sigDPModel_->getDPRate().value());
                }
        }

        // If we're reusing embedded events or if the generation is being
        // reset then clear the lists of used events
        if (reuseSignal_ || !genOK) {
                if (signalTree_) {
                        signalTree_->clearUsedList();
                }
        }
        for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
                LauEmbeddedData* data = bkgndTree_[bkgndID];
                if (reuseBkgnd_[bkgndID] || !genOK) {
                        if (data) {
                                data->clearUsedList();
                        }
                }
        }

        return genOK;
}

Bool_t LauSimpleFitModel::generateSignalEvent()
{
        // Generate signal event
        Bool_t genOK(kTRUE);
        Bool_t genSCF(kFALSE);

        if (this->useDP()) {
                if (this->enableEmbedding() && signalTree_) {
                        if (useReweighting_) {
                                // Select a (random) event from the generated data. Then store the
                                // reconstructed DP co-ords, together with other pdf information,
                                // as the embedded data.
                                genOK = signalTree_->getReweightedEvent(sigDPModel_);
                        } else {

                                signalTree_->getEmbeddedEvent(kinematics_);
                        }
                        if (signalTree_->haveBranch("mcMatch")) {
                                Int_t match = static_cast<Int_t>(signalTree_->getValue("mcMatch"));
                                if (match) {
                                        this->setGenNtupleIntegerBranchValue("genTMSig",1);
                                        this->setGenNtupleIntegerBranchValue("genSCFSig",0);
                                        genSCF = kFALSE;
                                } else {
                                        this->setGenNtupleIntegerBranchValue("genTMSig",0);
                                        this->setGenNtupleIntegerBranchValue("genSCFSig",1);
                                        genSCF = kTRUE;
                                }
                        }
                } else {
                        genOK = sigDPModel_->generate();
                        if ( genOK && useSCF_ ) {
                                Double_t frac(0.0);
                                if ( useSCFHist_ ) {
                                        frac = scfFracHist_->calcEfficiency( kinematics_ );
                                } else {
                                        frac = scfFrac_.genValue();
                                }
                                if ( frac < LauRandom::randomFun()->Rndm() ) {
                                        this->setGenNtupleIntegerBranchValue("genTMSig",1);
                                        this->setGenNtupleIntegerBranchValue("genSCFSig",0);
                                        genSCF = kFALSE;
                                } else {
                                        this->setGenNtupleIntegerBranchValue("genTMSig",0);
                                        this->setGenNtupleIntegerBranchValue("genSCFSig",1);
                                        genSCF = kTRUE;

                                        // Optionally smear the DP position
                                        // of the SCF event
                                        if ( scfMap_ != 0 ) {
                                                Double_t xCoord(0.0), yCoord(0.0);
                                                if ( kinematics_->squareDP() ) {
                                                        xCoord = kinematics_->getmPrime();
                                                        yCoord = kinematics_->getThetaPrime();
                                                } else {
                                                        xCoord = kinematics_->getm13Sq();
                                                        yCoord = kinematics_->getm23Sq();
                                                }

                                                // Find the bin number where this event is generated
                                                Int_t binNo = scfMap_->binNumber( xCoord, yCoord );

                                                // Retrieve the migration histogram
                                                TH2* histo = scfMap_->trueHist( binNo );

                                                const LauAbsEffModel * effModel = sigDPModel_->getEffModel();
                                                do {
                                                        // Get a random point from the histogram
                                                        histo->GetRandom2( xCoord, yCoord );

                                                        // Update the kinematics
                                                        if ( kinematics_->squareDP() ) {
                                                                kinematics_->updateSqDPKinematics( xCoord, yCoord );
                                                        } else {
                                                                kinematics_->updateKinematics( xCoord, yCoord );
                                                        }
                                                } while ( ! effModel->passVeto( kinematics_ ) );
                                        }
                                }
                        }
                }
        } else {
                if (this->enableEmbedding() && signalTree_) {
                        signalTree_->getEmbeddedEvent(0);
                        if (signalTree_->haveBranch("mcMatch")) {
                                Int_t match = static_cast<Int_t>(signalTree_->getValue("mcMatch"));
                                if (match) {
                                        this->setGenNtupleIntegerBranchValue("genTMSig",1);
                                        this->setGenNtupleIntegerBranchValue("genSCFSig",0);
                                        genSCF = kFALSE;
                                } else {
                                        this->setGenNtupleIntegerBranchValue("genTMSig",0);
                                        this->setGenNtupleIntegerBranchValue("genSCFSig",1);
                                        genSCF = kTRUE;
                                }
                        }
                } else if ( useSCF_ ) {
                        Double_t frac = scfFrac_.genValue();
                        if ( frac < LauRandom::randomFun()->Rndm() ) {
                                this->setGenNtupleIntegerBranchValue("genTMSig",1);
                                this->setGenNtupleIntegerBranchValue("genSCFSig",0);
                                genSCF = kFALSE;
                        } else {
                                this->setGenNtupleIntegerBranchValue("genTMSig",0);
                                this->setGenNtupleIntegerBranchValue("genSCFSig",1);
                                genSCF = kTRUE;
                        }
                }
        }
        if (genOK) {
                if ( useSCF_ ) {
                        if ( genSCF ) {
                                this->generateExtraPdfValues(&scfPdfs_, signalTree_);
                        } else {
                                this->generateExtraPdfValues(&signalPdfs_, signalTree_);
                        }
                } else {
                        this->generateExtraPdfValues(&signalPdfs_, signalTree_);
                }
        }
        // Check for problems with the embedding
        if (this->enableEmbedding() && signalTree_ && (signalTree_->nEvents() == signalTree_->nUsedEvents())) {
                std::cerr << "WARNING in LauSimpleFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events." << std::endl;
                signalTree_->clearUsedList();
        }
        return genOK;
}

Bool_t LauSimpleFitModel::generateBkgndEvent(UInt_t bkgndID)
{
        // Generate background event
        Bool_t genOK(kTRUE);

        LauAbsBkgndDPModel* model = bkgndDPModels_[bkgndID];
        LauEmbeddedData* embeddedData(0);
        if (this->enableEmbedding()) {
                embeddedData = bkgndTree_[bkgndID];
        }
        LauPdfList* extraPdfs = &bkgndPdfs_[bkgndID];
        if (this->useDP()) {
                if (embeddedData) {
                        embeddedData->getEmbeddedEvent(kinematics_);
                } else {
                        if (model == 0) {
                                std::cerr << "ERROR in LauSimpleFitModel::generateBkgndEvent : Can't find the DP model for background class \"" << this->bkgndClassName(bkgndID) << "\"." << std::endl;
                                gSystem->Exit(EXIT_FAILURE);
                        }
                        genOK = model->generate();
                }
        } else {
                if (embeddedData) {
                        embeddedData->getEmbeddedEvent(0);
                }
        }
        if (genOK) {
                this->generateExtraPdfValues(extraPdfs, embeddedData);
        }
        // Check for problems with the embedding
        if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) {
                std::cerr << "WARNING in LauSimpleFitModel::generateBkgndEvent : Source of embedded " << this->bkgndClassName(bkgndID) << " events used up, clearing the list of used events." << std::endl;
                embeddedData->clearUsedList();
        }
        return genOK;
}

void LauSimpleFitModel::setupGenNtupleBranches()
{
        // Setup the required ntuple branches
        this->addGenNtupleDoubleBranch("evtWeight");
        this->addGenNtupleIntegerBranch("genSig");
        this->addGenNtupleDoubleBranch("efficiency");
        if ( useSCF_ || ( this->enableEmbedding() && signalTree_ != 0 && signalTree_->haveBranch("mcMatch") ) ) {
                this->addGenNtupleIntegerBranch("genTMSig");
                this->addGenNtupleIntegerBranch("genSCFSig");
        }
        const UInt_t nBkgnds = this->nBkgndClasses();
        for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                TString name( this->bkgndClassName(iBkgnd) );
                name.Prepend("gen");
                this->addGenNtupleIntegerBranch(name);
        }
        if (this->useDP() == kTRUE) {
                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 (kinematics_->squareDP()) {
                        this->addGenNtupleDoubleBranch("mPrime");
                        this->addGenNtupleDoubleBranch("thPrime");
                }
        }
        for (LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter) {
                std::vector<TString> varNames = (*pdf_iter)->varNames();
                for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
                        if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
                                this->addGenNtupleDoubleBranch( (*var_iter) );
                        }
                }
        }
}

void LauSimpleFitModel::setDPBranchValues()
{
        // 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());
        }
}

void LauSimpleFitModel::generateExtraPdfValues(LauPdfList* extraPdfs, LauEmbeddedData* embeddedData)
{
        // 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 LauSimpleFitModel::propagateParUpdates()
{
        // Update the total normalisation for the signal likelihood
        if (this->useDP() == kTRUE) {
                this->updateCoeffs();
                sigDPModel_->updateCoeffs(coeffs_);
        }

        // Update the signal events from the background events if not doing an extended fit
        if ( !this->doEMLFit() && !signalEvents_->fixed() ) {
                this->updateSigEvents();
        }
}

void LauSimpleFitModel::updateSigEvents()
{
        // The background parameters will have been set from Minuit.
        // We need to update the signal events using these.
        Double_t nTotEvts = this->eventsPerExpt();

        signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts);
        for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
                LauAbsRValue* nBkgndEvents = (*iter);
                if ( nBkgndEvents->isLValue() ) {
                        LauParameter* yield = dynamic_cast<LauParameter*>( nBkgndEvents );
                        yield->range(-2.0*nTotEvts,2.0*nTotEvts);
                }
        }

        if ( signalEvents_->fixed() ) {
                return;
        }

        // Subtract background events (if any) from signal.
        Double_t signalEvents = nTotEvts;
        if ( usingBkgnd_ == kTRUE ) {
                for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
                        signalEvents -= (*iter)->value();
                }
        }
        signalEvents_->value(signalEvents);
}

void LauSimpleFitModel::cacheInputFitVars()
{
        // Fill the internal data trees of the signal and backgrond models.

        LauFitDataTree* inputFitData = this->fitData();

        // First the Dalitz plot variables (m_ij^2)
        if (this->useDP() == kTRUE) {

                // need to append SCF smearing bins before caching DP amplitudes
                if ( scfMap_ != 0 ) {
                        this->appendBinCentres( inputFitData );
                }
                sigDPModel_->fillDataTree(*inputFitData);

                if (usingBkgnd_ == kTRUE) {
                        for (LauBkgndDPModelList::iterator iter = bkgndDPModels_.begin(); iter != bkgndDPModels_.end(); ++iter) {
                                (*iter)->fillDataTree(*inputFitData);
                        }
                }
        }

        // ...and then the extra PDFs
        this->cacheInfo(signalPdfs_, *inputFitData);
        this->cacheInfo(scfPdfs_, *inputFitData);
        for (LauBkgndPdfsList::iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) {
                this->cacheInfo((*iter), *inputFitData);
        }

        // and finally the SCF fractions and jacobians
        if ( useSCF_ && useSCFHist_ ) {
                if ( !inputFitData->haveBranch( "m13Sq" )  || !inputFitData->haveBranch( "m23Sq" ) ) {
                        std::cerr << "ERROR in LauSimpleFitModel::cacheInputFitVars : Input data does not contain DP branches and so can't cache the SCF fraction." << std::endl;
                        gSystem->Exit(EXIT_FAILURE);
                }

                UInt_t nEvents = inputFitData->nEvents();
                recoSCFFracs_.clear();
                recoSCFFracs_.reserve( nEvents );
                if ( kinematics_->squareDP() ) {
                        recoJacobians_.clear();
                        recoJacobians_.reserve( nEvents );
                }
                for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
                        const LauFitData& dataValues = inputFitData->getData(iEvt);
                        LauFitData::const_iterator m13_iter = dataValues.find("m13Sq");
                        LauFitData::const_iterator m23_iter = dataValues.find("m23Sq");
                        kinematics_->updateKinematics( m13_iter->second, m23_iter->second );
                        Double_t scfFrac = scfFracHist_->calcEfficiency( kinematics_ );
                        recoSCFFracs_.push_back( scfFrac );
                        if ( kinematics_->squareDP() ) {
                                recoJacobians_.push_back( kinematics_->calcSqDPJacobian() );
                        }
                }
        }
}

void LauSimpleFitModel::appendBinCentres( LauFitDataTree* inputData )
{
        // We'll be caching the DP amplitudes and efficiencies of the centres of the true bins.
        // To do so, we attach some fake points at the end of inputData, the number of the entry
        // minus the total number of events corresponding to the number of the histogram for that
        // given true bin in the LauScfMap object. (What this means is that when Laura is provided with
        // the LauScfMap object by the user, it's the latter who has to make sure that it contains the
        // right number of histograms and in exactly the right order!) 

        // Get the x and y co-ordinates of the bin centres
        std::vector<Double_t> binCentresXCoords;
        std::vector<Double_t> binCentresYCoords;
        scfMap_->listBinCentres(binCentresXCoords, binCentresYCoords);

        // The SCF histograms could be in square Dalitz plot histograms.
        // The dynamics takes normal Dalitz plot coords, so we might have to convert them back.
        Bool_t sqDP = kinematics_->squareDP();
        UInt_t nBins = binCentresXCoords.size();
        fakeSCFFracs_.clear();
        fakeSCFFracs_.reserve( nBins );
        if ( sqDP ) {
                fakeJacobians_.clear();
                fakeJacobians_.reserve( nBins );
        }

        for (UInt_t iBin = 0; iBin < nBins; ++iBin) {

                if ( sqDP ) {
                        kinematics_->updateSqDPKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]);
                        binCentresXCoords[iBin] = kinematics_->getm13Sq();
                        binCentresYCoords[iBin] = kinematics_->getm23Sq();
                        fakeJacobians_.push_back( kinematics_->calcSqDPJacobian() );
                } else {
                        kinematics_->updateKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]);
                }

                fakeSCFFracs_.push_back( scfFracHist_->calcEfficiency( kinematics_ ) );
        }

        // Set up inputFitVars_ object to hold the fake events
        inputData->appendFakePoints(binCentresXCoords,binCentresYCoords);
}

Double_t LauSimpleFitModel::getTotEvtLikelihood(UInt_t iEvt)
{
        // Get the DP likelihood for signal and backgrounds
        this->getEvtDPLikelihood(iEvt);

        // Get the combined extra PDFs likelihood for signal and backgrounds
        this->getEvtExtraLikelihoods(iEvt);

        // If appropriate, combine the TM and SCF likelihoods
        Double_t sigLike = sigDPLike_ * sigExtraLike_;
        if ( useSCF_ ) {
                Double_t scfFrac(0.0);
                if (useSCFHist_) {
                        scfFrac = recoSCFFracs_[iEvt];
                } else {
                        scfFrac = scfFrac_.unblindValue();
                }
                sigLike *= (1.0 - scfFrac);
                if ( (scfMap_ != 0) && (this->useDP() == kTRUE) ) {
                        // if we're smearing the SCF DP PDF then the SCF frac
                        // is already included in the SCF DP likelihood
                        sigLike += (scfDPLike_ * scfExtraLike_);
                } else {
                        sigLike += (scfFrac * scfDPLike_ * scfExtraLike_);
                }
        }

        // Construct the total event likelihood
        Double_t likelihood = signalEvents_->unblindValue() * sigLike;
        const UInt_t nBkgnds = this->nBkgndClasses();
        for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
                likelihood += (bkgndEvents_[bkgndID]->unblindValue() * bkgndDPLike_[bkgndID] * bkgndExtraLike_[bkgndID]);
        }

        return likelihood;
}

Double_t LauSimpleFitModel::getEventSum() const
{
        Double_t eventSum(0.0);
        eventSum += signalEvents_->unblindValue();
        for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) {
                eventSum += (*iter)->unblindValue();
        }
        return eventSum;
}

void LauSimpleFitModel::getEvtDPLikelihood(UInt_t iEvt)
{
        // Function to return the signal and background likelihoods for the
        // Dalitz plot for the given event evtNo.

        if (this->useDP() == kTRUE) {
                sigDPModel_->calcLikelihoodInfo(iEvt);
                sigDPLike_ = sigDPModel_->getEvtLikelihood();

                if ( useSCF_ == kTRUE ) {
                        if ( scfMap_ == 0 ) {
                                // we're not smearing the SCF DP position
                                // so the likelihood is the same as the TM
                                scfDPLike_ = sigDPLike_;
                        } else {
                                // calculate the smeared SCF DP likelihood
                                scfDPLike_ = this->getEvtSCFDPLikelihood(iEvt);
                        }
                }

                const UInt_t nBkgnds = this->nBkgndClasses();
                for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
                        if (usingBkgnd_ == kTRUE) {
                                bkgndDPLike_[bkgndID] = bkgndDPModels_[bkgndID]->getLikelihood(iEvt);
                        } else {
                                bkgndDPLike_[bkgndID] = 0.0;
                        }
                }
        } else {
                // There's always going to be a term in the likelihood for the
                // signal, so we'd better not zero it.
                sigDPLike_ = 1.0;
                scfDPLike_ = 1.0;

                const UInt_t nBkgnds = this->nBkgndClasses();
                for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
                        if (usingBkgnd_ == kTRUE) {
                                bkgndDPLike_[bkgndID] = 1.0;
                        } else {
                                bkgndDPLike_[bkgndID] = 0.0;
                        }
                }
        }
}

Double_t LauSimpleFitModel::getEvtSCFDPLikelihood(UInt_t iEvt)
{
        Double_t scfDPLike(0.0);

        Double_t recoJacobian(1.0);
        Double_t xCoord(0.0);
        Double_t yCoord(0.0);
        Bool_t squareDP = kinematics_->squareDP();
        if ( squareDP ) {
                xCoord = sigDPModel_->getEvtmPrime();
                yCoord = sigDPModel_->getEvtthPrime();
                recoJacobian = recoJacobians_[iEvt];
        } else {
                xCoord = sigDPModel_->getEvtm13Sq();
                yCoord = sigDPModel_->getEvtm23Sq();
        }

        // Find the bin that our reco event falls in
        Int_t recoBin = scfMap_->binNumber( xCoord, yCoord );

        // Find out which true Bins contribute to the given reco bin
        const std::vector<Int_t>* trueBins = scfMap_->trueBins(recoBin);

        Int_t nDataEvents = this->eventsPerExpt();

        // Loop over the true bins
        for (std::vector<Int_t>::const_iterator iter = trueBins->begin(); iter != trueBins->end(); ++iter)
        {
                Int_t trueBin = (*iter);

                // prob of a true event in the given true bin migrating to the reco bin
                Double_t pRecoGivenTrue =  scfMap_->prob( recoBin, trueBin );

                // We've cached the DP amplitudes and the efficiency for the
                // true bin centres, just after the data points
                sigDPModel_->calcLikelihoodInfo( nDataEvents + trueBin );
                Double_t pTrue = sigDPModel_->getEvtLikelihood();

                // Get the cached SCF fraction (and jacobian if we're using the square DP)
                Double_t scfFraction = fakeSCFFracs_[ trueBin ];
                Double_t jacobian(1.0);
                if ( squareDP ) {
                        jacobian = fakeJacobians_[ trueBin ];
                }

                scfDPLike += pTrue * jacobian * scfFraction * pRecoGivenTrue;

        }

        // Divide by the reco jacobian
        scfDPLike /= recoJacobian;

        return scfDPLike;
}

void LauSimpleFitModel::getEvtExtraLikelihoods(UInt_t iEvt)
{
        // Function to return the signal and background likelihoods for the
        // extra variables for the given event evtNo.

        sigExtraLike_ = 1.0;

        for (LauPdfList::iterator iter = signalPdfs_.begin(); iter != signalPdfs_.end(); ++iter) {
                (*iter)->calcLikelihoodInfo(iEvt);
                sigExtraLike_ *= (*iter)->getLikelihood();
        }
        if (useSCF_) {
                scfExtraLike_ = 1.0;
                for (LauPdfList::iterator iter = scfPdfs_.begin(); iter != scfPdfs_.end(); ++iter) {
                        (*iter)->calcLikelihoodInfo(iEvt);
                        scfExtraLike_ *= (*iter)->getLikelihood();
                }
        }

        const UInt_t nBkgnds = this->nBkgndClasses();
        for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) {
                if (usingBkgnd_) {
                        bkgndExtraLike_[bkgndID] = 1.0;
                        LauPdfList& pdfList = bkgndPdfs_[bkgndID];
                        for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
                                (*pdf_iter)->calcLikelihoodInfo(iEvt);
                                bkgndExtraLike_[bkgndID] *= (*pdf_iter)->getLikelihood();
                        }
                } else {
                        bkgndExtraLike_[bkgndID] = 0.0;
                }
        }
}

void LauSimpleFitModel::updateCoeffs() 
{
        coeffs_.clear();
        coeffs_.reserve(nSigComp_);
        for (UInt_t i = 0; i < nSigComp_; i++) {
                coeffs_.push_back(coeffPars_[i]->particleCoeff());
        }
}

void LauSimpleFitModel::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");
                if ( sigDPModel_->usingScfModel() ) {
                        this->addSPlotNtupleDoubleBranch("scffraction");
                }
        }

        // Store the total event likelihood for each species.
        if (useSCF_) {
                this->addSPlotNtupleDoubleBranch("sigTMTotalLike");
                this->addSPlotNtupleDoubleBranch("sigSCFTotalLike");
                this->addSPlotNtupleDoubleBranch("sigSCFFrac");
        } else {
                this->addSPlotNtupleDoubleBranch("sigTotalLike");
        }
        if (usingBkgnd_) {
                const UInt_t nBkgnds = this->nBkgndClasses();
                for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                        TString name( this->bkgndClassName(iBkgnd) );
                        name += "TotalLike";
                        this->addSPlotNtupleDoubleBranch(name);
                }
        }

        // Store the DP likelihoods
        if (this->useDP()) {
                if (useSCF_) {
                        this->addSPlotNtupleDoubleBranch("sigTMDPLike");
                        this->addSPlotNtupleDoubleBranch("sigSCFDPLike");
                } else {
                        this->addSPlotNtupleDoubleBranch("sigDPLike");
                }
                if (usingBkgnd_) {
                        const UInt_t nBkgnds = this->nBkgndClasses();
                        for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                                TString name( this->bkgndClassName(iBkgnd) );
                                name += "DPLike";
                                this->addSPlotNtupleDoubleBranch(name);
                        }
                }
        }

        // Store the likelihoods for each extra PDF
        if (useSCF_) {
                this->addSPlotNtupleBranches(&signalPdfs_, "sigTM");
                this->addSPlotNtupleBranches(&scfPdfs_, "sigSCF");
        } else {
                this->addSPlotNtupleBranches(&signalPdfs_, "sig");
        }
        if (usingBkgnd_) {
                const UInt_t nBkgnds = this->nBkgndClasses();
                for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                        const TString& bkgndClass = this->bkgndClassName(iBkgnd);
                        const LauPdfList* pdfList = &(bkgndPdfs_[iBkgnd]);
                        this->addSPlotNtupleBranches(pdfList, bkgndClass);
                }
        }
}

void LauSimpleFitModel::addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix)
{
        if (extraPdfs) {
                // 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);
                        std::vector<TString> varNames = (*pdf_iter)->varNames();
                        for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_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 = varNames.begin(); var_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 LauSimpleFitModel::addSPlotNtupleBranches : Can't yet deal with 3D PDFs." << std::endl;
                        }
                }
        }
}

Double_t LauSimpleFitModel::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) {
                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);
                        std::vector<TString> varNames = (*pdf_iter)->varNames();
                        for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_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 = varNames.begin(); var_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 LauSimpleFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs." << std::endl;
                        }
                }
        }
        return totalLike;
}

LauSPlot::NameSet LauSimpleFitModel::variableNames() const
{
        LauSPlot::NameSet nameSet;
        if (this->useDP()) {
                nameSet.insert("DP");
        }
        // Loop through all the signal PDFs
        for (LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter) {
                // Loop over the variables involved in each PDF
                std::vector<TString> varNames = (*pdf_iter)->varNames();
                for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_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 LauSimpleFitModel::freeSpeciesNames() const
{
        LauSPlot::NumbMap numbMap;
        if (!signalEvents_->fixed() && this->doEMLFit()) {
                numbMap["sig"] = signalEvents_->genValue();
        }
        if ( usingBkgnd_ ) {
                const UInt_t nBkgnds = this->nBkgndClasses();
                for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                        const TString& bkgndClass = this->bkgndClassName(iBkgnd);
                        const LauAbsRValue* par = bkgndEvents_[iBkgnd];
                        if (!par->fixed()) {
                                numbMap[bkgndClass] = par->genValue();
                                if ( ! par->isLValue() ) {
                                        std::cerr << "WARNING in LauSimpleFitModel::freeSpeciesNames : \"" << par->name() << "\" is a LauFormulaPar, which implies it is perhaps not entirely free to float in the fit, so the sWeight calculation may not be reliable" << std::endl;
                                }
                        }
                }
        }
        return numbMap;
}

LauSPlot::NumbMap LauSimpleFitModel::fixdSpeciesNames() const
{
        LauSPlot::NumbMap numbMap;
        if (signalEvents_->fixed() && this->doEMLFit()) {
                numbMap["sig"] = signalEvents_->genValue();
        }
        if ( usingBkgnd_ ) {
                const UInt_t nBkgnds = this->nBkgndClasses();
                for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                        const TString& bkgndClass = this->bkgndClassName(iBkgnd);
                        const LauAbsRValue* par = bkgndEvents_[iBkgnd];
                        if (par->fixed()) {
                                numbMap[bkgndClass] = par->genValue();
                        }
                }
        }
        return numbMap;
}

LauSPlot::TwoDMap LauSimpleFitModel::twodimPDFs() const
{
        LauSPlot::TwoDMap twodimMap;

        for (LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter) {
                // Count the number of input variables that are not DP variables
                UInt_t nVars(0);
                std::vector<TString> varNames = (*pdf_iter)->varNames();
                for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
                        if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
                                ++nVars;
                        }
                }
                if ( nVars == 2 ) {
                        if (useSCF_) {
                                twodimMap.insert( std::make_pair( "sigTM", std::make_pair( varNames[0], varNames[1] ) ) );
                        } else {
                                twodimMap.insert( std::make_pair( "sig", std::make_pair( varNames[0], varNames[1] ) ) );
                        }
                }
        }

        if ( useSCF_ ) {
                for (LauPdfList::const_iterator pdf_iter = scfPdfs_.begin(); pdf_iter != scfPdfs_.end(); ++pdf_iter) {
                        // Count the number of input variables that are not DP variables
                        UInt_t nVars(0);
                        std::vector<TString> varNames = (*pdf_iter)->varNames();
                        for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
                                if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
                                        ++nVars;
                                }
                        }
                        if ( nVars == 2 ) {
                                twodimMap.insert( std::make_pair( "sigSCF", std::make_pair( varNames[0], varNames[1] ) ) );
                        }
                }
        }

        if (usingBkgnd_) {
                const UInt_t nBkgnds = this->nBkgndClasses();
                for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                        const TString& bkgndClass = this->bkgndClassName(iBkgnd);
                        const LauPdfList& pdfList = bkgndPdfs_[iBkgnd];
                        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);
                                std::vector<TString> varNames = (*pdf_iter)->varNames();
                                for ( std::vector<TString>::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) {
                                        if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) {
                                                ++nVars;
                                        }
                                }
                                if ( nVars == 2 ) {
                                        twodimMap.insert( std::make_pair( bkgndClass, std::make_pair( varNames[0], varNames[1] ) ) );
                                }
                        }
                }
        }

        return twodimMap;
}

void LauSimpleFitModel::storePerEvtLlhds()
{
        std::cout << "INFO in LauSimpleFitModel::storePerEvtLlhds : Storing per-event likelihood values..." << std::endl;

        // 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

        LauFitDataTree* inputFitData = this->fitData();

        if (!this->useDP() && this->storeDPEff()) {
                sigDPModel_->initialise(coeffs_);
                sigDPModel_->fillDataTree(*inputFitData);
        }

        UInt_t evtsPerExpt(this->eventsPerExpt());

        for (UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt) {

                this->setSPlotNtupleIntegerBranchValue("iExpt",this->iExpt());
                this->setSPlotNtupleIntegerBranchValue("iEvtWithinExpt",iEvt);

                // the DP information
                this->getEvtDPLikelihood(iEvt);
                if (this->storeDPEff()) {
                        if (!this->useDP()) {
                                sigDPModel_->calcLikelihoodInfo(iEvt);
                        }
                        this->setSPlotNtupleDoubleBranchValue("efficiency",sigDPModel_->getEvtEff());
                        if ( sigDPModel_->usingScfModel() ) {
                                this->setSPlotNtupleDoubleBranchValue("scffraction",sigDPModel_->getEvtScfFraction());
                        }
                }
                if (this->useDP()) {
                        sigTotalLike_ = sigDPLike_;
                        if (useSCF_) {
                                scfTotalLike_ = scfDPLike_;
                                this->setSPlotNtupleDoubleBranchValue("sigTMDPLike",sigDPLike_);
                                this->setSPlotNtupleDoubleBranchValue("sigSCFDPLike",scfDPLike_);
                        } else {
                                this->setSPlotNtupleDoubleBranchValue("sigDPLike",sigDPLike_);
                        }
                        if (usingBkgnd_) {
                                const UInt_t nBkgnds = this->nBkgndClasses();
                                for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                                        TString name = this->bkgndClassName(iBkgnd);
                                        name += "DPLike";
                                        this->setSPlotNtupleDoubleBranchValue(name,bkgndDPLike_[iBkgnd]);
                                }
                        }
                } else {
                        sigTotalLike_ = 1.0;
                        if (useSCF_) {
                                scfTotalLike_ = 1.0;
                        }
                        if (usingBkgnd_) {
                                const UInt_t nBkgnds = this->nBkgndClasses();
                                for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                                        bkgndTotalLike_[iBkgnd] = 1.0;
                                }
                        }
                }

                // the signal PDF values
                if ( useSCF_ ) {
                        sigTotalLike_ *= this->setSPlotNtupleBranchValues(&signalPdfs_, "sigTM", iEvt);
                        scfTotalLike_ *= this->setSPlotNtupleBranchValues(&scfPdfs_, "sigSCF", iEvt);
                } else {
                        sigTotalLike_ *= this->setSPlotNtupleBranchValues(&signalPdfs_, "sig", iEvt);
                }

                // the background PDF values
                if (usingBkgnd_) {
                        const UInt_t nBkgnds = this->nBkgndClasses();
                        for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                                const TString& bkgndClass = this->bkgndClassName(iBkgnd);
                                LauPdfList& pdfs = bkgndPdfs_[iBkgnd];
                                bkgndTotalLike_[iBkgnd] *= this->setSPlotNtupleBranchValues(&(pdfs), bkgndClass, iEvt);
                        }
                }

                // the total likelihoods
                if (useSCF_) {
                        Double_t scfFrac(0.0);
                        if ( useSCFHist_ ) {
                                scfFrac = recoSCFFracs_[iEvt];
                        } else {
                                scfFrac = scfFrac_.unblindValue();
                        }
                        this->setSPlotNtupleDoubleBranchValue("sigSCFFrac",scfFrac);
                        sigTotalLike_ *= ( 1.0 - scfFrac );
                        if ( scfMap_ == 0 ) {
                                scfTotalLike_ *= scfFrac;
                        }
                        this->setSPlotNtupleDoubleBranchValue("sigTMTotalLike",sigTotalLike_);
                        this->setSPlotNtupleDoubleBranchValue("sigSCFTotalLike",scfTotalLike_);
                } else {
                        this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_);
                }
                if (usingBkgnd_) {
                        const UInt_t nBkgnds = this->nBkgndClasses();
                        for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) {
                                TString name = this->bkgndClassName(iBkgnd);
                                name += "TotalLike";
                                this->setSPlotNtupleDoubleBranchValue(name,bkgndTotalLike_[iBkgnd]);
                        }
                }
                // fill the tree
                this->fillSPlotNtupleBranches();
        }
        std::cout << "INFO in LauSimpleFitModel::storePerEvtLlhds : Finished storing per-event likelihood values." << std::endl;
}

void LauSimpleFitModel::embedSignal(const TString& fileName, const TString& treeName,
                Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment,
                Bool_t useReweighting)
{
        if (signalTree_) {
                std::cerr << "ERROR in LauSimpleFitModel::embedSignal : Already embedding signal from a file." << std::endl;
                return;
        }

        signalTree_ = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);

        Bool_t dataOK = signalTree_->findBranches();
        if (!dataOK) {
                delete signalTree_; signalTree_ = 0;
                std::cerr << "ERROR in LauSimpleFitModel::embedSignal : Problem creating data tree for embedding." << std::endl;
                return;
        }

        reuseSignal_ = reuseEventsWithinEnsemble;
        useReweighting_ = useReweighting;

        if (this->enableEmbedding() == kFALSE) {
                this->enableEmbedding(kTRUE);
        }
}

void LauSimpleFitModel::embedBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName,
                Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment)
{
        if ( ! this->validBkgndClass( bkgndClass ) ) {
                std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl;
                std::cerr << "                                       : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
                return;
        }

        UInt_t bkgndID = this->bkgndClassID( bkgndClass );

        if (bkgndTree_[bkgndID]) {
                std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Already embedding background from a file." << std::endl;
                return;
        }

        bkgndTree_[bkgndID] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment);

        Bool_t dataOK = bkgndTree_[bkgndID]->findBranches();
        if (!dataOK) {
                delete bkgndTree_[bkgndID]; bkgndTree_[bkgndID] = 0;
                std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Problem creating data tree for embedding." << std::endl;
                return;
        }

        reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble;

        if (this->enableEmbedding() == kFALSE) {
                this->enableEmbedding(kTRUE);
        }
}

void LauSimpleFitModel::weightEvents( const TString& dataFileName, const TString& dataTreeName )
{
        // Routine to provide weights for events that are uniformly distributed
        // in the DP (or square DP) so as to reproduce the given DP model

        if ( kinematics_->squareDP() ) {
                std::cout << "INFO in LauSimpleFitModel::weightEvents : will create weights assuming events were generated flat in the square DP" << std::endl;
        } else {
                std::cout << "INFO in LauSimpleFitModel::weightEvents : will create weights assuming events were generated flat in phase space" << std::endl;
        }

        // This reads in the given dataFile and creates an input
        // fit data tree that stores them for all events and experiments.
        Bool_t dataOK = this->verifyFitData(dataFileName,dataTreeName);
        if (!dataOK) {
                std::cerr << "ERROR in LauSimpleFitModel::weightEvents : Problem caching the data." << std::endl;
                return;
        }

        LauFitDataTree* inputFitData = this->fitData();

        if ( ! inputFitData->haveBranch( "m13Sq_MC" ) || ! inputFitData->haveBranch( "m23Sq_MC" ) ) {
                std::cerr << "WARNING in LauSimpleFitModel::weightEvents : Cannot find MC truth DP coordinate branches in supplied data, aborting." << std::endl;
                return;
        }

        // Create the ntuple to hold the DP weights
        TString weightsFileName( dataFileName );
        Ssiz_t index = weightsFileName.Last('.');
        weightsFileName.Insert( index, "_DPweights" );
        LauGenNtuple * weightsTuple = new LauGenNtuple( weightsFileName, dataTreeName );
        weightsTuple->addIntegerBranch("iExpt");
        weightsTuple->addIntegerBranch("iEvtWithinExpt");
        weightsTuple->addDoubleBranch("dpModelWeight");

        UInt_t nExpmt = this->nExpt();
        UInt_t firstExpmt = this->firstExpt();
        for (UInt_t iExpmt = firstExpmt; iExpmt < (firstExpmt+nExpmt); ++iExpmt) {

                inputFitData->readExperimentData(iExpmt);
                UInt_t nEvents = inputFitData->nEvents();

                if (nEvents < 1) {
                        std::cerr << "WARNING in LauSimpleFitModel::weightEvents : Zero events in experiment " << iExpmt << ", skipping..." << std::endl;
                        continue;
                }

                weightsTuple->setIntegerBranchValue( "iExpt", iExpmt );

                // Calculate and store the weights for the events in this experiment
                for ( UInt_t iEvent(0); iEvent < nEvents; ++iEvent ) {

                        weightsTuple->setIntegerBranchValue( "iEvtWithinExpt", iEvent );

                        const LauFitData& evtData = inputFitData->getData( iEvent );

                        Double_t m13Sq_MC = evtData.find("m13Sq_MC")->second;
                        Double_t m23Sq_MC = evtData.find("m23Sq_MC")->second;

                        Double_t dpModelWeight(0.0);

                        if ( kinematics_->withinDPLimits( m13Sq_MC, m23Sq_MC ) ) {

                                kinematics_->updateKinematics( m13Sq_MC, m23Sq_MC );
                                dpModelWeight = sigDPModel_->getEventWeight();

                                if ( kinematics_->squareDP() ) {
                                        dpModelWeight *= kinematics_->calcSqDPJacobian();
                                }

                                dpModelWeight /= sigDPModel_->getDPNorm();
                        }

                        weightsTuple->setDoubleBranchValue( "dpModelWeight", dpModelWeight );
                        weightsTuple->fillBranches();
                }

        }

        weightsTuple->buildIndex( "iExpt", "iEvtWithinExpt" );
        weightsTuple->addFriendTree(dataFileName, dataTreeName);
        weightsTuple->writeOutGenResults();

        delete weightsTuple;
}


void LauSimpleFitModel::savePDFPlots(const TString& label)
{
   savePDFPlotsWave(label, 0);
   savePDFPlotsWave(label, 1);
   savePDFPlotsWave(label, 2);

   std::cout << "LauCPFitModel::plot" << std::endl;
//      ((LauIsobarDynamics*)sigDPModel_)->plot();

   

   //Double_t minm13 = sigDPModel_->getKinematics()->getm13Min();
   Double_t minm13 = 0.0;
   Double_t maxm13 = sigDPModel_->getKinematics()->getm13Max();
   //Double_t minm23 = sigDPModel_->getKinematics()->getm23Min();
   Double_t minm23 = 0.0;
   Double_t maxm23 = sigDPModel_->getKinematics()->getm23Max();

   Double_t mins13 = minm13*minm13;
   Double_t maxs13 = maxm13*maxm13;
   Double_t mins23 = minm23*minm23;
   Double_t maxs23 = maxm23*maxm23;

   Double_t s13, s23, chPdf;

   Int_t n13=200.00, n23=200.00;
   Double_t delta13, delta23;
   delta13 = (maxs13 - mins13)/n13;
   delta23 = (maxs23 - mins23)/n23;
   UInt_t nAmp = sigDPModel_->getnCohAmp();
   for (UInt_t resID = 0; resID <= nAmp; ++resID)
   {
        TGraph2D *dt = new TGraph2D();
        TString resName = "TotalAmp";
        if (resID != nAmp){
                TString tStrResID = Form("%d", resID);
                const LauIsobarDynamics* model = sigDPModel_;
                const LauAbsResonance* resonance = model->getResonance(resID);
                resName = resonance->getResonanceName();
                std::cout << "resName = " << resName << std::endl;
        }


        resName.ReplaceAll("(", "");
        resName.ReplaceAll(")", "");
        resName.ReplaceAll("*", "Star");

        TCanvas *c = new TCanvas("c"+resName+label,resName+" ("+label+")",0,0,600,400);
        dt->SetName(resName+label);
        dt->SetTitle(resName+" ("+label+")");
        Int_t count=0;
        for (Int_t i=0; i<n13; i++) {
                s13 = mins13 + i*delta13;
                for (Int_t j=0; j<n23; j++) {
                        s23 = mins23 + j*delta23;
                        if (sigDPModel_->getKinematics()->withinDPLimits2(s23, s13))
                        {
                                //if (s13 > 4) continue;
                                sigDPModel_->calcLikelihoodInfo(s13, s23);
                                LauComplex chAmp = sigDPModel_->getEvtDPAmp();
                                if (resID != nAmp){
                                        chAmp = sigDPModel_->getFullAmplitude(resID);
                                }
                                chPdf = chAmp.abs2();
                                //if ((z > 0.04)||(z < -0.03)) continue;
                                //z = TMath::Log(z);
                                if (sigDPModel_->getDaughters()->gotSymmetricalDP()){
                                        Double_t sLow = s13;
                                        Double_t sHigh = s23;
                                        if (sLow>sHigh) {
                                                continue;
                                                //sLow = s23;
                                                //sHigh = s13;
                                        }

                                        //if (sLow > 3.5) continue;
                                        //if (i<10) std::cout << "SymmetricalDP" << std::endl;
                                        //if (z>0.02) z = 0.02;
                                        //if (z<-0.02) z = -0.02;

                                        dt->SetPoint(count,sHigh,sLow,chPdf);
                                        count++;
                                }
                                else {
                                        dt->SetPoint(count,s13,s23,chPdf);
                                        count++;

                                }

                        }
                }
        }
        gStyle->SetPalette(1);
        dt->GetXaxis()->SetTitle("m_{K#pi}(low)");
        dt->GetYaxis()->SetTitle("m_{K#pi}(high)");
        dt->Draw("SURF1");
        dt->GetXaxis()->SetTitle("m_{K#pi}(low)");
        dt->GetYaxis()->SetTitle("m_{K#pi}(high)");
        c->SaveAs("plot_2D_"+resName + "_"+label+".C");

   }
}

void LauSimpleFitModel::savePDFPlotsWave(const TString& label, const Int_t& spin)
{

        std::cout << "label = "<< label <<  ", spin = "<< spin << std::endl;

        TString tStrResID = "S_Wave";
        if (spin == 1) tStrResID = "P_Wave";
        if (spin == 2) tStrResID = "D_Wave";

        std::cout << "LauSimpleFitModel::savePDFPlotsWave: "<< tStrResID << std::endl;

        TCanvas *c = new TCanvas("c"+tStrResID+label,tStrResID+" ("+label+")",0,0,600,400);

        Double_t minm13 = 0.0;
        Double_t maxm13 = sigDPModel_->getKinematics()->getm13Max();
        Double_t minm23 = 0.0;
        Double_t maxm23 = sigDPModel_->getKinematics()->getm23Max();

        Double_t mins13 = minm13*minm13;
        Double_t maxs13 = maxm13*maxm13;
        Double_t mins23 = minm23*minm23;
        Double_t maxs23 = maxm23*maxm23;

        Double_t s13, s23, chPdf;
        TGraph2D *dt = new TGraph2D();
        dt->SetName(tStrResID+label);
        dt->SetTitle(tStrResID+" ("+label+")");

        Int_t n13=200.00, n23=200.00;
        Double_t delta13, delta23;
        delta13 = (maxs13 - mins13)/n13;
        delta23 = (maxs23 - mins23)/n23;
        UInt_t nAmp = sigDPModel_->getnCohAmp();

        Int_t count=0;
        for (Int_t i=0; i<n13; i++) {
                s13 = mins13 + i*delta13;
                for (Int_t j=0; j<n23; j++) {
                        s23 = mins23 + j*delta23;
                        if (sigDPModel_->getKinematics()->withinDPLimits2(s23, s13))
                        {
                                //if (s13 > 4) continue;
                                LauComplex chAmp(0,0);
                                Bool_t noWaveRes = kTRUE;
                                sigDPModel_->calcLikelihoodInfo(s13, s23);
                                for (UInt_t resID = 0; resID < nAmp; ++resID)
                                {
                                        const LauIsobarDynamics* model = sigDPModel_;
                                        const LauAbsResonance* resonance = model->getResonance(resID);
                                        Int_t spin_res = resonance->getSpin();
                                        if (spin != spin_res) continue;
                                        noWaveRes = kFALSE;
                                        chAmp += sigDPModel_->getFullAmplitude(resID);
                                }

                                if (noWaveRes) return;

                                chPdf = chAmp.abs2();
                                if (sigDPModel_->getDaughters()->gotSymmetricalDP()){
                                        Double_t sLow = s13;
                                        Double_t sHigh = s23;
                                        if (sLow>sHigh) {
                                                continue;
                                                //sLow = s23;
                                                //sHigh = s13;
                                        }
                                        dt->SetPoint(count,sHigh,sLow,chPdf);
                                        count++;
                                }
                                else {
                                        dt->SetPoint(count,s13,s23,chPdf);
                                        count++;

                                }

                        }
                }
        }
        gStyle->SetPalette(1);
        dt->GetXaxis()->SetTitle("pipi");
        dt->GetYaxis()->SetTitle("pipi");
        dt->Draw("SURF1");
        dt->GetXaxis()->SetTitle("pipi");
        dt->GetYaxis()->SetTitle("pipi");
        c->SaveAs("plot_2D_"+tStrResID+"_"+label+".C");
}