No OneTemporary
Actions

Size

68 KB

Subscribers

None

View Options

	diff --git a/src/LauAbsFitModel.cc b/src/LauAbsFitModel.cc
	index 521eb7f..b07db1f 100644
	--- a/src/LauAbsFitModel.cc
	+++ b/src/LauAbsFitModel.cc
	@@ -1,1140 +1,1138 @@

	// Copyright University of Warwick 2004 - 2013.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauAbsFitModel.cc
	\brief File containing implementation of LauAbsFitModel class.
	*/

	#include <iostream>
	#include <limits>
	#include <vector>
	using std::cout;
	using std::cerr;
	using std::endl;
	using std::vector;

	#include "TMessage.h"
	#include "TMonitor.h"
	#include "TServerSocket.h"
	#include "TSocket.h"
	#include "TSystem.h"
	#include "TVectorD.h"
	#include "TVirtualFitter.h"

	#include "LauAbsFitModel.hh"
	#include "LauAbsFitter.hh"
	#include "LauAbsPdf.hh"
	#include "LauComplex.hh"
	#include "LauFitter.hh"
	#include "LauFitDataTree.hh"
	#include "LauFitNtuple.hh"
	#include "LauGenNtuple.hh"
	#include "LauParameter.hh"
	#include "LauParamFixed.hh"
	#include "LauPrint.hh"
	#include "LauSPlot.hh"

	ClassImp(LauAbsFitModel)


	LauAbsFitModel::LauAbsFitModel() :
	twoStageFit_(kFALSE),
	useAsymmFitErrors_(kFALSE),
	compareFitData_(kFALSE),
	writeLatexTable_(kFALSE),
	writeSPlotData_(kFALSE),
	storeDPEff_(kFALSE),
	randomFit_(kFALSE),
	emlFit_(kFALSE),
	poissonSmear_(kFALSE),
	enableEmbedding_(kFALSE),
	usingDP_(kTRUE),
	pdfsDependOnDP_(kFALSE),
	firstExpt_(0),
	nExpt_(0),
	evtsPerExpt_(0),
	iExpt_(0),
	inputFitData_(0),
	fitNtuple_(0),
	genNtuple_(0),
	sPlotNtuple_(0),
	fitStatus_(0),
	NLL_(0),
	numberOKFits_(0),
	numberBadFits_(0),
	nParams_(0),
	nFreeParams_(0),
	worstLogLike_(std::numeric_limits<Double_t>::max()),
	withinAsymErrorCalc_(kFALSE),
	nullString_(""),
	doSFit_(kFALSE),
	sWeightBranchName_(""),
	sWeightScaleFactor_(1.0),
	fitToyMCFileName_("fitToyMC.root"),
	fitToyMCTableName_("fitToyMCTable"),
	fitToyMCScale_(10),
	fitToyMCPoissonSmear_(kFALSE),
	sPlotFileName_(""),
	sPlotTreeName_(""),
	sPlotVerbosity_(""),
	sMaster_(0),
	messageFromMaster_(0),
	slaveId_(-1),
	nSlaves_(0),
	parValues_(0)
	{
	}


	LauAbsFitModel::~LauAbsFitModel()
	{
	delete inputFitData_; inputFitData_ = 0;
	delete fitNtuple_; fitNtuple_ = 0;
	delete genNtuple_; genNtuple_ = 0;
	delete sPlotNtuple_; sPlotNtuple_ = 0;
	delete sMaster_; sMaster_ = 0;
	delete[] parValues_; parValues_ = 0;
	}

	void LauAbsFitModel::run(const TString& applicationCode, const TString& dataFileName, const TString& dataTreeName,
	const TString& histFileName, const TString& tableFileName)
	{
	// Chose whether you want to generate or fit events in the Dalitz plot.
	// To generate events choose applicationCode = "gen", to fit events choose
	// applicationCode = "fit".

	TString runCode(applicationCode);
	runCode.ToLower();

	TString histFileNameCopy(histFileName);
	TString tableFileNameCopy(tableFileName);
	TString dataFileNameCopy(dataFileName);
	TString dataTreeNameCopy(dataTreeName);

	// Initialise the fit par vectors. Each class that inherits from this one
	// must implement this sensibly for all vectors specified in clearFitParVectors,
	// i.e. specify parameter names, initial, min, max and fixed values
	this->initialise();

	if (dataFileNameCopy == "") {dataFileNameCopy = "data.root";}
	if (dataTreeNameCopy == "") {dataTreeNameCopy = "genResults";}

	if (runCode.Contains("gen")) {

	if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
	if (tableFileNameCopy == "") {tableFileNameCopy = "genResults";}

	this->setGenValues();
	this->generate(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);

	} else if (runCode.Contains("fit")) {

	if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
	if (tableFileNameCopy == "") {tableFileNameCopy = "fitResults";}

	this->fit(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);

	} else if (runCode.Contains("reweight")) {

	this->weightEvents(dataFileNameCopy, dataTreeNameCopy);
	}
	}

	void LauAbsFitModel::runSlave(const TString& dataFileName, const TString& dataTreeName,
	const TString& histFileName, const TString& tableFileName,
	const TString& addressMaster, const UInt_t portMaster)
	{
	if ( sMaster_ != 0 ) {
	std::cerr << "ERROR in LauAbsFitModel::runSlave : master socket already present" << std::endl;
	return;
	}

	// Open connection to master
	sMaster_ = new TSocket(addressMaster, portMaster);
	sMaster_->Recv( messageFromMaster_ );
	messageFromMaster_->ReadUInt( slaveId_ );
	messageFromMaster_->ReadUInt( nSlaves_ );

	delete messageFromMaster_;
	messageFromMaster_ = 0;

	std::cout << "INFO in LauAbsFitModel::runSlave : Established connection to master on port " << portMaster << std::endl;
	std::cout << " : We are slave " << slaveId_ << " of " << nSlaves_ << std::endl;

	// Initialise the fit par vectors. Each class that inherits from this one
	// must implement this sensibly for all vectors specified in clearFitParVectors,
	// i.e. specify parameter names, initial, min, max and fixed values
	this->initialise();

	nParams_ = fitVars_.size();
	parValues_ = new Double_t[nParams_];
	for ( UInt_t iPar(0); iPar < nParams_; ++iPar ) {
	parValues_[iPar] = fitVars_[iPar]->initValue();
	}

	TString dataFileNameCopy(dataFileName);
	TString dataTreeNameCopy(dataTreeName);
	TString histFileNameCopy(histFileName);
	TString tableFileNameCopy(tableFileName);

	if (dataFileNameCopy == "") {dataFileNameCopy = "data.root";}
	if (dataTreeNameCopy == "") {dataTreeNameCopy = "genResults";}
	if (histFileNameCopy == "") {histFileNameCopy = "parInfo.root";}
	if (tableFileNameCopy == "") {tableFileNameCopy = "fitResults";}

	this->fitSlave(dataFileNameCopy, dataTreeNameCopy, histFileNameCopy, tableFileNameCopy);

	std::cout << "Fit slave " << slaveId_ << " has finished successfully" << std::endl;
	}

	void LauAbsFitModel::doSFit( const TString& sWeightBranchName, Double_t scaleFactor )
	{
	if ( sWeightBranchName == "" ) {
	cerr << "WARNING in LauAbsFitModel::doSFit : sWeight branch name is empty string, not setting-up sFit." << endl;
	return;
	}

	doSFit_ = kTRUE;
	sWeightBranchName_ = sWeightBranchName;
	sWeightScaleFactor_ = scaleFactor;
	}

	void LauAbsFitModel::setBkgndClassNames( const std::vector<TString>& names )
	{
	if ( !bkgndClassNames_.empty() ) {
	cerr << "WARNING in LauAbsFitModel::setBkgndClassNames : Names already stored, not changing them." << endl;
	return;
	}

	UInt_t nBkgnds = names.size();
	for ( UInt_t i(0); i < nBkgnds; ++i ) {
	bkgndClassNames_.insert( std::make_pair( i, names[i] ) );
	}

	this->setupBkgndVectors();
	}

	Bool_t LauAbsFitModel::validBkgndClass( const TString& className ) const
	{
	if ( bkgndClassNames_.empty() ) {
	return kFALSE;
	}

	Bool_t found(kFALSE);
	for ( LauBkgndClassMap::const_iterator iter = bkgndClassNames_.begin(); iter != bkgndClassNames_.end(); ++iter ) {
	if ( iter->second == className ) {
	found = kTRUE;
	break;
	}
	}

	return found;
	}

	UInt_t LauAbsFitModel::bkgndClassID( const TString& className ) const
	{
	if ( ! this->validBkgndClass( className ) ) {
	cerr << "ERROR in LauAbsFitModel::bkgndClassID : Request for ID for invalid background class \"" << className << "\"." << endl;
	return (bkgndClassNames_.size() + 1);
	}

	UInt_t bgID(0);
	for ( LauBkgndClassMap::const_iterator iter = bkgndClassNames_.begin(); iter != bkgndClassNames_.end(); ++iter ) {
	if ( iter->second == className ) {
	bgID = iter->first;
	break;
	}
	}

	return bgID;
	}

	const TString& LauAbsFitModel::bkgndClassName( UInt_t classID ) const
	{
	LauBkgndClassMap::const_iterator iter = bkgndClassNames_.find( classID );

	if ( iter == bkgndClassNames_.end() ) {
	cerr << "ERROR in LauAbsFitModel::bkgndClassName : Request for name of invalid background class ID " << classID << "." << endl;
	return nullString_;
	}

	return iter->second;
	}

	void LauAbsFitModel::clearFitParVectors()
	{
	cout << "Clearing fit vectors" << endl;
	fitVars_.clear();
	}

	void LauAbsFitModel::clearExtraVarVectors()
	{
	cout << "Clearing extra ntuple variable vectors" << endl;
	extraVars_.clear();
	}

	void LauAbsFitModel::setGenValues()
	{
	// makes sure each parameter holds its genValue as its current value
	for (LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter) {
	(iter)->value((iter)->genValue());
	}
	this->propagateParUpdates();
	}

	void LauAbsFitModel::writeSPlotData(const TString& fileName, const TString& treeName, Bool_t storeDPEfficiency, const TString& verbosity)
	{
	if (this->writeSPlotData()) {
	cerr << "ERROR in LauAbsFitModel::writeSPlotData : Already have an sPlot ntuple setup, not doing it again." << endl;
	return;
	}
	writeSPlotData_ = kTRUE;
	sPlotFileName_ = fileName;
	sPlotTreeName_ = treeName;
	sPlotVerbosity_ = verbosity;
	storeDPEff_ = storeDPEfficiency;
	}

	// TODO : histFileName isn't used here at the moment but could be used for
	// storing the values of the parameters used in the generation.
	// These could then be read and used for setting the "true" values
	// in a subsequent fit.
	void LauAbsFitModel::generate(const TString& dataFileName, const TString& dataTreeName, const TString& /histFileName/, const TString& tableFileNameBase)
	{
	// Create the ntuple for storing the results
	cout << "Creating generation ntuple." << endl;
	if (genNtuple_ != 0) {delete genNtuple_; genNtuple_ = 0;}
	genNtuple_ = new LauGenNtuple(dataFileName,dataTreeName);

	// add branches for storing the experiment number and the number of
	// the event within the current experiment
	this->addGenNtupleIntegerBranch("iExpt");
	this->addGenNtupleIntegerBranch("iEvtWithinExpt");
	this->setupGenNtupleBranches();

	// Start the cumulative timer
	cumulTimer_.Start();

	Bool_t genOK(kTRUE);
	do {
	// Loop over the number of experiments
	for (iExpt_ = firstExpt_; iExpt_ < (firstExpt_+nExpt_); ++iExpt_) {

	// Start the timer to see how long each experiment takes to generate
	timer_.Start();

	// Store the experiment number in the ntuple
	this->setGenNtupleIntegerBranchValue("iExpt",iExpt_);

	// Do the generation for this experiment
	cout << "Generating experiment number " << iExpt_ << endl;
	genOK = this->genExpt();

	// Stop the timer and see how long the program took so far
	timer_.Stop();
	timer_.Print();

	if (!genOK) {
	// delete and recreate an empty tree
	genNtuple_->deleteAndRecreateTree();

	// then break out of the experiment loop
	cerr << "ERROR in LauAbsFitModel::generate : Problem in toy MC generation. Starting again with updated parameters..." << endl;
	break;
	}

	if (this->writeLatexTable()) {
	TString tableFileName(tableFileNameBase);
	tableFileName += "_";
	tableFileName += iExpt_;
	tableFileName += ".tex";
	this->writeOutTable(tableFileName);
	}

	} // Loop over number of experiments
	} while (!genOK);

	// Print out total timing info.
	cumulTimer_.Stop();
	cout << "Finished generating all experiments." << endl;
	cout << "Cumulative timing:" << endl;
	cumulTimer_.Print();

	// Build the event index
	cout << "Building experiment:event index." << endl;
	// TODO - can test this return value?
	//Int_t nIndexEntries =
	genNtuple_->buildIndex("iExpt","iEvtWithinExpt");

	// Write out toy MC ntuple
	cout << "Writing data to file " << dataFileName << "." << endl;
	genNtuple_->writeOutGenResults();
	}

	void LauAbsFitModel::addGenNtupleIntegerBranch(const TString& name)
	{
	genNtuple_->addIntegerBranch(name);
	}

	void LauAbsFitModel::addGenNtupleDoubleBranch(const TString& name)
	{
	genNtuple_->addDoubleBranch(name);
	}

	void LauAbsFitModel::setGenNtupleIntegerBranchValue(const TString& name, Int_t value)
	{
	genNtuple_->setIntegerBranchValue(name,value);
	}

	void LauAbsFitModel::setGenNtupleDoubleBranchValue(const TString& name, Double_t value)
	{
	genNtuple_->setDoubleBranchValue(name,value);
	}

	Int_t LauAbsFitModel::getGenNtupleIntegerBranchValue(const TString& name) const
	{
	return genNtuple_->getIntegerBranchValue(name);
	}

	Double_t LauAbsFitModel::getGenNtupleDoubleBranchValue(const TString& name) const
	{
	return genNtuple_->getDoubleBranchValue(name);
	}

	void LauAbsFitModel::fillGenNtupleBranches()
	{
	genNtuple_->fillBranches();
	}

	void LauAbsFitModel::addSPlotNtupleIntegerBranch(const TString& name)
	{
	sPlotNtuple_->addIntegerBranch(name);
	}

	void LauAbsFitModel::addSPlotNtupleDoubleBranch(const TString& name)
	{
	sPlotNtuple_->addDoubleBranch(name);
	}

	void LauAbsFitModel::setSPlotNtupleIntegerBranchValue(const TString& name, Int_t value)
	{
	sPlotNtuple_->setIntegerBranchValue(name,value);
	}

	void LauAbsFitModel::setSPlotNtupleDoubleBranchValue(const TString& name, Double_t value)
	{
	sPlotNtuple_->setDoubleBranchValue(name,value);
	}

	void LauAbsFitModel::fillSPlotNtupleBranches()
	{
	sPlotNtuple_->fillBranches();
	}

	void LauAbsFitModel::fit(const TString& dataFileName, const TString& dataTreeName, const TString& histFileName, const TString& tableFileNameBase)
	{
	// Routine to perform the total fit.

	cout << "INFO in LauAbsFitModel::fit : First experiment = " << firstExpt_ << endl;
	cout << "INFO in LauAbsFitModel::fit : Number of experiments = " << nExpt_ << endl;

	// Start the cumulative timer
	cumulTimer_.Start();

	numberOKFits_ = 0, numberBadFits_ = 0;
	fitStatus_ = -1;

	// Create and setup the fit results ntuple
	cout << "INFO in LauAbsFitModel::fit : Creating fit ntuple." << endl;
	if (fitNtuple_ != 0) {delete fitNtuple_; fitNtuple_ = 0;}
	fitNtuple_ = new LauFitNtuple(histFileName, this->useAsymmFitErrors());

	// Create and setup the sPlot ntuple
	if (this->writeSPlotData()) {
	cout << "INFO in LauAbsFitModel::fit : Creating sPlot ntuple." << endl;
	if (sPlotNtuple_ != 0) {delete sPlotNtuple_; sPlotNtuple_ = 0;}
	sPlotNtuple_ = new LauGenNtuple(sPlotFileName_,sPlotTreeName_);
	this->setupSPlotNtupleBranches();
	}

	// 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->cacheFitData(dataFileName,dataTreeName);
	if (!dataOK) {
	cerr << "ERROR in LauAbsFitModel::fit : Problem caching the fit data." << endl;
	gSystem->Exit(EXIT_FAILURE);
	}

	// Loop over the number of experiments
	for (iExpt_ = firstExpt_; iExpt_ < (firstExpt_+nExpt_); ++iExpt_) {

	// Start the timer to see how long each fit takes
	timer_.Start();

	inputFitData_->readExperimentData(iExpt_);
	this->eventsPerExpt(inputFitData_->nEvents());

	if (this->eventsPerExpt() < 1) {
	cerr << "ERROR in LauAbsFitModel::fit : Zero events in experiment " << iExpt_ << ", skipping..." << endl;
	timer_.Stop();
	continue;
	}

	// Now the sub-classes must implement whatever they need to do
	// to cache any more input fit data they need in order to
	// calculate the likelihoods during the fit.
	// They need to use the inputFitData_ tree as input. For example,
	// inputFitData_ contains m13Sq and m23Sq. The appropriate fit model
	// then caches the resonance dynamics for the signal model, as
	// well as the background likelihood values in the Dalitz plot
	this->cacheInputFitVars();
	if ( this->doSFit() ) {
	this->cacheInputSWeights();
	}

	// Do the fit for this experiment
	this->fitExpt();

	// Write the results into the ntuple
	this->finaliseFitResults(tableFileNameBase);

	// Stop the timer and see how long the program took so far
	timer_.Stop();
	timer_.Print();

	// Store the per-event likelihood values
	if ( this->writeSPlotData() ) {
	this->storePerEvtLlhds();
	}

	// Create a toy MC sample for the 1st successful experiment
	// using the fitted parameters so that the user can compare
	// the fit to the actual data. The toy MC stats are a
	// factor of 10 higher than the number of events specified
	// via setNEvGen. This is to reduce the statistical
	// fluctuations for the toy MC data.
	if (compareFitData_ == kTRUE && fitStatus_ == 3) {
	this->createFitToyMC(fitToyMCFileName_, fitToyMCTableName_);
	compareFitData_ = kFALSE; // only do this for the first successful experiment
	}

	// Keep track of how many fits worked or failed
	// NB values of fitStatus_ now indicate the status of the error matrix:
	// 0= not calculated at all
	// 1= approximation only, not accurate
	// 2= full matrix, but forced positive-definite
	// 3= full accurate covariance matrix
	if (fitStatus_ == 3) {
	numberOKFits_++;
	} else {
	numberBadFits_++;
	}

	} // Loop over number of experiments

	// Print out total timing info.
	cumulTimer_.Stop();
	cout << "INFO in LauAbsFitModel::fit : Cumulative timing:" << endl;
	cumulTimer_.Print();

	// Print out stats on OK fits.
	cout << "INFO in LauAbsFitModel::fit : Number of OK Fits = " << numberOKFits_ << endl;
	cout << "INFO in LauAbsFitModel::fit : Number of Failed Fits = " << numberBadFits_ << endl;
	Double_t fitEff(0.0);
	if (nExpt_ != 0) {fitEff = numberOKFits_/(1.0*nExpt_);}
	cout << "INFO in LauAbsFitModel::fit : Fit efficiency = " << fitEff*100.0 << "%." << endl;

	// Write out any fit results (ntuples etc...).
	this->writeOutAllFitResults();
	if ( this->writeSPlotData() ) {
	this->calculateSPlotData();
	}
	}

	void LauAbsFitModel::fitSlave(const TString& dataFileName, const TString& dataTreeName, const TString& histFileName, const TString& tableFileNameBase)
	{
	// Create and setup the fit results ntuple
	cout << "INFO in LauAbsFitModel::fitSlave : Creating fit ntuple." << endl;
	if (fitNtuple_ != 0) {delete fitNtuple_; fitNtuple_ = 0;}
	fitNtuple_ = new LauFitNtuple(histFileName, this->useAsymmFitErrors());

	// 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->cacheFitData(dataFileName,dataTreeName);
	if (!dataOK) {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Problem caching the fit data." << std::endl;
	gSystem->Exit(EXIT_FAILURE);
	}

	// Now process the various requests from the master

	TMessage messageToMaster(kMESS_ANY);

	while ( kTRUE ) {

	sMaster_->Recv( messageFromMaster_ );

	if ( messageFromMaster_->What() == kMESS_STRING ) {

	TString msgStr;
	messageFromMaster_->ReadTString( msgStr );

	std::cout << "INFO in LauAbsFitModel::fitSlave : Received message from master: " << msgStr << std::endl;

	if ( msgStr == "Send Parameters" ) {

	// Update initial fit parameters if required (e.g. if using random numbers).
	this->checkInitFitParams();

	// Send the fit parameters
	TObjArray array;
	for ( LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter ) {
	array.Add( *iter );
	}

	messageToMaster.Reset( kMESS_OBJECT );
	messageToMaster.WriteObject( &array );
	sMaster_->Send( messageToMaster );

	} else if ( msgStr == "Read Expt" ) {

	// Read the data for this experiment
	messageFromMaster_->ReadUInt( iExpt_ );

	inputFitData_->readExperimentData( iExpt_ );
	UInt_t nEvent = inputFitData_->nEvents();
	this->eventsPerExpt( nEvent );

	if ( nEvent < 1 ) {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Zero events in experiment " << firstExpt_ << ", the master should skip this experiment..." << std::endl;
	}

	messageToMaster.Reset( kMESS_ANY );
	messageToMaster.WriteUInt( slaveId_ );
	messageToMaster.WriteUInt( nEvent );
	sMaster_->Send( messageToMaster );

	} else if ( msgStr == "Cache" ) {

	// Perform the caching

	this->cacheInputFitVars();

	messageToMaster.Reset( kMESS_ANY );
	messageToMaster.WriteUInt( slaveId_ );
	messageToMaster.WriteBool( kTRUE );
	sMaster_->Send( messageToMaster );

	} else if ( msgStr == "Write Results" ) {

	this->writeOutAllFitResults();

	messageToMaster.Reset( kMESS_ANY );
	messageToMaster.WriteUInt( slaveId_ );
	messageToMaster.WriteBool( kTRUE );
	sMaster_->Send( messageToMaster );

	} else if ( msgStr == "Finish" ) {

	std::cout << "INFO in LauAbsFitModel::fitSlave : Message from master to finish" << std::endl;
	break;
	} else {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Unexpected message from master" << std::endl;
	gSystem->Exit( EXIT_FAILURE );
	}

	} else if ( messageFromMaster_->What() == kMESS_OBJECT ) {

	std::cout << "INFO in LauAbsFitModel::fitSlave : Received message from master: Finalise" << std::endl;

	messageFromMaster_->ReadInt( fitStatus_ );
	messageFromMaster_->ReadDouble( NLL_ );

	TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromMaster_->ReadObject( messageFromMaster_->GetClass() ) );
	if ( ! objarray ) {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Error reading parameters from master" << std::endl;
	gSystem->Exit( EXIT_FAILURE );
	}

	TMatrixD * covMat = dynamic_cast<TMatrixD*>( messageFromMaster_->ReadObject( messageFromMaster_->GetClass() ) );
	if ( ! covMat ) {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Error reading covariance matrix from master" << std::endl;
	gSystem->Exit( EXIT_FAILURE );
	}
	covMatrix_.Clear();
	covMatrix_.ResizeTo( covMat->GetNrows(), covMat->GetNcols() );
	covMatrix_.SetMatrixArray( covMat->GetMatrixArray() );
	delete covMat; covMat = 0;


	UInt_t nPars = objarray->GetEntries();

	if ( nPars != nParams_ ) {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Unexpected number of parameters received from master" << std::endl;
	std::cerr << " ::fitSlave : Received " << nPars << " when expecting " << nParams_ << std::endl;
	gSystem->Exit( EXIT_FAILURE );
	}

	for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
	LauParameter* parameter = dynamic_cast<LauParameter>( (objarray)[iPar] );
	if ( ! parameter ) {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Error reading parameter from master" << std::endl;
	gSystem->Exit( EXIT_FAILURE );
	}

	if ( parameter->name() != fitVars_[iPar]->name() ) {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Error reading parameter from master" << std::endl;
	gSystem->Exit( EXIT_FAILURE );
	}

	(fitVars_[iPar]) = parameter;
	}

	this->finaliseFitResults( tableFileNameBase );

	// Send the finalised parameters
	TObjArray array;
	for ( LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter ) {
	array.Add( *iter );
	}

	messageToMaster.Reset( kMESS_ANY );
	messageToMaster.WriteUInt( slaveId_ );
	messageToMaster.WriteBool( kTRUE );
	messageToMaster.WriteObject( &array );
	sMaster_->Send( messageToMaster );

	} else if ( messageFromMaster_->What() == kMESS_ANY ) {

	UInt_t nPars(0);
	messageFromMaster_->ReadUInt( nPars );

	if ( nPars != nParams_ ) {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Unexpected number of parameters received from master" << std::endl;
	std::cerr << " ::fitSlave : Received " << nPars << " when expecting " << nParams_ << std::endl;
	gSystem->Exit( EXIT_FAILURE );
	}

	messageFromMaster_->ReadFastArray( parValues_, nPars );

	for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
	if ( ! fitVars_[iPar]->fixed() ) {
	fitVars_[iPar]->value( parValues_[iPar] );
	}
	}
	this->propagateParUpdates();

	Double_t negLogLike = this->getTotNegLogLikelihood();

	messageToMaster.Reset( kMESS_ANY );
	messageToMaster.WriteDouble( negLogLike );
	sMaster_->Send( messageToMaster );

	} else {
	std::cerr << "ERROR in LauAbsFitModel::fitSlave : Unexpected message type" << std::endl;
	gSystem->Exit( EXIT_FAILURE );
	}

	delete messageFromMaster_;
	messageFromMaster_ = 0;
	}
	}

	Bool_t LauAbsFitModel::cacheFitData(const TString& dataFileName, const TString& dataTreeName)
	{
	// From the input data stream, store the variables into the
	// internal tree inputFitData_ that can be used by the sub-classes
	// in calculating their likelihood functions for the fit
	delete inputFitData_;
	inputFitData_ = new LauFitDataTree(dataFileName,dataTreeName);
	Bool_t dataOK = inputFitData_->findBranches();

	if (!dataOK) {
	delete inputFitData_; inputFitData_ = 0;
	}

	return dataOK;
	}

	void LauAbsFitModel::cacheInputSWeights()
	{
	Bool_t hasBranch = inputFitData_->haveBranch( sWeightBranchName_ );
	if ( ! hasBranch ) {
	std::cerr << "ERROR in LauAbsFitModel::cacheInputSWeights : Input data does not contain variable \"" << sWeightBranchName_ << "\".\n";
	std::cerr << " : Turning off sFit!" << std::endl;
	doSFit_ = kFALSE;
	return;
	}

	UInt_t nEvents = this->eventsPerExpt();
	sWeights_.clear();
	sWeights_.reserve( nEvents );

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

	const LauFitData& dataValues = inputFitData_->getData(iEvt);

	LauFitData::const_iterator iter = dataValues.find( sWeightBranchName_ );

	sWeights_.push_back( iter->second * sWeightScaleFactor_ );
	}
	}

	void LauAbsFitModel::fitExpt()
	{
	// Routine to perform the actual fit for the given experiment

	// Reset the worst likelihood found to its catch-all value
	worstLogLike_ = std::numeric_limits<Double_t>::max();

	// Update initial fit parameters if required (e.g. if using random numbers).
	this->checkInitFitParams();

	// Initialise the fitter
	LauFitter::fitter()->useAsymmFitErrors( this->useAsymmFitErrors() );
	LauFitter::fitter()->twoStageFit( this->twoStageFit() );
	LauFitter::fitter()->initialise( this, fitVars_ );

	nParams_ = LauFitter::fitter()->nParameters();
	nFreeParams_ = LauFitter::fitter()->nFreeParameters();

	// Now ready for minimisation step
	cout << "\nINFO in LauAbsFitModel::fitExpt : Start minimisation...\n";
	std::pair<Int_t,Double_t> fitResult = LauFitter::fitter()->minimise();

	fitStatus_ = fitResult.first;
	NLL_ = fitResult.second;

	// If we're doing a two stage fit we can now release (i.e. float)
	// the 2nd stage parameters and re-fit
	if (this->twoStageFit()) {

	if ( fitStatus_ != 3 ) {
	cerr << "ERROR in LauAbsFitModel:fitExpt : Not running second stage fit since first stage failed." << endl;
	LauFitter::fitter()->releaseSecondStageParameters();
	} else {
	LauFitter::fitter()->fixFirstStageParameters();
	LauFitter::fitter()->releaseSecondStageParameters();
	nParams_ = LauFitter::fitter()->nParameters();
	nFreeParams_ = LauFitter::fitter()->nFreeParameters();
	fitResult = LauFitter::fitter()->minimise();
	}
	}

	fitStatus_ = fitResult.first;
	NLL_ = fitResult.second;
	const TMatrixD& covMat = LauFitter::fitter()->covarianceMatrix();
	covMatrix_.Clear();
	covMatrix_.ResizeTo( covMat.GetNrows(), covMat.GetNcols() );
	covMatrix_.SetMatrixArray( covMat.GetMatrixArray() );

	// Store the final fit results and errors into protected internal vectors that
	// all sub-classes can use within their own finalFitResults implementation
	// used below (e.g. putting them into an ntuple in a root file)
	LauFitter::fitter()->updateParameters();
	LauFitter::fitter()->releaseFirstStageParameters();
	}

	void LauAbsFitModel::writeOutAllFitResults()
	{
	// Write out histograms at end
	if (fitNtuple_ != 0) {
	fitNtuple_->writeOutFitResults();
	}
	}

	void LauAbsFitModel::calculateSPlotData()
	{
	if (sPlotNtuple_ != 0) {
	sPlotNtuple_->addFriendTree(inputFitData_->fileName(), inputFitData_->treeName());
	sPlotNtuple_->writeOutGenResults();
	LauSPlot splot(sPlotNtuple_->fileName(), sPlotNtuple_->treeName(), this->firstExpt(), this->nExpt(),
	this->variableNames(), this->freeSpeciesNames(), this->fixdSpeciesNames(), this->twodimPDFs(),
	this->splitSignal(), this->scfDPSmear());
	splot.runCalculations(sPlotVerbosity_);
	splot.writeOutResults();
	}
	}

	void LauAbsFitModel::compareFitData(UInt_t toyMCScale, const TString& mcFileName, const TString& tableFileName, Bool_t poissonSmearing)
	{
	compareFitData_ = kTRUE;
	fitToyMCScale_ = toyMCScale;
	fitToyMCFileName_ = mcFileName;
	fitToyMCTableName_ = tableFileName;
	fitToyMCPoissonSmear_ = poissonSmearing;
	}

	void LauAbsFitModel::createFitToyMC(const TString& mcFileName, const TString& tableFileName)
	{
	// Create a toy MC sample so that the user can eventually
	// compare the "fitted" result with the data
	// Generate more toy MC to reduce statistical fluctuations. Use the
	// rescaling value fitToyMCScale_.

	// Store the info on the number of experiments, first expt and current expt
	UInt_t oldNExpt(this->nExpt());
	UInt_t oldFirstExpt(this->firstExpt());
	UInt_t oldIExpt(iExpt_);

	// Turn off Poisson smearing if required
	Bool_t poissonSmearing(this->doPoissonSmearing());
	this->doPoissonSmearing(fitToyMCPoissonSmear_);

	// Turn off embedding, since we need toy MC, not reco'ed events
	Bool_t enableEmbeddingOrig(this->enableEmbedding());
	this->enableEmbedding(kFALSE);

	// Need to make sure that the generation of the DP co-ordinates is
	// switched on if any of our PDFs depend on it
	Bool_t origUseDP = this->useDP();
	if ( this->pdfsDependOnDP() && !origUseDP ) {
	this->useDP( kTRUE );
	this->initialiseDPModels();
	}

	// Generate the toy MC
	cout << "Generating toy MC in " << mcFileName << " to compare fit with data..." << endl;
	cout << "Number of experiments to generate = " << this->nExpt() << ", which is a factor of "
	<<fitToyMCScale_ << " higher than that originally specified. "
	<<"This is to allow the toy MC to be made with reduced statistical "
	<<"fluctuations." << endl;

	// Set the genValue of each parameter to its current (fitted) value
	// but first store the original genValues for restoring later
	vector<Double_t> origGenValues; origGenValues.reserve(nParams_);
	for (LauParameterPList::iterator iter = fitVars_.begin(); iter != fitVars_.end(); ++iter) {
	origGenValues.push_back((*iter)->genValue());
	(iter)->genValue((iter)->value());
	}

	// If we're asked to generate more than 100 experiments then split it
	// up into multiple files since otherwise can run into memory issues
	// when building the index

	UInt_t totalExpts = oldNExpt * fitToyMCScale_;
	if ( totalExpts > 100 ) {
	UInt_t nFiles = totalExpts/100;
	if ( totalExpts%100 ) {
	nFiles += 1;
	}
	for ( UInt_t iFile(0); iFile < nFiles; ++iFile ) {

	UInt_t firstExp( iFile*100 );

	// Set number of experiments and first experiment to generate
	UInt_t nExp = ((firstExp + 100)>totalExpts) ? totalExpts-firstExp : 100;
	this->setNExpts(nExp, firstExp);

	// Create a unique filename and generate the events
	TString extraname = "_file";
	extraname += iFile;
	TString filename(mcFileName);
	filename.Insert( filename.Index("."), extraname );
	this->generate(filename, "genResults", "dummy.root", tableFileName);
	}
	} else {
	// Set number of experiments to new value
	this->setNExpts(oldNExpt*fitToyMCScale_, 0);
	// Generate the toy
	this->generate(mcFileName, "genResults", "dummy.root", tableFileName);
	}

	// Reset number of experiments to original value
	iExpt_ = oldIExpt;
	this->setNExpts(oldNExpt, oldFirstExpt);

	// Restore the Poisson smearing to its former value
	this->doPoissonSmearing(poissonSmearing);

	// Restore the embedding status
	this->enableEmbedding(enableEmbeddingOrig);

	// Restore "useDP" to its former status
	this->useDP( origUseDP );

	// Restore the original genValue to each parameter
	for (UInt_t i(0); i<nParams_; ++i) {
	fitVars_[i]->genValue(origGenValues[i]);
	}

	cout << "Finished in createFitToyMC." << endl;
	}

	Double_t LauAbsFitModel::getTotNegLogLikelihood()
	{
	// Calculate the total negative log-likelihood over all events.
	// This function assumes that the fit parameters and data tree have
	// already been set-up correctly.

	// Loop over the data points to calculate the log likelihood
	Double_t logLike = this->getLogLikelihood( 0, this->eventsPerExpt() );

	// Include the Poisson term in the extended likelihood if required
	if (this->doEMLFit()) {
	logLike -= this->getEventSum();
	}

	Double_t totNegLogLike = -logLike;
	return totNegLogLike;
	}

	Double_t LauAbsFitModel::getLogLikelihood( UInt_t iStart, UInt_t iEnd )
	{
	// Calculate the total negative log-likelihood over all events.
	// This function assumes that the fit parameters and data tree have
	// already been set-up correctly.

	// Loop over the data points to calculate the log likelihood
	Double_t logLike(0.0);

	// Loop over the number of events in this experiment
	Bool_t ok(kTRUE);
	for (UInt_t iEvt = iStart; iEvt < iEnd; ++iEvt) {

	Double_t likelihood = this->getTotEvtLikelihood(iEvt);

	if (likelihood > DBL_MIN) { // Is the likelihood zero?
	Double_t evtLogLike = TMath::Log(likelihood);
	if ( doSFit_ ) {
	evtLogLike *= sWeights_[iEvt];
	}
	logLike += evtLogLike;
	} else {
	ok = kFALSE;
	cerr << "WARNING in LauAbsFitModel::getLogLikelihood : Strange likelihood value for event " << iEvt << ": " << likelihood << "\n";
	this->printEventInfo(iEvt);
	this->printVarsInfo(); //Write the values of the floated variables for which the likelihood is zero
	break;
	}
	}

	if (!ok) {
	cerr << " : Returning worst NLL found so far to force MINUIT out of this region." << endl;
	logLike = worstLogLike_;
	} else if (logLike < worstLogLike_) {
	worstLogLike_ = logLike;
	}

	return logLike;
	}

	void LauAbsFitModel::setParsFromMinuit(Double_t* par, Int_t npar)
	{
	// This function sets the internal parameters based on the values
	// that Minuit is using when trying to minimise the total likelihood function.

	- // MINOS reports one fewer free parameter than there actually is,
	- // so increment npar so the following check doesn't fail
	- if ( withinAsymErrorCalc_ ) {
	- ++npar;
	- }
	-
	- if (static_cast<UInt_t>(npar) != nFreeParams_) {
	- cerr << "ERROR in LauAbsFitModel::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n";
	- cerr << " Expected: " << nFreeParams_ << ".\n" << endl;
	- gSystem->Exit(EXIT_FAILURE);
	+ // MINOS reports different numbers of free parameters depending on the
	+ // situation, so disable this check
	+ if ( ! withinAsymErrorCalc_ ) {
	+ if (static_cast<UInt_t>(npar) != nFreeParams_) {
	+ cerr << "ERROR in LauAbsFitModel::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n";
	+ cerr << " Expected: " << nFreeParams_ << ".\n" << endl;
	+ gSystem->Exit(EXIT_FAILURE);
	+ }
	}

	// Despite npar being the number of free parameters
	// the par array actually contains all the parameters,
	// free and floating...
	// Update all the floating ones with their new values.
	for (UInt_t i(0); i<nParams_; ++i) {
	if (!fitVars_[i]->fixed()) {
	fitVars_[i]->value(par[i]);
	}
	}

	this->propagateParUpdates();
	}

	UInt_t LauAbsFitModel::addFitParameters(LauPdfList& pdfList)
	{
	UInt_t nParsAdded(0);
	for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
	LauAbsPdf* pdf = (*pdf_iter);
	if ( pdf->isDPDependent() ) {
	this->pdfsDependOnDP( kTRUE );
	}
	LauParameterPList& pars = pdf->getParameters();
	for (LauParameterPList::iterator pars_iter = pars.begin(); pars_iter != pars.end(); ++pars_iter) {
	if ( !(pars_iter)->clone() && ( !(pars_iter)->fixed() \|\|
	(this->twoStageFit() && ( (pars_iter)->secondStage() \|\| (pars_iter)->firstStage())) ) ) {
	fitVars_.push_back(*pars_iter);
	++nParsAdded;
	}
	}
	}
	return nParsAdded;
	}

	void LauAbsFitModel::updateFitParameters(LauPdfList& pdfList)
	{
	for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
	(*pdf_iter)->updatePulls();
	}
	}

	void LauAbsFitModel::printFitParameters(const LauPdfList& pdfList, std::ostream& fout) const
	{
	LauPrint print;
	for (LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
	const LauParameterPList& pars = (*pdf_iter)->getParameters();
	for (LauParameterPList::const_iterator pars_iter = pars.begin(); pars_iter != pars.end(); ++pars_iter) {
	if (!(*pars_iter)->clone()) {
	fout << (*pars_iter)->name() << " & $";
	print.printFormat(fout, (*pars_iter)->value());
	if ((*pars_iter)->fixed() == kTRUE) {
	fout << "$ (fixed) \\\\";
	} else {
	fout << " \\pm ";
	print.printFormat(fout, (*pars_iter)->error());
	fout << "$ \\\\" << endl;
	}
	}
	}
	}
	}

	void LauAbsFitModel::cacheInfo(LauPdfList& pdfList, const LauFitDataTree& theData)
	{
	for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
	(*pdf_iter)->cacheInfo(theData);
	}
	}

	Double_t LauAbsFitModel::prodPdfValue(LauPdfList& pdfList, UInt_t iEvt)
	{
	Double_t pdfVal = 1.0;
	for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) {
	(*pdf_iter)->calcLikelihoodInfo(iEvt);
	pdfVal = (pdf_iter)->getLikelihood();
	}
	return pdfVal;
	}

	void LauAbsFitModel::printEventInfo(UInt_t iEvt) const
	{
	const LauFitData& data = inputFitData_->getData(iEvt);
	cerr << " : Input data values for this event:" << endl;
	for (LauFitData::const_iterator iter = data.begin(); iter != data.end(); ++iter) {
	cerr << " " << iter->first << " = " << iter->second << endl;
	}
	}

	void LauAbsFitModel::printVarsInfo() const
	{
	cerr << " : Current values of fit parameters:" << endl;
	for (UInt_t i(0); i<nParams_; ++i) {
	cerr << " " << (fitVars_[i]->name()).Data() << " = " << fitVars_[i]->value() << endl;
	}
	}

	diff --git a/src/LauSimFitMaster.cc b/src/LauSimFitMaster.cc
	index 99d3e17..9dac301 100644
	--- a/src/LauSimFitMaster.cc
	+++ b/src/LauSimFitMaster.cc
	@@ -1,897 +1,894 @@

	// Copyright University of Warwick 2013 - 2013.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauSimFitMaster.cc
	\brief File containing implementation of LauSimFitMaster class.
	*/

	#include <cstdlib>
	#include <iostream>

	#include "TMessage.h"
	#include "TMonitor.h"
	#include "TObjArray.h"
	#include "TObjString.h"
	#include "TServerSocket.h"
	#include "TSocket.h"
	#include "TSystem.h"
	#include "TVectorD.h"

	#include "LauAbsFitter.hh"
	#include "LauFitNtuple.hh"
	#include "LauFitter.hh"
	#include "LauParameter.hh"
	#include "LauParamFixed.hh"
	#include "LauSimFitMaster.hh"


	ClassImp(LauSimFitMaster)


	LauSimFitMaster::LauSimFitMaster( UInt_t numSlaves, UInt_t port ) :
	nSlaves_(numSlaves),
	reqPort_(port),
	nParams_(0),
	nFreeParams_(0),
	withinAsymErrorCalc_(kFALSE),
	numberOKFits_(0),
	numberBadFits_(0),
	fitStatus_(0),
	iExpt_(0),
	NLL_(0.0),
	socketMonitor_(0),
	messageFromSlave_(0),
	fitNtuple_(0)
	{
	messagesToSlaves_.resize( nSlaves_ );
	for ( UInt_t iSlave(0); iSlave < nSlaves_; ++iSlave ) {
	messagesToSlaves_[iSlave] = new TMessage();
	}
	}

	LauSimFitMaster::~LauSimFitMaster()
	{
	delete socketMonitor_; socketMonitor_ = 0;
	TString msgStr("Finish");
	TMessage message( kMESS_STRING );
	message.WriteTString(msgStr);
	for ( std::vector<TSocket*>::iterator iter = sSlaves_.begin(); iter != sSlaves_.end(); ++iter ) {
	(*iter)->Send(message);
	(*iter)->Close();
	delete (*iter);
	}
	sSlaves_.clear();
	for ( std::vector<LauParameter*>::iterator iter = params_.begin(); iter != params_.end(); ++iter ) {
	delete *iter;
	}
	params_.clear();
	for ( std::vector<Double_t*>::iterator iter = vectorPar_.begin(); iter != vectorPar_.end(); ++iter ) {
	delete[] (*iter);
	}
	vectorPar_.clear();
	delete messageFromSlave_; messageFromSlave_ = 0;
	for ( std::vector<TMessage*>::iterator iter = messagesToSlaves_.begin(); iter != messagesToSlaves_.end(); ++iter ) {
	delete (*iter);
	}
	messagesToSlaves_.clear();
	delete fitNtuple_;
	}

	void LauSimFitMaster::initSockets()
	{
	if ( socketMonitor_ != 0 ) {
	std::cerr << "ERROR in LauSimFitMaster::initSockets : Sockets already initialised." << std::endl;
	return;
	}

	//initialise socket connection, then accept a connection and return a full-duplex communication socket.
	socketMonitor_ = new TMonitor();

	TServerSocket *ss = new TServerSocket( reqPort_, kFALSE );
	UInt_t actual_port = ss->GetLocalPort();

	std::cout << "INFO in LauSimFitMaster::initSockets : Waiting for connection with " << nSlaves_ << " workers on port " << actual_port << std::endl;

	sSlaves_.resize(nSlaves_);
	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
	sSlaves_[iSlave] = ss->Accept();
	std::cout << " : Added slave " << iSlave << std::endl;
	}

	// tell the clients to start
	std::cout << "INFO in LauSimFitMaster::initSockets : Initialising slaves" << std::endl;
	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {

	TMessage message( kMESS_ANY );
	message.WriteUInt(iSlave);
	message.WriteUInt(nSlaves_);

	sSlaves_[iSlave]->Send(message);

	socketMonitor_->Add(sSlaves_[iSlave]);
	}
	std::cout << " : Now start fit\n" << std::endl;

	ss->Close();
	delete ss;
	}

	/*
	* OLD VERSION THAT JUST GETS THE NAMES - COULD HAVE A SERIES OF EXCHANGES TO GET THE NAMES, INIT VALUES, RANGES, ETC. INSTEAD OF PASSING PARAMETERS
	* THIS INCREASES THE GENERALITY OF THE CODE, I.E. THERE IS NO NEED FOR THE SLAVES TO KNOW ANY LAURA++ CLASS BUT THIS ONE, BUT MAKES IT RATHER MORE DENSE
	* FOR THE MOMENT I WILL STICK WITH THE METHOD OF PASSING LAUPARAMETER OBJECTS AROUND AND CONSIDER GOING BACK TO THIS GENERAL METHOD ONCE EVERYTHING IS WORKING
	*
	void LauSimFitMaster::getParametersFromSlavesFirstTime()
	{
	slaveIndices_.resize( nSlaves_ );

	TSocket* sActive(0);
	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
	// Send a message to the slave, requesting the list of parameter names
	TString msgStr = "Parameter Names";
	TMessage message( kMESS_STRING );
	message.WriteTString( msgStr );
	sSlaves_[iSlave]->Send(message);

	// Wait to receive the response and check that it has come from the slave we just requested from
	sActive = socketMonitor_->Select();
	if ( sActive != sSlaves_[iSlave] ) {
	std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Received message from a different slave than expected!" << std::endl;
	gSystem->Exit(1);
	}

	// Read the object and extract the parameter names
	sSlaves_[iSlave]->Recv( messageFromSlave_ );
	TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
	if ( ! objarray ) {
	std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameter names from slave" << std::endl;
	gSystem->Exit(1);
	}

	Int_t nPars = objarray->GetEntries();
	for ( Int_t iPar(0); iPar < nPars; ++iPar ) {
	TObjString* objstring = dynamic_cast<TObjString>( (objarray)[iPar] );
	if ( ! objstring ) {
	std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameter names from slave" << std::endl;
	gSystem->Exit(1);
	}
	TString parname = objstring->GetString();

	std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
	if ( iter != parIndices_.end() ) {
	UInt_t index = iter->second;
	slaveIndices_[iSlave].push_back( index );
	} else {
	UInt_t index = parIndices_.size();
	parIndices_.insert( std::make_pair( parname, index ) );
	parNames_.insert( std::make_pair( index, parname ) );
	slaveIndices_[iSlave].push_back( index );
	}
	}

	delete objarray; objarray = 0;
	delete messageFromSlave_; messageFromSlave_ = 0;
	}

	UInt_t nPars = parNames_.size();
	parValues_.resize( nPars );
	}
	*/

	void LauSimFitMaster::getParametersFromSlaves()
	{
	if ( socketMonitor_ == 0 ) {
	std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlaves : Sockets not initialised." << std::endl;
	return;
	}

	if ( params_.empty() ) {
	this->getParametersFromSlavesFirstTime();
	} else {
	this->updateParametersFromSlaves();
	}
	}

	void LauSimFitMaster::updateParametersFromSlaves()
	{
	TSocket* sActive(0);

	// Construct a message, requesting the list of parameter names
	TString msgStr = "Send Parameters";
	TMessage message( kMESS_STRING );
	message.WriteTString( msgStr );

	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
	// Send the message to the slave
	sSlaves_[iSlave]->Send(message);

	// Wait to receive the response and check that it has come from the slave we just requested from
	sActive = socketMonitor_->Select();
	if ( sActive != sSlaves_[iSlave] ) {
	std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Received message from a different slave than expected!" << std::endl;
	gSystem->Exit(1);
	}

	// Read the object and extract the parameter names
	sSlaves_[iSlave]->Recv( messageFromSlave_ );
	TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
	if ( ! objarray ) {
	std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Error reading parameter names from slave" << std::endl;
	gSystem->Exit(1);
	}

	// We want to auto-delete the supplied parameters since we only copy their values in this case
	objarray->SetOwner(kTRUE);

	const UInt_t nPars = objarray->GetEntries();
	if ( nPars != slaveIndices_[iSlave].size() ) {
	std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Unexpected number of parameters received from slave" << std::endl;
	gSystem->Exit(1);
	}

	for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
	LauParameter* parameter = dynamic_cast<LauParameter>( (objarray)[iPar] );
	if ( ! parameter ) {
	std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Error reading parameter from slave" << std::endl;
	gSystem->Exit(1);
	}

	TString parname = parameter->name();
	Double_t parvalue = parameter->initValue();

	std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
	if ( iter == parIndices_.end() ) {
	std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Unexpected parameter name received from slave" << std::endl;
	gSystem->Exit(1);
	}

	const UInt_t index = iter->second;
	if ( slaveIndices_[iSlave][iPar] != index ) {
	std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Unexpected parameter received from slave" << std::endl;
	gSystem->Exit(1);
	}

	params_[index]->initValue( parvalue );
	parValues_[index] = parvalue;
	vectorPar_[iSlave][iPar] = parvalue;
	this->checkParameter( parameter, index );
	}

	delete objarray; objarray = 0;
	delete messageFromSlave_; messageFromSlave_ = 0;
	}
	}

	void LauSimFitMaster::getParametersFromSlavesFirstTime()
	{
	slaveIndices_.resize( nSlaves_ );
	vectorPar_.resize( nSlaves_ );
	vectorRes_.resize( nSlaves_ );

	TSocket* sActive(0);

	// Construct a message, requesting the list of parameter names
	TString msgStr = "Send Parameters";
	TMessage message( kMESS_STRING );
	message.WriteTString( msgStr );

	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
	// Send the message to the slave
	sSlaves_[iSlave]->Send(message);

	// Wait to receive the response and check that it has come from the slave we just requested from
	sActive = socketMonitor_->Select();
	if ( sActive != sSlaves_[iSlave] ) {
	std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Received message from a different slave than expected!" << std::endl;
	gSystem->Exit(1);
	}

	// Read the object and extract the parameter names
	sSlaves_[iSlave]->Recv( messageFromSlave_ );
	TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
	if ( ! objarray ) {
	std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameters from slave" << std::endl;
	gSystem->Exit(1);
	}

	const UInt_t nPars = objarray->GetEntries();

	vectorPar_[iSlave] = new Double_t[nPars];

	for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
	LauParameter* parameter = dynamic_cast<LauParameter>( (objarray)[iPar] );
	if ( ! parameter ) {
	std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameter from slave" << std::endl;
	gSystem->Exit(1);
	}

	TString parname = parameter->name();
	Double_t parvalue = parameter->initValue();

	std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
	if ( iter != parIndices_.end() ) {
	UInt_t index = iter->second;
	slaveIndices_[iSlave].push_back( index );
	this->checkParameter( parameter, index );
	} else {
	UInt_t index = parIndices_.size();
	parIndices_.insert( std::make_pair( parname, index ) );
	parNames_.insert( std::make_pair( index, parname ) );
	slaveIndices_[iSlave].push_back( index );
	params_.push_back( parameter );
	parValues_.push_back( parvalue );
	}
	vectorPar_[iSlave][iPar] = parvalue;
	}

	delete objarray; objarray = 0;
	delete messageFromSlave_; messageFromSlave_ = 0;
	}

	nParams_ = params_.size();
	}

	void LauSimFitMaster::printParInfo() const
	{
	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
	const std::vector<UInt_t>& indices = slaveIndices_[iSlave];

	std::cout << "INFO in LauSimFitMaster::printParInfo : Slave " << iSlave << " has the following parameters:\n";
	for ( std::vector<UInt_t>::const_iterator iter = indices.begin(); iter != indices.end(); ++iter ) {
	const TString& parName = parNames_.find(*iter)->second;
	Double_t parValue = parValues_[*iter];
	const LauParameter* par = params_[*iter];
	if ( par->name() != parName \|\| par->initValue() != parValue ) {
	std::cerr << "ERROR in LauSimFitMaster::printParInfo : Discrepancy in parameter name and value records, this is very strange!!" << std::endl;
	}

	std::cout << " : " << parName << " = " << parValue << " and has index " << *iter << "\n";
	}

	std::cout << std::endl;
	}

	std::cout << "INFO in LauSimFitMaster::printParInfo : " << "There are " << nParams_ << " parameters in total" << std::endl;
	}

	void LauSimFitMaster::checkParameter( const LauParameter* param, UInt_t index ) const
	{
	const LauParameter* storedPar = params_[index];

	TString parName = storedPar->name();
	if ( param->name() != parName ) {
	std::cerr << "ERROR in LauSimFitMaster::checkParameter : Parameter name is different!! This shouldn't happen!!" << std::endl;
	}
	if ( param->initValue() != storedPar->initValue() ) {
	std::cerr << "WARNING in LauSimFitMaster::checkParameter : Initial value for parameter " << parName << " is different, will use the value first set: " << storedPar->initValue() << std::endl;
	}
	if ( param->minValue() != storedPar->minValue() ) {
	std::cerr << "WARNING in LauSimFitMaster::checkParameter : Minimum allowed value for parameter " << parName << " is different, will use the value first set: " << storedPar->minValue() << std::endl;
	}
	if ( param->maxValue() != storedPar->maxValue() ) {
	std::cerr << "WARNING in LauSimFitMaster::checkParameter : Maximum allowed value for parameter " << parName << " is different, will use the value first set: " << storedPar->maxValue() << std::endl;
	}
	if ( param->fixed() != storedPar->fixed() ) {
	std::cerr << "WARNING in LauSimFitMaster::checkParameter : Fixed/floating property of parameter " << parName << " is different, will use the value first set: " << (storedPar->fixed() ? "fixed" : "floating") << std::endl;
	}
	if ( param->firstStage() != storedPar->firstStage() ) {
	std::cerr << "WARNING in LauSimFitMaster::checkParameter : First stage property of parameter " << parName << " is different, will use the value first set: " << (storedPar->firstStage() ? "true" : "false") << std::endl;
	}
	if ( param->secondStage() != storedPar->secondStage() ) {
	std::cerr << "WARNING in LauSimFitMaster::checkParameter : Second stage property of parameter " << parName << " is different, will use the value first set: " << (storedPar->secondStage() ? "true" : "false") << std::endl;
	}
	}

	void LauSimFitMaster::initialise()
	{
	this->initSockets();
	}

	void LauSimFitMaster::runSimFit( const TString& fitNtupleFileName, UInt_t nExpt, UInt_t firstExpt, Bool_t useAsymmErrors, Bool_t twoStageFit )
	{
	// Routine to perform the total fit.

	// First, initialise
	this->initialise();

	std::cout << "INFO in LauSimFitMaster::runSimFit : First experiment = " << firstExpt << std::endl;
	std::cout << "INFO in LauSimFitMaster::runSimFit : Number of experiments = " << nExpt << std::endl;

	// Start the cumulative timer
	cumulTimer_.Start();

	numberOKFits_ = 0, numberBadFits_ = 0;
	fitStatus_ = -1;

	// Create and setup the fit results ntuple
	std::cout << "INFO in LauSimFitMaster::runSimFit : Creating fit ntuple." << std::endl;
	if (fitNtuple_ != 0) {delete fitNtuple_; fitNtuple_ = 0;}
	fitNtuple_ = new LauFitNtuple(fitNtupleFileName, useAsymmErrors);

	// Loop over the number of experiments
	for (iExpt_ = firstExpt; iExpt_ < (firstExpt+nExpt); ++iExpt_) {

	// Start the timer to see how long each fit takes
	timer_.Start();

	// Instruct the slaves to read the data for this experiment
	Bool_t readOK = this->readData();
	if ( ! readOK ) {
	std::cerr << "ERROR in LauSimFitMaster::runSimFit : One or more slaves reported problems with reading data for experiment " << iExpt_ << ", skipping..." << std::endl;
	timer_.Stop();
	continue;
	}

	// Instruct the slaves to perform the caching
	this->cacheInputData();

	// Do the fit
	this->fitExpt( useAsymmErrors, twoStageFit );

	// Stop the timer and see how long the program took so far
	timer_.Stop();
	timer_.Print();

	// Instruct the slaves to finalise the results
	this->finalise();

	// Keep track of how many fits succeeded or failed
	if (fitStatus_ == 3) {
	++numberOKFits_;
	} else {
	++numberBadFits_;
	}
	}

	// Print out total timing info.
	std::cout << "INFO in LauSimFitMaster::runSimFit : Cumulative timing:" << std::endl;
	cumulTimer_.Stop();
	cumulTimer_.Print();

	// Print out stats on OK fits.
	std::cout << "INFO in LauSimFitMaster::runSimFit : Number of OK Fits = " << numberOKFits_ << std::endl;
	std::cout << "INFO in LauSimFitMaster::runSimFit : Number of Failed Fits = " << numberBadFits_ << std::endl;
	Double_t fitEff(0.0);
	if (nExpt != 0) {fitEff = numberOKFits_/(1.0*nExpt);}
	std::cout << "INFO in LauSimFitMaster::runSimFit : Fit efficiency = " << fitEff*100.0 << "%." << std::endl;

	// Instruct the slaves to write out any fit results (ntuples etc...).
	this->writeOutResults();
	}

	Bool_t LauSimFitMaster::readData()
	{
	if ( socketMonitor_ == 0 ) {
	std::cerr << "ERROR in LauSimFitMaster::readData : Sockets not initialised." << std::endl;
	return kFALSE;
	}

	// Construct a message, requesting to read the data for the given experiment
	TString msgStr("Read Expt");
	TMessage message( kMESS_STRING );
	message.WriteTString( msgStr );
	message.WriteUInt( iExpt_ );

	// Send the message to the slaves
	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
	sSlaves_[iSlave]->Send(message);
	}

	TSocket* sActive(0);
	UInt_t responsesReceived(0);
	Bool_t ok(kTRUE);
	while ( responsesReceived != nSlaves_ ) {

	// Get the next queued response
	sActive = socketMonitor_->Select();

	// Extract from the message the ID of the slave and the number of events read
	sActive->Recv( messageFromSlave_ );
	UInt_t iSlave(0);
	UInt_t nEvents(0);
	messageFromSlave_->ReadUInt( iSlave );
	messageFromSlave_->ReadUInt( nEvents );

	if ( nEvents <= 0 ) {
	std::cerr << "ERROR in LauSimFitMaster::readData : Slave " << iSlave << " reports no events found for experiment " << iExpt_ << std::endl;
	ok = kFALSE;
	} else {
	std::cerr << "INFO in LauSimFitMaster::readData : Slave " << iSlave << " reports " << nEvents << " events found for experiment " << iExpt_ << std::endl;
	}

	++responsesReceived;
	}

	return ok;
	}

	Bool_t LauSimFitMaster::cacheInputData()
	{
	if ( socketMonitor_ == 0 ) {
	std::cerr << "ERROR in LauSimFitMaster::cacheInputData : Sockets not initialised." << std::endl;
	return kFALSE;
	}

	// Construct a message, requesting it to read the data for the given experiment
	TString msgStr("Cache");
	TMessage message( kMESS_STRING );
	message.WriteTString( msgStr );

	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
	// Send the message to the slave
	sSlaves_[iSlave]->Send(message);
	}

	TSocket* sActive(0);
	UInt_t responsesReceived(0);
	Bool_t allOK(kTRUE);
	while ( responsesReceived != nSlaves_ ) {

	// Get the next queued response
	sActive = socketMonitor_->Select();

	// Extract from the message the ID of the slave and the success/failure flag
	sActive->Recv( messageFromSlave_ );
	UInt_t iSlave(0);
	Bool_t ok(kTRUE);
	messageFromSlave_->ReadUInt( iSlave );
	messageFromSlave_->ReadBool( ok );

	if ( ! ok ) {
	std::cerr << "ERROR in LauSimFitMaster::cacheInputData : Slave " << iSlave << " reports an error performing caching" << std::endl;
	allOK = kFALSE;
	}

	++responsesReceived;
	}

	return allOK;
	}

	void LauSimFitMaster::checkInitFitParams()
	{
	this->getParametersFromSlaves();
	this->printParInfo();
	}

	void LauSimFitMaster::fitExpt( Bool_t useAsymmErrors, Bool_t twoStageFit )
	{
	// Routine to perform the actual fit for the given experiment

	// Instruct the salves to update initial fit parameters if required (e.g. if using random numbers).
	this->checkInitFitParams();

	// Initialise the fitter
	LauFitter::fitter()->useAsymmFitErrors( useAsymmErrors );
	LauFitter::fitter()->twoStageFit( twoStageFit );
	LauFitter::fitter()->initialise( this, params_ );

	nParams_ = LauFitter::fitter()->nParameters();
	nFreeParams_ = LauFitter::fitter()->nFreeParameters();

	// Now ready for minimisation step
	std::cout << "\nINFO in LauSimFitMaster::fitExpt : Start minimisation...\n";
	std::pair<Int_t,Double_t> fitResult = LauFitter::fitter()->minimise();

	fitStatus_ = fitResult.first;
	NLL_ = fitResult.second;

	// If we're doing a two stage fit we can now release (i.e. float)
	// the 2nd stage parameters and re-fit
	if ( twoStageFit ) {

	if ( fitStatus_ != 3 ) {
	std::cerr << "ERROR in LauSimFitMaster:fitExpt : Not running second stage fit since first stage failed." << std::endl;
	LauFitter::fitter()->releaseSecondStageParameters();
	} else {
	LauFitter::fitter()->fixFirstStageParameters();
	LauFitter::fitter()->releaseSecondStageParameters();
	nParams_ = LauFitter::fitter()->nParameters();
	nFreeParams_ = LauFitter::fitter()->nFreeParameters();
	fitResult = LauFitter::fitter()->minimise();
	}
	}

	fitStatus_ = fitResult.first;
	NLL_ = fitResult.second;
	const TMatrixD& covMat = LauFitter::fitter()->covarianceMatrix();
	covMatrix_.Clear();
	covMatrix_.ResizeTo( covMat.GetNrows(), covMat.GetNcols() );
	covMatrix_.SetMatrixArray( covMat.GetMatrixArray() );

	// Store the final fit results and errors into protected internal vectors that
	// all sub-classes can use within their own finalFitResults implementation
	// used below (e.g. putting them into an ntuple in a root file)
	LauFitter::fitter()->updateParameters();
	LauFitter::fitter()->releaseFirstStageParameters();
	}

	void LauSimFitMaster::setParsFromMinuit(Double_t* par, Int_t npar)
	{
	// This function sets the internal parameters based on the values
	// that Minuit is using when trying to minimise the total likelihood function.

	- // MINOS reports one fewer free parameter than there nominally is
	- // (since it is doing a scan over one parameter at a time),
	- // so increment npar so the following check doesn't fail
	- if ( withinAsymErrorCalc_ ) {
	- ++npar;
	- }
	-
	- if (static_cast<UInt_t>(npar) != nFreeParams_) {
	- std::cerr << "ERROR in LauSimFitMaster::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n";
	- std::cerr << " Expected: " << nFreeParams_ << ".\n" << std::endl;
	- gSystem->Exit(EXIT_FAILURE);
	+ // MINOS reports different numbers of free parameters depending on the
	+ // situation, so disable this check
	+ if ( ! withinAsymErrorCalc_ ) {
	+ if (static_cast<UInt_t>(npar) != nFreeParams_) {
	+ std::cerr << "ERROR in LauSimFitMaster::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n";
	+ std::cerr << " Expected: " << nFreeParams_ << ".\n" << std::endl;
	+ gSystem->Exit(EXIT_FAILURE);
	+ }
	}

	// Despite npar being the number of free parameters
	// the par array actually contains all the parameters,
	// free and floating...
	// Update all the parameters with their new values.
	for (UInt_t i(0); i<nParams_; ++i) {
	parValues_[i] = par[i];
	}
	}

	Double_t LauSimFitMaster::getTotNegLogLikelihood()
	{
	if ( socketMonitor_ == 0 ) {
	std::cerr << "ERROR in LauSimFitMaster::getTotNegLogLikelihood : Sockets not initialised." << std::endl;
	return 0.0;
	}

	// Send current values of the parameters to the slaves.
	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {

	std::vector<UInt_t>& indices = slaveIndices_[iSlave];
	UInt_t nPars = indices.size();
	for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
	vectorPar_[iSlave][iPar] = parValues_[ indices[iPar] ];
	}

	TMessage* message = messagesToSlaves_[iSlave];
	message->Reset( kMESS_ANY );
	message->WriteUInt( nPars );
	message->WriteFastArray( vectorPar_[iSlave], nPars );

	sSlaves_[iSlave]->Send(*message);
	}

	Double_t negLogLike(0.0);
	TSocket *sActive(0);
	UInt_t responsesReceived(0);
	while ( responsesReceived != nSlaves_ ) {

	sActive = socketMonitor_->Select();
	sActive->Recv(messageFromSlave_);

	messageFromSlave_->ReadDouble( vectorRes_[responsesReceived] );

	negLogLike += vectorRes_[responsesReceived];

	++responsesReceived;
	}

	return negLogLike;
	}

	Bool_t LauSimFitMaster::finalise()
	{
	if ( socketMonitor_ == 0 ) {
	std::cerr << "ERROR in LauSimFitMaster::finalise : Sockets not initialised." << std::endl;
	return kFALSE;
	}

	// Prepare the covariance matrices
	covMatrices_.resize( nSlaves_ );

	LauParamFixed pred;

	std::map<UInt_t,UInt_t> freeParIndices;

	UInt_t counter(0);
	for ( UInt_t iPar(0); iPar < nParams_; ++iPar ) {
	const LauParameter* par = params_[iPar];
	if ( ! pred(par) ) {
	freeParIndices.insert( std::make_pair(iPar,counter) );
	++counter;
	}
	}

	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
	const UInt_t nPar = slaveIndices_[iSlave].size();

	std::vector<UInt_t> freeIndices;
	freeIndices.reserve( nPar );

	for ( UInt_t iPar(0); iPar < nPar; ++iPar ) {
	UInt_t index = slaveIndices_[iSlave][iPar];
	std::map<UInt_t,UInt_t>::iterator freeIter = freeParIndices.find(index);
	if ( freeIter == freeParIndices.end() ) {
	continue;
	}
	UInt_t freeIndex = freeIter->second;
	freeIndices.push_back( freeIndex );
	}

	const UInt_t nFreePars = freeIndices.size();
	TMatrixD& covMat = covMatrices_[iSlave];
	covMat.ResizeTo( nFreePars, nFreePars );

	for ( UInt_t iPar(0); iPar < nFreePars; ++iPar ) {
	for ( UInt_t jPar(0); jPar < nFreePars; ++jPar ) {
	UInt_t i = freeIndices[iPar];
	UInt_t j = freeIndices[jPar];
	covMat( iPar, jPar ) = covMatrix_( i, j );
	}
	}
	}

	// The array to hold the parameters
	TObjArray array;

	// Send messages to all slaves containing the final parameters and fit status, NLL
	// TODO - at present we lose the information on the correlations between the parameters that are unique to each slave
	// - so should we store the full correlation matrix in an ntuple? along with all the parameters?
	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {

	array.Clear();

	std::vector<UInt_t>& indices = slaveIndices_[iSlave];
	UInt_t nPars = indices.size();
	for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
	array.Add( params_[ indices[iPar] ] );
	}

	TMatrixD& covMat = covMatrices_[iSlave];

	TMessage* message = messagesToSlaves_[iSlave];
	message->Reset( kMESS_OBJECT );
	message->WriteInt( fitStatus_ );
	message->WriteDouble( NLL_ );
	message->WriteObject( &array );
	message->WriteObject( &covMat );

	sSlaves_[iSlave]->Send(*message);
	}

	TSocket *sActive(0);
	UInt_t responsesReceived(0);
	Bool_t allOK(kTRUE);
	while ( responsesReceived != nSlaves_ ) {

	// Get the next queued response
	sActive = socketMonitor_->Select();

	// Extract from the message the ID of the slave and the number of events read
	sActive->Recv( messageFromSlave_ );
	UInt_t iSlave(0);
	Bool_t ok(kTRUE);
	messageFromSlave_->ReadUInt( iSlave );
	messageFromSlave_->ReadBool( ok );

	if ( ok ) {
	TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
	if ( ! objarray ) {
	std::cerr << "ERROR in LauSimFitMaster::finalise : Error reading finalised parameters from slave" << std::endl;
	allOK = kFALSE;
	} else {
	// We want to auto-delete the supplied parameters since we only copy their values in this case
	objarray->SetOwner(kTRUE);

	const UInt_t nPars = objarray->GetEntries();
	if ( nPars != slaveIndices_[iSlave].size() ) {
	std::cerr << "ERROR in LauSimFitMaster::finalise : Unexpected number of finalised parameters received from slave" << std::endl;
	allOK = kFALSE;
	} else {
	for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
	LauParameter* parameter = dynamic_cast<LauParameter>( (objarray)[iPar] );
	if ( ! parameter ) {
	std::cerr << "ERROR in LauSimFitMaster::finalise : Error reading parameter from slave" << std::endl;
	allOK = kFALSE;
	continue;
	}

	TString parname = parameter->name();

	std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
	if ( iter == parIndices_.end() ) {
	std::cerr << "ERROR in LauSimFitMaster::finalise : Unexpected parameter name received from slave" << std::endl;
	allOK = kFALSE;
	continue;
	}

	const UInt_t index = iter->second;
	if ( slaveIndices_[iSlave][iPar] != index ) {
	std::cerr << "ERROR in LauSimFitMaster::finalise : Unexpected parameter received from slave" << std::endl;
	allOK = kFALSE;
	continue;
	}

	Double_t parvalue = parameter->value();
	params_[index]->value( parvalue );
	parValues_[index] = parvalue;
	vectorPar_[iSlave][iPar] = parvalue;
	}
	}
	delete objarray;
	}
	} else {
	std::cerr << "ERROR in LauSimFitMaster::finalise : Slave " << iSlave << " reports an error performing finalisation" << std::endl;
	allOK = kFALSE;
	}

	++responsesReceived;
	}

	// Fill our ntuple as well
	if ( fitNtuple_ != 0 ) {
	for ( std::vector<LauParameter*>::iterator iter = params_.begin(); iter != params_.end(); ++iter ) {
	if (!(*iter)->fixed()) {
	(*iter)->updatePull();
	}
	}
	std::vector<LauParameter> extraVars;
	fitNtuple_->storeParsAndErrors( params_, extraVars );
	fitNtuple_->storeCorrMatrix( iExpt_, NLL_, fitStatus_, covMatrix_ );
	fitNtuple_->updateFitNtuple();
	}

	return allOK;
	}

	Bool_t LauSimFitMaster::writeOutResults()
	{
	if ( socketMonitor_ == 0 ) {
	std::cerr << "ERROR in LauSimFitMaster::writeOutResults : Sockets not initialised." << std::endl;
	return kFALSE;
	}

	// Construct a message, requesting to write out the fit results
	TString msgStr("Write Results");
	TMessage message( kMESS_STRING );
	message.WriteTString( msgStr );

	// Send the message to the slaves
	for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
	sSlaves_[iSlave]->Send(message);
	}

	TSocket *sActive(0);
	UInt_t responsesReceived(0);
	Bool_t allOK(kTRUE);
	while ( responsesReceived != nSlaves_ ) {

	// Get the next queued response
	sActive = socketMonitor_->Select();

	// Extract from the message the ID of the slave and the number of events read
	sActive->Recv( messageFromSlave_ );
	UInt_t iSlave(0);
	Bool_t ok(kTRUE);
	messageFromSlave_->ReadUInt( iSlave );
	messageFromSlave_->ReadBool( ok );

	if ( ! ok ) {
	std::cerr << "ERROR in LauSimFitMaster::writeOutResults : Slave " << iSlave << " reports an error performing finalisation" << std::endl;
	allOK = kFALSE;
	}

	++responsesReceived;
	}

	// Write out our ntuple as well
	if (fitNtuple_ != 0) {
	fitNtuple_->writeOutFitResults();
	}

	return allOK;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:17 PM (8 h, 6 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486514
Default Alt Text: (68 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions